Programação multi-thread em python

Esse talvez possa parecer um post meio maluco, para o que comumente discutimos aqui nesse blog, mas é muito interessante e curioso, para saber o que está acontecendo sobre o capo do carro, ou do computador.

Uma das coisas mais diferentes que já vi, foi quando fiz a disciplina de sistemas operacionais, isso porque tudo que eu achava que sabia sobre computador, vi que estava errado. Nem tudo é tão simples como um script continuo que fazemos aqui no R, e uma prova disso é a programação Multi-thread.

thread

Vamos pensar na figura acima, tomando como exemplo os MCMC que vemos comumente aqui, veja alguma exemplos de implementação aqui de um gibbs samples e aqui do mais simples metropolis hastings.

Nesses dois casos fizemos algo como em single-threaded, isso porque fizemos uma cadeia e fomos realizando vários cálculos e salvando resultados. Acontece que, todas as analises de dados que fazemos, não usamos apenas uma cadeia, usamos várias cadeias, no último post mesmo, usando o openbugs, usamos três cadeias, que podem ser processadas separadamente, podem ser multithreaded.
Os valores das amostras, parâmetros, modelos, tudo é igual, então não compensa gravar isso varias vezes na memoria, é melhor compartilhar, mas para uma cadeia de markov, a gente sempre precisa do valor anterior para continuar ela, então as cadeias não precisam dividir os cálculos que estão fazendo entre si, apenas os dados e parâmetros e modelos, agora olha no modelinho ali em cima, não é exatamente isso que diz respeito a programação multi-thread?
Veja que com multi-thread, cada thread tem sua pilha, seus registradores, seus dados particulares, que é sua cadeia de markov que está computando, mas também divide dados, que são as amostras, modelo, parâmetros entre todas as threads.

E qual a relevância disso? Ai que entre aqueles processadores de muitos núcleos, sem isso, sem paralelizar, dividir a computação, a gente teria que fazer uma cadeia no processador, terminar, fazer outra, terminar, pense que leva um minuto por cadeia, com três então levamos três minutos para fazer tudo, mas se fizermos algo multi-thread, da para mandar cada cadeia para cada núcleo do processador (mais ou menos isso, simplificando, sei que não é assim), como cada um leva 1 minuto e da para fazer em paralelo, ao invés de esperar 3 minutos pelo resultado, esperamos um.

image0021225707175745

Podemos fazer um exemplo em python que é bem simples, usando a biblioteca threading.

Para esse exemplo simples, vamos criar uma classe bem bobinha

1
2
3
4
5
6
7
8
9
10
class minhaThread (threading.Thread):
    def __init__(self, threadID, nome, contador):
        threading.Thread.__init__(self)
        self.threadID = threadID
        self.nome = nome
        self.contador = contador
    def run(self):
        print "Iniciando thread %s com %d processos" % (self.name,self.contador)
        processo(self.nome, self.contador)
        print "Finalizando " + self.nome

A gente ja falou de orientação a objetos em R aqui e aqui, uma das coisas que a orientação a objetos permite é a herança, que é herdar tudo que a classe pai tem, aqui estamos herdando a classe Thread definida na biblioteca threading, isso é o que esta escrito aqui:

1
class minhaThread (threading.Thread):

Por isso o parenteses, lembro que ja usei classe em algum problema do rosalind aqui, mas não lembro em qual post agora, mas beleza, o resto é bem besta a classe vai ter 3 atributos, um ID, um nome e um contador, que seria algo como os parâmetros para fazer um mcmc, vamos abstrair de uma atividade útil agora, vamos apenas pensar em termos de computação.

1
2
3
4
5
    def __init__(self, threadID, nome, contador):
        threading.Thread.__init__(self)
        self.threadID = threadID
        self.nome = nome
        self.contador = contador

Então nos construtor da thread definimos esses atributos e ela vai ter apenas um método, uma função chamada run.

1
2
3
4
    def run(self):
        print "Iniciando thread %s com %d processos" % (self.name,self.contador)
        processo(self.nome, self.contador)
        print "Finalizando " + self.nome

Tem que chamar run, por causa que é assim que funciona o pacote, mas o run vai so imprimir que estamos começando, processar algo, pensando que processar pode ser qualquer coisa e depois avisar que terminamos.

Aqui o processo definimos assim

1
2
3
4
def processo(nome, contador):
    while contador:
        print "Thread %s fazendo o processo %d" % (nome, contador)
        contador -= 1

Ele é uma função externa, não é da classe, mas ele apenas imprimi uma mensagem, falando que esta processando, e diminui o contador em 1, até zerar, e terminar.

Depois disso está tudo pronto, so precisamos criar nossas threads

1
2
3
# Criando as threads
thread1 = minhaThread(1, "Alice", 8)
thread2 = minhaThread(2, "Bob", 8)

Criamos elas, simples assim, aqui no caso são duas, chamadas de alice e bob, hehe, alice e bob são exemplos mais comuns em física, criptografia e teoria dos jogos, ao invés de falar A e B, os caras falam Alice e Bob, para ficar mais fácil de entender.

Ok, mas apos criado, nada foi processado ainda, dai então começamos elas.

1
2
3
# Comecando novas Threads
thread1.start()
thread2.start()

O resto do programa é meramente para fazer a thread mãe esperar alice e bob para terminar, não convêm explicar agora, mas essa parte so faz isso, esperar alice e bob terminarem

1
2
3
4
5
6
threads = []
threads.append(thread1)
threads.append(thread2)
 
for t in threads:
    t.join()

Agora é legal quando a gente roda esse programa.

As vezes ele roda uma thread seguida da outra

Screenshot from 2015-07-23 16:21:42

E as vezes tudo vira uma loucura.

Screenshot from 2015-07-23 16:23:00

Isso porque uma vez que você inicio o processamento, o sistema operacional que vai escalonar as threads, dai a gente não manda mais, só o sistema operacional, então se o processador ta liberado, ele manda elas para la para fazer conta, mas as vezes ta ocupado, então a thread tem que esperar sua vez, agora uma thread pode ir para um núcleo e ir até o final la, enquanto a outra entrou num outro e teve que esperar alguém terminar de processar, como no segundo exemplo, Alice começou, o Bob esperou um pouco, Alice ja tinha processado 2 vezes, ai o Bob começou a trabalhar, mas terminou depois.

Um exemplo bem besta, mas não é preciso saber fazer muitas contas complexas para imaginar que isso é melhor que primeiro alice fazer todo o trabalho dela e depois bob fazer todo o trabalho dela. Pena que isso é extremamente complexo, ainda mais quando a gente tem que compartilhar algum dado e mexer nele, porque ai a vez de quem mexeu pode ser importante, mas isso é outra historia.

Claro que programação não é a meta de muitas pessoas que visitam aqui, mas ter uma ideia que isso existe, pode ajudar a entender comentários de manuais de programa entre outras coisas.

Por exemplo, um programa comum em biologia é o mrBayes, basta olhar o manual que você vai notar alguma notas de instalação, veja essa parte aqui:

instructions for compiling and running the parallel version of the program

Ou seja, se você instala as coisas so dando dois cliques, não le o manual, pode estar perdendo essa possibilidade de usar multi-thread e algo que pode rodar em uma hora vai levar um dia. Ja que ainda no exemplo do mrBayes, da para usar até o processador da placa de video para ajudar a fazer cálculos e agilizar o processamento. E pode abrir o olho, que tudo que você ler MCMC é um forte candidato a usar multi-thread.

O R em si não é multi-thread, mas existem vários pacotes que tentam incorporar isso, como visto no taskview de high performance. Além de que o R possui interfaces que ligam eles a muitos outros programas, como é o caso do Openbugs mesmos, que os calculos são feitos neles, mas depois devolvemos os dados no R para fazer figuras e outras analises, mas temos muitos outros exemplo.

Apesar que na maioria das vezes isso não é muito relevante, para fazer regressões e muitos modelos, tudo é tão rápido que a gente nem sente diferença, mas algumas coisas mais avançadas (talvez “avançando” seja um termo ruim aqui, coisas que precisam fazer mais contas, comum em evolução), isso começa a ser importante.

Bem é isso ai, a primeira vez que ouvi sobre esse negocio de multi-thread, a primeira coisa que veio na minha cara é MCMC, acho que ter ouvido falar de MCMC antes de fazer Sistemas Operacionais foi um diferencial para tornar tudo mais interessante e belo, mas tive um excelente professor em sistemas operacionais também para ajudar. O script vai estar la no repositório recologia, tem mais alguma exemplos da minha aula de sistemas operacionais no github em outro repositório aqui mas em linguagem C, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e até mais.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
import threading
 
class minhaThread (threading.Thread):
    def __init__(self, threadID, nome, contador):
        threading.Thread.__init__(self)
        self.threadID = threadID
        self.nome = nome
        self.contador = contador
    def run(self):
        print "Iniciando thread %s com %d processos" % (self.name,self.contador)
        processo(self.nome, self.contador)
        print "Finalizando " + self.nome
 
def processo(nome, contador):
    while contador:
        print "Thread %s fazendo o processo %d" % (nome, contador)
        contador -= 1
 
# Criando as threads
thread1 = minhaThread(1, "Alice", 8)
thread2 = minhaThread(2, "Bob", 8)
 
# Comecando novas Threads
thread1.start()
thread2.start()
 
threads = []
threads.append(thread1)
threads.append(thread2)
 
for t in threads:
    t.join()
 
print "Saindo da main"

Mergesort em C++ usando o pacote Rcpp

A muito tempo atrás eu comecei a ver como dava para usar código em linguagem C do R neste post aqui. Depois disso, no comentário alguém mais esperto que eu falou sobre um tal de pacote Rcpp, que permitia de maneira mais simples, usar código em linguagem C++ do R, depois disso eu fiz mais alguns testes, tentando implementar alguns algorítimos de ordenação, o bubblesort aqui e o insertion sort aqui.

Esse é o tipo de coisa inútil de se fazer, mas é divertido. E esses algorítimos de ordenação tem em qualquer curso introdutório de programação, então são implementações básicas para qualquer linguagem. Depois do mergesort eu ainda tenho o quicksort, mas no mergesort tudo começa a ficar mais complexo.

Por curiosidade, o método padrão na função sort() do R é o Shellsort, criado por Donald Shell mas o que está implementado no R é uma variante do Robert Sedgewick, esse cara da quatro cursos sobre programação gratuitos no Coursera, so olhar aqui.

Mas na função sort do R também temos o Quicksort, mas ele vai ficar pra próxima.

Agora vamos ver o mergesort. Ele é um algorítimo recursivo, ou seja uma função que chama ela mesma, e a ideia básica é, melhor que ordenar todo um vetor, é pegar dois vetores ordenados e então simplesmente intercalar eles. Mas a gente precisa de dois vetores ordenados pra poder usar esse esquema certo?
Pra isso é so a gente quebrar o vetor no meio, e de novo no meio para as duas metades, e de novo para os quatro quartos, e de novo até que vai chegar um momento que so vamos ter mini vetores de um elemento, um vetor de um elemento está necessariamente ordenado, logo agora é so voltar intercalando esses vetores que fica tudo bem.

MergeSort

É até estranho como isso funciona, mas funciona. Bem o legal que diferente do InsertionSort e o BubbleSort, aqui a gente vai precisar de uma função auxiliar, então a gente vai escrever o código em um arquivo, que vamos chamar de merge.cpp, e depois vamos chamar esse arquivo do R usando a função sourceCpp() do pacote Rcpp.

Então o código do mergesort é assim:

#include <Rcpp.h>
using namespace Rcpp;
 
void intercala(int p, int q, int r, NumericVector v,  NumericVector w)
{
  int i, j, k;
   i = p;
   j = q;
   k = 0;
   while (i < q && j < r) {
      if (v[i] < v[j]) {
         w[k] = v[i];
         i++;
      }
      else {
         w[k] = v[j];
         j++;
      }
      k++;
   }
   while (i < q) {
      w[k] = v[i];
      i++;
      k++;
   }
   while (j < r) {
      w[k] = v[j];
      j++;
      k++;
   }
   for (i = p; i < r; i++)
      v[i] = w[i-p];
}
 
void mergesort(int p, int r, NumericVector v, NumericVector aux)
{
   int q;
   if (p < r - 1) {
      q = (p + r) / 2;
      mergesort(p, q, v,aux);
      mergesort(q, r, v,aux);
      intercala(p, q, r, v,aux);
   }
}
 
// [[Rcpp::export]]
NumericVector mergesortC(NumericVector vetor) {
  Rcpp::NumericVector saida = Rcpp::clone(vetor);
  Rcpp::NumericVector aux = Rcpp::clone(vetor);
  int n = saida.size();
  mergesort(0,n,saida,aux);
  return saida;
}

Então, tudo isso está dentro de um arquivo que eu chamei de merge.cpp
Basicamente o código você acha em qualquer esquina da rede mundial de computadores, mas especificamente aqui, temos que importar a biblioteca Rcpp.h, já que estamos usando um objeto específico dela, que é o NumericVector. Outra coisa é que o R manda um ponteiro pra cá, que aponta aonde estão os dados, desse objeto. Se a gente tentar mexer direto nele, não vai dar certo, eu tentei e fiz o R travar com isso heheh, ai perguntando no forum aqui, o cara que criou o pacote me mandou ir estudar, ai lendo o livro dele eu vi que mexer direto nos dados, muitas vezes não é uma boa ideia, dai a gente usa o método clone de NumericVector, a gente clona ele para outro vetor, e então manipula esse novo vetor, pra não mexer na memoria que não devemos.
Outra coisa é que, talvez um ponto fraco do mergeSort, é que para intercalar dois vetores, precisamos de outro lugar para ir guardando eles em ordem, ou seja, diferente do insertionsort, não da para ir arrumando dentro do próprio vetor, por isso precisamos de um vetor auxiliar, aqui eu ja criei um vetor auxiliar chamado aux, para ser usado para esse propósito. Eu até tentei criar vetores dentro da função intercala, mas isso deixou tudo bem lento, desse jeito a velocidade ficou um pouco melhor. Por último, mas bem legal, é que nesse caso, veja que temos três funções, uma que recebe e retorna tudo para o R, e outras duas que fazem o trabalho propriamente dito. Mas quando a gente olha no R, desse jeito a gente só vai ver a mergesortC, isso porque em cima dela tem escrito

// [[Rcpp::export]]

Isso faz com que a gente só veja ela, e não as funções auxiliares.

Concluído a construção da função em C++, agora é simples no R

> library(Rcpp) > sourceCpp("merge.cpp") > mergesortC(sample(1:10)) [1] 1 2 3 4 5 6 7 8 9 10

Basicamente, só damos um sourceCpp no arquivo com o código, podemos usar o argumento verbose=TRUE se quisermos acompanhar a compilação e ta tudo pronto para usar, e podemos ver que o mergesortC funciona certinho.

Agora a gente pode até medir o tempo que ele precisa para organizar, por exemplo, 1000 números inteiros, usando a função system.time.

> system.time(mergesortC(sample(1:1000))) usuário sistema decorrido 0.000 0.000 0.001

Alias, para ter uma noção de porque as pessoas continuaram inventando algorítimos de ordenação, a gente pode fazer uma comparação, entre o mergesort e o insertionsort por exemplo, e ver qual deles é o mais rápido, para organizar um vetor de um tamanho qualquer, é so a gente ficar bagunçando um vetor, com sample, e registar o tempo que ele leva para organizar, e repetir isso varias vezes.

dados<-data.frame(Tempo=rep(NA,1000),Algoritimo=rep(c("InsertionSort","MergeSort"),each=500)) for(i in 1:1000){ if(i<=500){ dados[i,1]<-system.time(insertionsortC(sample(1:10000)))[3] }else{ dados[i,1]<-system.time(mergesortC(sample(1:10000)))[3] } print(i) }

A gente cria então um dataframe para receber os dados, e faz o experimento, pegando o tempo total gasto no processo, o insertionsortC eu peguei do post passado.

O resultado disso é o seguinte.

01

É, o mergesort pode ser mais complicado de entender, e implementar, mas é muito mais rápido.

> aggregate(dados$Tempo, list(dados$Algoritimo), mean) Group.1 x 1 InsertionSort 0.043494 2 MergeSort 0.005362 > aggregate(dados$Tempo, list(dados$Algoritimo), sd) Group.1 x 1 InsertionSort 0.0044236978 2 MergeSort 0.0004933993

Ele é muitas vezes mais rápido, e varia muito menos no tempo de execução, não é a toa que o sort do R não usa o InsertionSort, mas ele também não usa o mergesort, ou seja, tem como melhorar ainda, mas essa vai ficar para a próxima porque ja está na hora de dormir.

Bem o script está aqui embaixo além do repositório recologia no github. Esse foi o primeiro código no Rcpp que envolvia mais de uma função, e no final muitas coisas envolvem varias funções, então acho que assim vamos abrindo portas para possibilidades legais. Outra coisa é que esse tipo de estratégia para resolver problemas, de partir o problema em partes menores, e resolver essas partes antes de resolver um problema maior, pode parecer algo longe para biologia, mas a gente vai pela mesma estrada em problemas como o de alinhamento de sequências. E vamos dormir.

library(Rcpp)
 
sourceCpp("merge.cpp",verbose=T)
 
mergesortC(sample(1:10))
 
system.time(mergesortC(sample(1:1000)))
 
cppFunction("
    NumericVector insertionsortC(NumericVector vetor) {
        int n = vetor.size();
 
        double aux;
        int i , j;
 
        for(i=1;i<n;i++) {
            aux=vetor[i];
            j=i-1;
            while(j>=0 && vetor[j]>aux) {
                vetor[j+1]=vetor[j];
                j=j-1;
                }
            vetor[j+1]=aux;
            }
        return vetor;
        }
")
 
dados<-data.frame(Tempo=rep(NA,1000),Algoritimo=rep(c("InsertionSort","MergeSort"),each=500))
 
for(i in 1:1000){
    if(i<=500){
        dados[i,1]<-system.time(insertionsortC(sample(1:10000)))[3]
    }else{
        dados[i,1]<-system.time(mergesortC(sample(1:10000)))[3]
    }
    print(i)
}
 
#Figura 1
boxplot(Tempo~Algoritimo,data=dados)
 
aggregate(dados$Tempo, list(dados$Algoritimo), mean)
aggregate(dados$Tempo, list(dados$Algoritimo), sd)

Optimização: O método de Newton

A gente viu a alguns posts atrás, como achar a raiz de uma função, ou raízes, usando o método de Newton-Raphson.

Como a gente viu naquele post, ele é a base para um método muito semelhante, mas que visa optimizar um valor em relação a uma função, sendo optimizar no sentido de encontrar o minimo ou máximo. Existem muitos algoritimos para isso, a gente inclusive ja viu um, chamado gradient descent, mas vamos ver aqui porque o método de Newton-Raphson pode servir para esse propósito com uma pequena modificação.

Vamos ver aqui como encontrar o máximo de uma função, mas encontrar o minimo é so multiplicar a função por -1, que ela inverte e pronto, se a gente achar o máximo, a gente acha o mínimo dela. E vice versa.

Suponha que temos uma função f:\mathbb{R} \to \mathbb{R} com a primeira e segunda derivativa contínua. f então tem um máximo global em x_i se f(x) \leq f(x_i) para todo f(x). A função f tem uma máxima local em x_i se f(x) \leq f(x_i) para todo x dentro de um certo intervalo \epsilon, (ou melhor, todo x tal que \vert x - x_i \vert \leq \epsilon), trocando em miúdos, quando a gente olha a função, se a gente ver ela oscilando, ela vai formar montanhas, então máxima global é o maior pico de montanha da função, e maxima local é o maior pico de montanha do intervalo da função que a gente ta olhando no gráfico, e particularmente, a gente olha um bocado a função de verossimilhança, e a gente sempre procura a máxima verosimilhança!

Uma condição necessária para x_i ser uma máxima local, é que f\prime(x) = 0 e f\prime\prime(x) \leq 0. Mas é suficiente que f\prime(x) = 0 e f\prime\prime(x) \le 0

Encontrar a máxima local é muito mais fácil que encontrar a máxima global, já que a gente para muitas funções, elas vão do menos infinito ao mais infinito, e a gente nunca vai conseguir percorrer toda ela. Mas dizer que encontramos um máximo local é garantido, dado a precisão que queremos.

Então, suponha que x_i é o máximo local da função f, Supondo que x_n \to x_i se n \to \infty, nos precisamos de um critério de parada para o momento em que \vert x_n - x_i \vert \leq \epsilon para uma tolerância predeterminada \epsilon. Infelizmente, no geral, isso não é tão simples, então a gente usa combinações dos seguintes critérios, na maioria dos algorítimos.

\vert x_n - x_{n-1} \vert \leq \epsilon
\vert f(x_n) - f(x_{n-1}) \vert \leq \epsilon
\vert f\prime(x_n) \vert \leq \epsilon

Se a sequência \{ x_n \}_n^\infty converge para uma máxima local, então todos os três critérios tem que ser satisfeitos. Mas o inverso não é verdadeiro. Logo mesmo que uma técnica para encontrar o máximo local convirja,nos ainda temos que checar que no final realmente se chegamos a um máximo local.

Outro problema é que \{ x_n \}_n^\infty pode simplesmente não convergir, a função, como o crescimento exponencial pode crescer para sempre, então temos que sempre definir um limite de iterações, senão a gente pode acabar presos em um loop infinito.

Agora, com nosso problema definido, como diabos Newton otimizava uma função?

Bem, se temos que f:[a,b] \to \mathbb{R} tem uma derivativa f\prime contínua, então o problema de achar o máximo local é equivalente a achar o máximo de f(a), f(b) e f(x_1), f(x_2), \dots , f(x_n) , onde x_1, x_2, \dots, x_n são as raízes de f\prime. Se nos aplicarmos o método de Newton-Raphson para achar a raiz de f\prime, nos então temos o método de Newton para otimizar a função f.

{x_{n+1}}= {{x_n} - {f\prime(x_n)\over{f\prime\prime(x_n)}}}

Por algum motivo, o Newton divide o crédito para o algorítimo de encontrar a raiz de uma função, mas não quando é aplicado a optimização. Agora para implementar, veja que teremos praticamente o mesmo algorítimo anterior, precisamos de uma função que devolva a valor de x para a função, mais a primeira e segunda derivada. O resto é mais ou menos igual. O algorítimo então no R fica assim.

newton <- function(ftn, x0, tol = 1e-9, n.max = 100) {
    x <- x0
    f.x <- ftn(x)
    n <- 0
    while ((abs(f.x[2]) > tol) & (n < n.max)) {
        x <- x - f.x[2]/f.x[3]
        f.x <- ftn(x)
        n <- n + 1
        cat("Na iteração", n, "o valor de x foi de", x, "\n")
    }
    if (n == n.max) {
        cat('Método de newton falhou em convergir\n')
    } else {
        cat("O algoritimo convergiu\n")
        return(x)
    }
}

Agora vamos pegar uma função de exemplo.
Vamos pegar a função gamma, da distribuição gamma, como esse é um método para duas dimensões, deixamos definidos os dois parâmetros dela como 2 e 3. de um ?gamma se quiser curiosar nessa função.

funcao <- function(x) {
    if (x < 0) return(c(0, 0, 0))
    if (x == 0) return(c(0, 0, NaN))
    y <- exp(-2*x)
    return(c(4*x^2*y, 8*x*(1-x)*y, 8*(1-2*x^2)*y))
}

Graficamente, ela tem essa forma.

01

Bem, a gente pode não saber o número, mas já sabemos onde é o máximo dela. É naquele topo ali.
O problema que antes do zero ela tem aquela área reta, e é só pensar, aquela área, vai dar um monte de zero dividido por zero, e o algorítimo vai estacionar ali e dali não vai sair, mesmo não sendo o máximo da função. Vem que quando a primeira derivada chegar a zero, vai ser um máximo, a não ser nesses casos especiais, por isso o ponto onde iniciamos o algorítimo pode ser importante, então varias tentativas, podem ser importantes, se você inicia em vários ponto diferentes e chega no mesmo lugar, temos mais certeza.

> newton(funcao, 0.25) Na iteração 1 o valor de x foi de 0.03571429 Na iteração 2 o valor de x foi de 0.001187431 Na iteração 3 o valor de x foi de 1.406649e-06 Na iteração 4 o valor de x foi de 1.978656e-12 O algoritimo convergiu [1] 1.978656e-12

Olha ai quando começamos perto daquela linha reta, a gente estaciona la, por causa desse espaço “especial” da função.

> newton(funcao, 0.5) Na iteração 1 o valor de x foi de 0 O algoritimo convergiu [1] 0

Mesmo começando em 0.5, a gente cai naquele mesmo lugar. Veja como o fato de a máxima e a minima ter a derivada igual a zero é um problema.

> newton(funcao, 0.75) Na iteração 1 o valor de x foi de 2.25 Na iteração 2 o valor de x foi de 1.941781 Na iteração 3 o valor de x foi de 1.662202 Na iteração 4 o valor de x foi de 1.418995 Na iteração 5 o valor de x foi de 1.222585 Na iteração 6 o valor de x foi de 1.085797 Na iteração 7 o valor de x foi de 1.017193 Na iteração 8 o valor de x foi de 1.000839 Na iteração 9 o valor de x foi de 1.000002 Na iteração 10 o valor de x foi de 1 O algoritimo convergiu [1] 1

Agora se a gente começa em 0.75, a gente consegue achar a máxima certinho.
Vamos olhar isso no gráfico para ver se esses números estão funcionando mesmo.

02

Bem vamos pegar uma função mais simples, vamos pegar aquela mesma função de segundo grau que vimos no método de Newton-Raphson. So vamos multiplicar ela por -1 para inverter ela.

funcao <- function(x) {
    fx <- -1*(2*x^2+4*x-4)
    dfx <- -1*(4*x+4)
    ddfx<- -1*(4)
    return(c(fx, dfx, ddfx))
}

Aqui, só existe um ponto na função onde a derivada é zero, então não tem erro.

> newton(funcao, -2) Na iteração 1 o valor de x foi de -1 O algoritimo convergiu [1] -1

E se você olhar esse ponto no gráfico, para conferir.

03

Aqui, veja que de qualquer lugar que a gente começa, a gente cai no mesmo lugar.

> newton(funcao, 0) Na iteração 1 o valor de x foi de -1 O algoritimo convergiu [1] -1 > newton(funcao, 1) Na iteração 1 o valor de x foi de -1 O algoritimo convergiu [1] -1

Agora vamos ver o que acontece quando a gente tenta esse método, com essa mesma função, sem inverter ela.

funcao <- function(x) {
    fx <- 2*x^2+4*x-4
    dfx <- 4*x+4
    ddfx<- 4
    return(c(fx, dfx, ddfx))
}

Quando a gente usa o nosso método de Newton

> newton(funcao, -2) Na iteração 1 o valor de x foi de -1 O algoritimo convergiu [1] -1

A gente cai no mesmo ponto. Isso porque essa função só tem um ponto onde a derivada é zero, lembre-se que a primeira derivada dela é uma reta, a segunda derivada é uma constante, então a gente fica no mesmo ponto, até porque para os dois outros lados, vamos para o infinito, ou seja o algorítimo não teria como parar, nunca.

04

E de todo lugar que a gente começa, cai no mesmo lugar.

> newton(funcao, 0) Na iteração 1 o valor de x foi de -1 O algoritimo convergiu [1] -1 > newton(funcao, 1) Na iteração 1 o valor de x foi de -1 O algoritimo convergiu [1] -1

Bem, da para ver, que quem não entende muito de matemática, o jeito é garantir sempre olhando graficamente, vendo figurinhas é fácil garantir que as coisas estão funcionando. Mas quando a gente não pega uns exemplos muito complexos, o algorítimo funciona bem. Se você pensar em como são os likelihoods, eles vão ser parecido talvez com a função de segundo grau invertido, então neles o algorítimo vai funcionar bem.

E é isso ai, o próximo passo, é adaptar esse algorítimo para superfícies de verossimilhança, quando a gente vai ter vários parâmetros. Scripts no repositório recologia e aqui em baixo, e até o próximo post.

Referencia:
Owen Jones, Robert Maillardet & Andrew Robinson 2009 – Introduction to Scientific Programming and Simulation Using R. CRC Press 499pp

newton <- function(ftn, x0, tol = 1e-9, n.max = 100) {
    x <- x0
    f.x <- ftn(x)
    n <- 0
    while ((abs(f.x[2]) > tol) & (n < n.max)) {
        x <- x - f.x[2]/f.x[3]
        f.x <- ftn(x)
        n <- n + 1
        cat("Na iteração", n, "o valor de x foi de", x, "\n")
    }
    if (n == n.max) {
        cat('Método de newton falhou em convergir\n')
    } else {
        cat("O algoritimo convergiu\n")
        return(x)
    }
}
 
###################
#Exemplo 1
#Função gamma(x,2,3)
###################
 
funcao <- function(x) {
    if (x < 0) return(c(0, 0, 0))
    if (x == 0) return(c(0, 0, NaN))
    y <- exp(-2*x)
    return(c(4*x^2*y, 8*x*(1-x)*y, 8*(1-2*x^2)*y))
}
 
entrada<-seq(-5,10,0.01)
saida<-vector()
for(i in 1:length(entrada)) {
    saida[i]<-funcao(entrada[i])[1]
}
 
#Figura 1 e 2
plot(entrada,saida,type="l",frame=F,xlab="x",ylab="f(x)")
abline(v=newton(funcao, 0.25),lty=2,lwd=2,col="red")
abline(v=newton(funcao, 0.5),lty=3,lwd=2,col="blue")
abline(v=newton(funcao, 0.75),lty=3,lwd=2,col="green",h=funcao(newton(funcao, 0.75))[1])
legend("right",lty=c(2,3,3),lwd=2,col=c("red","blue","green"),legend=c(0.25,0.5,0.75),title="Início",bty="n")
 
###################################
#Exemplo 2                        #
#Função de segundo grau invertida #
###################################
 
funcao <- function(x) {
    fx <- -1*(2*x^2+4*x-4)
    dfx <- -1*(4*x+4)
    ddfx<- -1*(4)
    return(c(fx, dfx, ddfx))
}
 
 
entrada<-seq(-6,6,0.01)
saida<-vector()
for(i in 1:length(entrada)) {
    saida[i]<-funcao(entrada[i])[1]
}
 
#Figura 3
plot(entrada,saida,type="l",frame=F,xlab="x",ylab="f(x)")
abline(v=newton(funcao, -2),lty=2,lwd=2,col="red")
abline(h=funcao(newton(funcao, -2))[1],lty=3,lwd=1,col="red")
newton(funcao, 0)
newton(funcao, 1)
 
 
#########################
#Exemplo 3              #
#Função de segundo grau #
#########################
funcao <- function(x) {
    fx <- 2*x^2+4*x-4
    dfx <- 4*x+4
    ddfx<- 4
    return(c(fx, dfx, ddfx))
}
 
 
entrada<-seq(-6,6,0.01)
saida<-vector()
for(i in 1:length(entrada)) {
    saida[i]<-funcao(entrada[i])[1]
}
 
#Figura 4
plot(entrada,saida,type="l",frame=F,xlab="x",ylab="f(x)")
abline(v=newton(funcao, -2),lty=2,lwd=2,col="red")
abline(h=funcao(newton(funcao, -2))[1],lty=3,lwd=1,col="red")
newton(funcao, 0)
newton(funcao, 1)

Root-Finding: O método de Newton-Raphson

Acho que a primeira coisa vem a cabeça, depois de ler o título desse post e lembrando que a gente está num blog de um biólogo é…

“Porque diabos eu quero saber isso?”

É uma boa pergunta, mas veja, quanto mais a gente vai estudando alguma coisa, algum sistema, a gente vai refinando nosso conhecimento, a gente vai fazendo predições melhores, mais específicas, e predições melhores implicam em melhores modelos, e algumas vezes ainda a gente vai cair em modelos não lineares, infelizmente, nem tudo é uma reta. Eu mesmo ja cai aqui nesse caso.

Mas beleza, e qual a relação com esse negócio de achar raiz?

Se você for usar a função nls do R por exemplo, ela tem um argumento chamado algorithm, de um ?nls para ver o help.

Usage nls(formula, data, start, control, algorithm, trace, subset, weights, na.action, model, lower, upper, ...)

Se a gente olhar esse argumento , a gente ve que

algorithm - character string specifying the algorithm to use. The default algorithm is a Gauss-Newton algorithm. Other possible values are "plinear" for the Golub-Pereyra algorithm for partially linear least-squares models and "port" for the ‘nl2sol’ algorithm from the Port library

Temos esse Gauss-Newton como o método default, e se você olhar sobre ele no wikipedia aqui, vai ver a primeira frase.

“It can be seen as a modification of Newton’s method for finding a minimum of a function.”

Na minha tradução nórdica livre, “ele pode ser visto como uma modificação do método de Newton para achar o mínimo de uma função.”

Esse método de achar por sua vez é uma extensão do método de Newton-Raphson.

Newton acho que não precisa de introduções, simplesmente o cara do da teoria da gravidade, e temos esse tal de Raphson, que publicou um ano antes do Newton o método, mas o Newton desenvolveu o método independentemente, e para uso mais genérico.

Então estamos no ponto de partida aqui, entender aqui vai nos permitir entender melhor os métodos de optimização mais complexos, o nosso objetivo, e a hora que a gente pegar o jeito aqui, usar o optmin do R pode ficar mais fácil, e plausível. E usar o optmin, assim como métodos de otimização em geral abrem várias oportunidades de fazer coisas legais. E aqui você pode ler, testar as hipóteses que realmente nos interessam em biologia, e não testar o que da com os dados! Não podemos prender as idéias das pessoas a modelos pre-definidos meramente porque a gente não sabe expandir ou melhorar estes.

Mas vamos começar a embrenhar nesses algorítimos, que no fim das contas são idéias bem legais e úteis.

Para começar, vamos tentar entender o método de Newton-Raphson para achar a raiz de uma função.

Antes de tudo vamos definir melhor nosso objetivo.

Nosso objetivo é achar a raiz de uma função, mas o que é a raiz de uma função?
A raiz de uma determinada função f(x) é o valor de x para o qual f(x)=0, simples assim.

Então existe uma função f qualquer a qual entra um valor real e sai um valor real, f:\mathbb{R} \to \mathbb{R}, e para essa função pode existir um a \in \mathbb{R} tal que f(a)=0, ou sempre existe, não tenho certeza se posso afirmar isso, mas não vamos parar aqui.

Então quando a gente faz aquele velho gráfico, colocando a resposta no eixo y e o preditor no eixo x, a raiz é o lugar do preditor onde a resposta fica 0.

Agora eu concordo que até agora falamos somente blablabla. Porque o interesse nisso?
Suponha que temos uma dívida que tem um valor P, um taxa de juros mensal de r e uma duração de n meses pagar, e um pagamento mensal no valor de A. O valor da dívida após n pagamentos, ou n meses, será P_n.

Então, pensando em tudo isso como uma formula, temos o seguinte:

P_0=P

Primeiro, a gente tem uma dívida inicial P dai a recorrência é a seguinte.

P_{n+1}=P_n \cdot (1+r)-A

Isso é, a cada mês, você paga o juros no referente a quanto você devia mês passado, e então reduz da quantidade da dívida a quantidade A que você sempre paga, essa é uma equação de recorrência de primeira ordem, e tem a seguinte solução:

{P_{n}}={{P \cdot(1+r)^n} -  {{A\cdot((1+r)^n-1 )}\over{r}}}

Colocando P_N=0, ou seja, a divida acabou, temos o seguinte.

{A\over{P}}={r \cdot (1+r)^n \over{(1+r)^n-1}}

Dessa equação, podemos isolar parâmetros, como a taxa de juros (r) por exemplo, suponha que sabemos P, n e A, podemos achar r achando a raiz da seguinte função.

f(x) ={{A\over{P}}-{x \cdot (1+x)^n \over{(1+x)^n-1}}}

Agora acho que ja da para ter uma idéia, que achar raízes tem a ver com otimizar parâmetros certo? Tem a ver não, na verdade, mas pode servir de caminho.
Primeiro a gente vai ver como o método de Newton-Raphson funciona achando raízes, e depois podemos modificar ele um pouco para começar estimar parâmetros em modelos.

Mas vamos la. Agora vamos ver como é que o Newton-Raphson funciona.

Suponha que temos uma função f que é diferenciável com uma derivativa continua f\prime e uma raiz a. Agora se temos um x_0 \in \mathbb{R} como nosso chute inicial de a, a linha entre o ponto (x_0,f(x_0)) com inclinação f\prime(x_0) é a melhor linha aproximação de linha reta da função f(x) para o ponto x_0, esse é o significado da derivativa, a equação para essa linha reta é dada por

{f\prime(x_0)}={{f(x_0)-y}\over{x_0-x}}

Agora, essa linha reta cruza o eixo-x até o ponto x_1, que deve ser uma melhor aproximação que x_0 para a, para achar x_1 nos fazemos

{f\prime(x_0)}={{f(x_0)-0}\over{x_0-x_1}}

logo, a gente tem para x_1 que

 {x_1} = {x_0 - {f(x_0) \over{f\prime(x_0)}}}

De forma mais simples, ja da para sacar o que o algorítimo vai fazer, a gente vai ficar subtraindo  f(x) \over{f\prime(x)} do nosso chute para gerar outro chute, até não aguentar mais.

Ou seja, apos achar x_1, a gente vai atrás do x_2, que vai ser:

 {x_2} = {x_1 - {f(x_1) \over{f\prime(x_1)}}}

Ou de forma geral, podemos escrever esse método da seguinte forma.

 {x_{n+1}} = {x_n - {f(x_n) \over{f\prime(x_n)}}}

Uma relação de recorrência.

Pode ser demonstrado (não que eu consiga) que se f é “bem comportado” em a (de uma olhada aqui para ver o caos, mas se  f\prime(a) = 0 e  f\prime\prime(a) é finito e contínuo em a, vai dar certo) e você começar com x_0 perto o suficiente de a, então x_n vai convergir rápido. Infelizmente, quase sempre a gente não conhece o comportamento de f, a gente pensa nos modelos primeiro, depois vai estudar eles, então a gente não sabe a princípio se esse método vai funcionar ou não.

O método de Newton-Raphson não é garantido de te dar uma resposta, você não tem garantias de uma solução, entretanto, se x_n \to a, então f e f\prime são contínuos e nos temos

{a} = {\lim \limits_{n\to\inf} x_{n+1}} = {\lim \limits_{n\to\inf} ({x_n - {f(x_n) \over{f\prime(x_n)}}})} = {\lim \limits_{n\to\inf} x_n - {f(\lim \limits_{n\to\inf} x_{n}) \over{f\prime(\lim \limits_{n\to\inf} x_{n})}}} = {{a} - {f(a)\over{f\prime(a)}}}

Veja que pode parecer um monte de formula complexa, mas a gente ta substituindo a formula dali de cima dentro dela, dentro dela denovo e denovo, recorrentemente.
Então, desde que f\prime(a) \ne \pm \inf, nos devemos ter f(a) = 0.
Como nos estamos esperando f(x_n) \to 0, um boa condição de parada para o algorítimo é quando \vert f(x_n) \vert \leq Q para uma tolerância Q que nos determinarmos.
Se a sequência \{x_n\}_{n}^{\inf} = {0} está convergindo para a raiz a então para x próximo de a temos que f(x) \approx f(a)\cdot(x-a), então se \vert f(x_n) \vert \leq \epsilon então {\vert (x-a) \vert } \leq {\epsilon\over{f\prime(a)}} aproximadamente.

Agora, depois de todas essa formulas, vamos ver o algorítimo funcionando, e talvez vendo ele funcionando, a gente melhore a nossa crença que as coisas ali em cima funcionam, ou melhor que isso, vendo funcionando a gente melhora nosso entendimento.

newtonraphson <- function(ftn, x0, tol = 1e-9, max.iter = 100) {
    x <- x0
    fx <- ftn(x)
    iter <- 0
 
    while ((abs(fx[1]) > tol) && (iter < max.iter)) {
        x <- x - fx[1]/fx[2]
        fx <- ftn(x)
        iter <- iter + 1
        cat("Na iteração", iter, "o valor de x foi de", x, "\n")
    }
    # A saida depende do que acontece
    if (abs(fx[1]) > tol) {
        cat("O algoritimo falhou em convergir\n")
        return(NULL)
    } else {
        cat("O algoritimo convergiu\n")
        return(x)
    }
}

E é isso ai, pequenininho. O que é nossa entrada? A função, oras bolas, um valor inicial pra começar a recorrência, uma tolerância, pra gente saber onde parar e um número máximo de iterações, isso porque como a gente disse, a gente não tem certeza se toda função tem raiz. Então a gente vai tentar esse tanto de vezes ai, se num deu, a gente para o algorítimo e abandona, vamos pra outro método.

Ai olha a nossa função, a gente começa com

x <- x0
fx <- ftn(x)
iter <- 0

A gente começa o valor de x como x_0
Salvamos em fx, a saída da nossa função, que vai ser uma função que vai calcular o valor de f(x) e o valor de f\prime(x) e devolver o resultado num vetor de dois valores. Lembrando que o R pode facilitar encontrar a derivada com a função deriv, de uma olhada aqui para ver como faz.
E a entrada é a iteração zero.

Dai pra frente é exatamente a mesma coisa.

while ((abs(fx[1]) > tol) && (iter < max.iter)) {
        x <- x - fx[1]/fx[2]
        fx <- ftn(x)
        iter <- iter + 1
        cat("Na iteração", iter, "o valor de x foi de", x, "\n")
    }

Olha so a gente calcula o próximo x com a nossa formulinha de recorrência, dai usando esse valor novo de x e a gente usa de novo nossa função para pegar novos valores de f(x) e f\prime(x) e vamos para uma nova iteração.
A gente vai fazer isso até que ou tenhamos ultrapassado o valor máximo de iterações (para não correr o risco de fazer um loop que não para nunca) ou até não ultrapassarmos nossa tolerância entre as iterações.

Agora vamos fazer um exemplo aqui pra ver se vai dar certo.

Primeiro pegamos uma equação.

f(x)=log(x) - exp(-x)

Bem, depois de uma aulinha de calculo 1 na faculdade ou a matéria de calculo do coursera ou olhar muitos videos do khan academy, a gente sabe que

{f\prime(x)}={{1\over{x}} + {exp(-x)}}

E passando para a linguagem do R a gente faz o seguinte.

funcao <- function(x) {
    fx <- log(x) - exp(-x)
    dfx <- 1/x + exp(-x)
    return(c(fx, dfx))
}

Veja que a entrada é o valor de x, e a saída é um vetor de dois valores, o f(x) e f\prime(x) respectivamente.

Usando nossa o método do Newton-Raphson temos o seguinte.

> newtonraphson(funcao, 2, 1e-06) Na iteração 1 o valor de x foi de 1.12202 Na iteração 2 o valor de x foi de 1.294997 Na iteração 3 o valor de x foi de 1.309709 Na iteração 4 o valor de x foi de 1.3098 O algoritimo convergiu [1] 1.3098

Bang, f(1.3098) \approx 0

A gente pode olhar num gráfico isso pra ficar mais fácil.

01

Olha ai, quando o x é 1.3098, no y a gente ta no zero.

A gente pode colocar esse número na função que fizemos também.

> funcao(1.3098) [1] 4.280091e-07 1.033349e+00

Lembrando que o f(x) é o primeiro número e esta em notação científica, ou seja

> format(4.280091e-07,scientific=F) [1] "0.0000004280091"

E isso é bem próximo de zero, se quiser mais precisão, é só diminuir a tolerância, mais ai precisara de mais iterações, para mais precisão, e isso depende na verdade da aplicação que a gente vai dar para esse número. E em ultima instancia, estamos limitados ao menor número que o computador pode representar, e isso não é infinito.

Agora veja como é a derivativa da função

02

Veja que a derivativa nos da a direção na qual mudar o valor de x, a gente sabe se o valor de x ta mudando, e como está mudando.

Veja que imprimimos os valores que fomos usando, então podemos olhar como foram as iterações.
Ou é so modificar a função e guardar em algum lugar essa informação, transformando a saída numa lista por exemplo.

03

Veja que começamos no 2, ai diminuímos bastante, e fomos voltando, cada vez mais devagarinho, veja que as setas vão diminuindo, até chegar num momento que não da mais nem para desenhar setinhas.

Certo, mas vai que foi sorte, vamos tentar para outra função, a boa e velha x^2, pra ela temos o seguinte

{f(x)}={x^2}

e a derivada é fácil, é

{f\prime(x)}={2 \cdot x}

Passamos essa função pro R, pra poder usar o método de Newton-Raphson

funcao <- function(x) {
    fx <- x^2
    dfx <- 2*x
    return(c(fx, dfx))
}

e vamos usar ela

> newtonraphson(funcao, 10, 1e-06) Na iteração 1 o valor de x foi de 5 Na iteração 2 o valor de x foi de 2.5 Na iteração 3 o valor de x foi de 1.25 . . . Na iteração 12 o valor de x foi de 0.002441406 Na iteração 13 o valor de x foi de 0.001220703 Na iteração 14 o valor de x foi de 0.0006103516 O algoritimo convergiu [1] 0.0006103516

Veja que a diferença da iteração 13 para 14 é menor que a nossa tolerância

> format(1e-06,scientific=F) [1] "0.000001"

Veja o último número que muda durante as iterações.
Agora o lado bom dessa função é que a gente tem outro método pra achar raiz, para testa o que estamos usando, podemos usar a boa e velha bhaskara.

Se a gente usa bhaskara, a gente acha as duas raízes, e ambas são zero, veja que a entrada da função são os coeficientes, a, b e c.

> bhaskara(1,0,0) [1] 0 0

A gente pode olhar a solução também num gráfico como no exemplo anterior.

04

Mas e se tivesse duas raízes? Vamos pegar uma segunda função de segundo grau aqui, que dai a gente pode conferir com bhaskara.

Vamos usar essa função

f(x)=2x^2+4x-4

a derivada dela é bem fácil também

f\prime(x)=4x+4

E no R fica assim

funcao <- function(x) {
    fx <- 2*x^2+4*x-4
    dfx <- 4*x+4
    return(c(fx, dfx))
}

Primeiro vamos ver onde queremos chegar usando bhaskara

> bhaskara(2,4,-4) [1] 0.7320508 -2.7320508

Dai vamos usar o método de Newton-Raphson.

> newtonraphson(funcao, 10, 1e-06) Na iteração 1 o valor de x foi de 4.636364 Na iteração 2 o valor de x foi de 2.084311 Na iteração 3 o valor de x foi de 1.028488 Na iteração 4 o valor de x foi de 0.753711 Na iteração 5 o valor de x foi de 0.7321846 Na iteração 6 o valor de x foi de 0.7320508 O algoritimo convergiu [1] 0.7320508

E veja que conseguimos a primeira raiz, mas veja que começamos chutando o valor de 10, dai encontramos a raiz mais próxima desse chute inicial, isso quer dizer que se chutarmos um número mais próximo a outra raiz, acharemos ela

> newtonraphson(funcao,-10, 1e-06) Na iteração 1 o valor de x foi de -5.666667 Na iteração 2 o valor de x foi de -3.654762 Na iteração 3 o valor de x foi de -2.892403 Na iteração 4 o valor de x foi de -2.738845 Na iteração 5 o valor de x foi de -2.732064 Na iteração 6 o valor de x foi de -2.732051 O algoritimo convergiu [1] -2.732051

E pimba, temos as duas raízes, e com a certeza que o algorítimo funciona beleza, já que já sabíamos o que esperar.
Podemos também ver o método no gráfico.

05

Mas então porque não usar sempre bhaskara? Ai a gente cai na mesma pergunta, de porque não usar sempre funções lineares nos modelos, porque as vezes não da.
Vamos a um último exemplo, para as funções esquisitas.

Pense na função

f(x)=exp(x)

espantosamente, a derivada dela é

f\prime(x)=exp(x)

Não vamos debater isso, mas é legal né, a derivada é a própria função, aqui bhaskara não funciona, mas Newton-Raphson pode dar certo.

Estabelecemos nossa função no R

funcao <- function(x) {
    fx <- exp(x)
    dfx <- exp(x)
    return(c(fx, dfx))
}
> newtonraphson(funcao, 10, 1e-06,100) Na iteração 1 o valor de x foi de 9 Na iteração 2 o valor de x foi de 8 Na iteração 3 o valor de x foi de 7 . . . Na iteração 23 o valor de x foi de -13 Na iteração 24 o valor de x foi de -14 O algoritimo convergiu [1] -14

E, pimba, deu certo também.

Podemos olhar no gráfico para ver se esse número faz sentido

06

E parece que está tudo certinho.

E esse é o método Newton-Raphson. Como a gente já viu no algorítimo do gradient descent, daqui, com algumas modificações, a gente vai parametrizar modelos também.

Agora imagine que o Newton, um cara que olhava planeta, tinha dados de posição de planeta e parametrizava as coisas fazendo contas na mão, ainda bem que esse método não precisa de sei la, 10 mil iterações, senão ele tava ferrado hehehe.

Mas é isso ai, o primeiro algorítimo de encontrar raízes que vemos, num próximo post nos progredimos ele um pouco para determinar parâmetros.

Bem, eu vou ficar por aqui, que já faz muito tempo que estou escrevendo aqui. Scripts no repositório recologia e aqui em baixo, e até o próximo post, este aqui espero dar sequência e usar o método de Newton para achar o mínimo de funções, ou seja, estimar parâmetros pela máxima verossimilhança.

Referencia:
Owen Jones, Robert Maillardet & Andrew Robinson 2009 – Introduction to Scientific Programming and Simulation Using R. CRC Press 499pp

bhaskara<-function(a,b,c) {
    delta = b*b - 4*a*c
    if(delta < 0) {
        print("A equacão não possui raizes reais.\n");
        saida<-c(NA,NA)
        return(saida)
    }else{
        x1 = (-b + sqrt(delta)) / (2*a)
        x2 = (-b - sqrt(delta)) / (2*a)
        saida<-c(x1,x2)
        return(saida)
    }
}
 
newtonraphson <- function(ftn, x0, tol = 1e-9, max.iter = 100) {
    x <- x0
    fx <- ftn(x)
    iter <- 0
 
    while ((abs(fx[1]) > tol) && (iter < max.iter)) {
        x <- x - fx[1]/fx[2]
        fx <- ftn(x)
        iter <- iter + 1
        cat("Na iteração", iter, "o valor de x foi de", x, "\n")
    }
    # A saida depende do que acontece
    if (abs(fx[1]) > tol) {
        cat("O algoritimo falhou em convergir\n")
        return(NULL)
    } else {
        cat("O algoritimo convergiu\n")
        return(x)
    }
}
 
##########################
#Exemplo 01
##########################
funcao <- function(x) {
    fx <- log(x) - exp(-x)
    dfx <- 1/x + exp(-x)
    return(c(fx, dfx))
}
 
newtonraphson(funcao, 2, 1e-06)
 
funcao(1.3098)
 
format(4.280091e-07,scientific=F)
 
#figura 1
curve(log(x) - exp(-x),0,3,frame=F)
points(1.3098,funcao(1.3098)[1],cex=2)
abline(v=1.3098,h=funcao(1.3098)[1],lty=3,col=2)
legend("topleft",legend="Raiz",pch=1,bty="n")
 
#figura 2
jpeg("02.jpg")
par(mfrow=c(2,1))
curve(log(x) - exp(-x),0,3,frame=F,main="Função")
abline(v=1.3098,h=funcao(1.3098)[1],lty=2)
curve(1/x + exp(-x),0,3,frame=F,main="Primeira derivada")
abline(v=1.3098,lty=2)
 
#figura 3
pontos<-c(2,1.12202,1.294997,1.309709,1.3098)
plot(pontos,log(pontos)-exp(-pontos),type="p",pch=19,cex=0.5,frame=F,xlab="x",ylab="f(x)")
text(pontos[1]-0.1,(log(pontos)-exp(-pontos))[1],"inicio")
#ps, a ultima seta não tem tamanho suficiente para o gráfico
for(i in 1:5) {
    arrows(pontos[i],(log(pontos)-exp(-pontos))[i],pontos[i+1],(log(pontos)-exp(-pontos))[i+1])
}
 
###########################3
#Exemplo 2
############################
funcao <- function(x) {
    fx <- x^2
    dfx <- 2*x
    return(c(fx, dfx))
}
 
newtonraphson(funcao, 10, 1e-06)
format(1e-06,scientific=F)
 
bhaskara(1,0,0)
 
#figura 4
par(mfrow=c(2,1))
curve(x^2,-5,5,frame=F,main="Função")
abline(v=0,h=0.0006103516,lty=3,col=2)
curve(2*x,-5,5,frame=F,main="Primeira derivada")
abline(v=0.0006103516,lty=3,col=2)
 
 
###################
#Exemplo 3
###################
funcao <- function(x) {
    fx <- 2*x^2+4*x-4
    dfx <- 4*x+4
    return(c(fx, dfx))
}
 
bhaskara(2,4,-4)
 
newtonraphson(funcao, 10, 1e-06)
newtonraphson(funcao,-10, 1e-06)
 
#figura 5
par(mfrow=c(2,1))
curve(2*x^2+4*x-4,-10,10,frame=F,main="Função")
abline(h=0,lty=3,col=2)
abline(v=0.7320508,lty=3,col=2)
abline(v=-2.732051,lty=3,col=2)
curve(4*x+4,-10,10,frame=F,main="Primeira derivada")
abline(v=0.7320508,lty=3,col=2)
abline(v=-2.732051,lty=3,col=2)
 
###############
#exemplo 4
###############
 
funcao <- function(x) {
    fx <- exp(x)
    dfx <- exp(x)
    return(c(fx, dfx))
}
 
newtonraphson(funcao, 10, 1e-06,100)
 
#figura 6
par(mfrow=c(2,1))
curve(exp(x),-15,-5,main="Função exp(x)",frame=F)
abline(h=0,v=-14,lty=3,col=2)
curve(exp(x),-15,-5,main="Primeira derivada de exp(x)",frame=F)
abline(v=-14,lty=3,col=2)

Insertionsort em R e C++ usando o pacote rcpp

Bem, depois de vermos o algorítimo mais simples que tem para ordenar um vetor, o bubble sort. Agora a gente vai olhar outro um pouco mais legal, esse se chama insertion sort.

A ideia é a seguinte:

Então a gente tem um vetor qualquer, com n elementos.
 C = \{X_1 , X_2 , X_3 , \dots , X_n  \}

E a gente quer que esses elementos fiquem assim:

 Regra = \{x_1 \leq x_2 \leq x_3 \leq \dots \leq x_n  \}

Como a gente faz?
Vamos supor um conjunto aqui.

 C = \{6, 5, 3, 1, 8, 7, 2, 4\}

Ele esta desorganizado porque ele não respeita nossa Regra_1. So olhar o primeiro elemento que é 6.

So que a estratégia aqui muda. O que a gente faz?
Primeiro pega um sub-vetor, um pedaço desse la em cima.

 C = \{6\}

Esse vetor está organizado? Sim, não temos ninguém para comparar. Um vetor de 1 elemento sempre está organizado.

Agora pegamos o segundo elemento, o 5.
Ai fazemos a pergunta onde ele deveria está?
Antes ou depois do 6?
Bem ele é menor que o 6 então tem que estar depois, então movemos o 6 um índice pra frente e dai colocamos o 5.

 C = \{5 ,6\}

Agora pegamos outro elemento, o 3, e repetimos o os passos anteriores.

Perguntamos, onde é a posição do 3?
Ele é menor que o 6? Sim então movemos o 6 um índice para a frente.
Ele é menor que o 5? Sim, então movemos o 5 um índice para a frente, como estamos na primeira casa do vetor, paramos. E ficamos com o seguinte.

 C = \{3 ,5 ,6\}

Notem que o que temos é que o inicio do vetor sempre esta organizado dessa forma, e sempre adicionamos um elemento a um vetor organizado. Mas se o elemento que formos adicionar deve ficar na ultima posição, a gente não vai mover nada de lugar, ou seja, vai estar tudo certo e teremos organizado com menos trocas de posição de elementos. Diferente do bubblesort, que tínhamos que comparar o vetor inteiro para garantir que o ultimo elemento é o maior, mas depois temos que refazer varias comparações para garantir que o segundo maior elemento é o penúltimo, ou seja nunca tinha como escapar de comparar tudo com tudo, aqui tem.

Mas talvez vendo uma animaçãozinha fica mais fácil de entender.

Insertion-sort-example-300px

Apesar de parecer complicado, esse é o jeito que a maioria das pessoas organiza cartas na mão, quando a gente quer deixar em alguma ordem. Veja so.

Ok, mas então como fica esse processo no R?

insertionsort<-function(vetor){
    n<-length(vetor)
 
    for(i in 2:n) {
        aux=vetor[i]
        j=i-1
        while(j>=1 && vetor[j]>aux) {
            vetor[j+1]<-vetor[j]
            j=j-1
            }
        vetor[j+1]=aux
        }
    return(vetor)
    }

Bem simples, um loop usando i diz qual o tamanho do vetor que já está organizado, veja que a gente começa do 2, ja que o índice 1 já está organizado, e ai é só salvar a ultima posição, usamos um while para quanto o auxiliar for menor que o valor no vetor ir empurrando ele para a frente e quando isso não é verdade, colocamos o valor que salvamos no auxiliar na posição desejada.

E esse processo funciona que uma beleza.

> vetor<-sample(100) > vetor [1] 78 2 25 43 10 56 8 64 49 41 32 60 59 5 3 38 69 46 13 77 70 57 4 [24] 26 51 45 53 73 85 55 36 66 83 80 61 79 48 50 62 42 35 23 97 19 27 20 [47] 65 94 54 81 1 82 75 93 86 24 71 63 100 34 39 28 6 7 30 33 31 11 72 [70] 18 99 16 17 96 95 52 76 14 89 22 84 98 91 67 37 87 40 88 21 58 92 29 [93] 12 68 15 90 74 47 9 44 > insertionsort(vetor) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 [47] 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 [70] 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 [93] 93 94 95 96 97 98 99 100

Podemos fazer isso também usando o Rcpp, mas ai é preciso lembrar, que na linguagem C e C++, diferente do R, os índices começam no 0 e não em 1, então o primeiro elemento de um vetor está na posição vetor[1] no R, mas esta na posição vetor[0] no C++.
Com isso em mente, o código é praticamente o mesmo.

library(Rcpp)
cppFunction("
    NumericVector insertionsortC(NumericVector vetor) {
        int n = vetor.size();
 
        double aux;
        int i , j;
 
        for(i=1;i<n;i++) {
            aux=vetor[i];
            j=i-1;
            while(j>=0 && vetor[j]>aux) {
                vetor[j+1]=vetor[j];
                j=j-1;
                }
            vetor[j+1]=aux;
            }
        return vetor;
        }
")

E funciona legal também.

> insertionsortC(vetor) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 [47] 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 [70] 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 [93] 93 94 95 96 97 98 99 100

Bem a gente vai ver daqui dois algorítimos de ordenação, se eu continuar fazendo isso, que o que a gente manda para o Rcpp, o código que está entre aspas, pode ser um arquivo texto com o código. E poderíamos mandar mais funções juntas, mas no quicksort a gente vai ver como organizar isso, por enquanto deixa desse jeito que é mais simples.

Mas outra forma ainda de implementar isso é com recursão. A gente ja usou recursão para resolver problemas no Rosalind. Mas basicamente é so chamar a mesma função dentro dela mesma. Um exemplo bem simples é o fatorial.

Fatorial é aquele número que multiplicamos por ele menos 1 até chegar a 1.
Assim,  3!=3 \cdot 2 \cdot 1 e 2 fatorial seria  2!=2 \cdot 1

Então outra forma de escrever isso é  3!=3 \cdot 2!. Fazendo isso até chegar no 1 fatorial, logo no R isso fica assim.

fatorial<-function(n) {
    if(n==1) {
        return(1)
        } else {
            return(n*fatorial(n-1))
            }
    }

Que da.

> fatorial(3) [1] 6

Ou seja, temos um caso base, que é o 1 e vamos fazendo alguma coisa até chegar nessa base.
O insertionsort pode seguir esse raciocínio também da seguinte forma:

cppFunction("
    NumericVector insertionsortRC(NumericVector vetor, int n) {
 
        double aux;
        int i;
 
        if(n>1) {
            insertionsortRC(vetor,n-1);
            aux=vetor[n-1];
            i=n-1;
            while(vetor[i-1]>aux && i>=0 ) {
                vetor[i]=vetor[i-1];
                i--;
                }
            vetor[i]=aux;
            }
 
        return vetor;
        }
    ")

Que funciona também

> insertionsortRC(vetor,length(vetor)) [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 [24] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 [47] 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 [70] 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 [93] 93 94 95 96 97 98 99 100

Agora eu so conseguir fazer o insertionsort recursivo fazendo uma função que recebe dois argumentos. Eu até perguntei no stack-overflow aqui se tinha como fazer de outra forma, mas não sei se é possivel, mas tudo funciona de qualquer forma.
Que descobri que não sabia fazer foi fazer o insertion sort de forma recursiva no R. Tem uma pergunta aberta aqui na lista brasilera de r r-br. Se alguém puder me explicar como faz, eu fico grato.

E bem eu acho que é isso, agora temos dois algorítimos de ordenação. Vamos adicionar mais alguns a essa lista e depois será divertido tentar analisar o quão bom cada um é, e quanto…

Bem o script está aqui embaixo além do repositório recologia no github. Até a próxima.

vetor<-sample(100)
vetor
 
insertionsort<-function(vetor){
    n<-length(vetor)
 
    for(i in 2:n) {
        aux=vetor[i]
        j=i-1
        while(j>=1 && vetor[j]>aux) {
            vetor[j+1]<-vetor[j]
            j=j-1
            }
        vetor[j+1]=aux
        }
    return(vetor)
    }
 
insertionsort(vetor)
 
library(Rcpp)
 
cppFunction("
    NumericVector insertionsortC(NumericVector vetor) {
        int n = vetor.size();
 
        double aux;
        int i , j;
 
        for(i=1;i<n;i++) {
            aux=vetor[i];
            j=i-1;
            while(j>=0 && vetor[j]>aux) {
                vetor[j+1]=vetor[j];
                j=j-1;
                }
            vetor[j+1]=aux;
            }
        return vetor;
        }
")
 
 
 
insertionsortC(vetor)
 
cppFunction("
    NumericVector insertionsortRC(NumericVector vetor, int n) {
 
        double aux;
        int i;
 
        if(n>1) {
            insertionsortRC(vetor,n-1);
            aux=vetor[n-1];
            i=n-1;
            while(vetor[i-1]>aux && i>=0 ) {
                vetor[i]=vetor[i-1];
                i--;
                }
            vetor[i]=aux;
            }
 
        return vetor;
        }
    ")
 
vetor
insertionsortRC(vetor,length(vetor))
 
fatorial<-function(n) {
    if(n==1) {
        return(1)
        } else {
            return(n*fatorial(n-1))
            }
    }
 
fatorial(3)