As vezes, menos é melhor.

A gente ja viu vários posts com vários tipos de modelos, de poisson, binomial, mixtos e até expansões como overdispersion e zeros inflados.
Mas antes de mais nada, antes de sair tentando tudo e ver o que da mais certo. A gente tem que pensar, será que os parâmetros que estamos estimando fazem sentido. A gente consegue discutir aquele monte de número? Eles tem algum sentido biológico (para quem é da biologia)?

Vamos olhar um conjunto de dados, contagens de uma população ao longo do tempo, pensa que vemos o seguinte.

01

Não precisamos mentir. A maioria das pessoas (inclua eu nessa) pensaria em tacar uma regressão linear nesse conjunto de dados, um pensamento instantâneo, como se estivéssemos condicionados a fazer isso. Não que isso seja totalmente errado. Vamos ver como esse modelo ficaria.

02

Olha, não ficaria tão ruim, até que o modelo ajustaria. Mas e ai, o que estaríamos fazendo? Vendo como a população dessa espécie muda ao longo do tempo? E a inclinação é a intensidade como a população está mudando e o intercepto é como ela começou?

Eu acho que não é um interpretação tão ruim assim. Mas podemos melhorar mais e mais esse modelo, usando polinômios por exemplo.
Polinômios de segundo grau são uma parábola, de terceiro grau fica com uma curvinha, e cada vez o modelo vai ficando melhor, até que com 6 graus temos o seguinte.

03

Agora ta perfeito. Mas infelizmente, perfeito somente para esse grupo de dados, conforme a gente aumenta o número de parâmetros, a gente sempre melhora o ajuste, diminui o espaço para erro de medida, o problema é que o que os parâmetros fazem? O que significam? Nada mais, a gente so tem um ajuste bom para um modelo ininterpretável.
Olhe o que criamos, ele pode parecer bonito no grafico, mas não significa muito. E pense que fazemos modelos para pensar nos dados que coletamos, mas nos dados também que não coletamos, que outras pessoas vão coletar, e mais que isso, para entender como o mundo funciona, e será que temos muita capacidade de predição com esse modelo?
Se parar para pensar por um momento, a gente lembra que existem modelos de crescimento populacional que vão além de retas e linhas. Que cada letrinha, cada parametro tem um sentido especial que discutimos um bocado por aqui. Então se a gente passa tanto tempo lendo sobre eles, deveríamos gastar eles. até porque dai fica fácil discutir os dados.

04

Olha ae, talvez ele não seja tão bonito quanto ao ajuste, como o modelo polinimial, talvez até pior que o ajuste da regressão linear, mas é muito mais o que pensamos, e temos que pensar mais para frente, no tempo a frente.

05

Veja o que preveríamos para o futuro de acordo com cada modelo. O modelo de regressão linear preveria um crescimento constante, e o pior de tudo, o modelo polinomial preveria uma queda brusca na população. Agora você acha que tem algum evento obvio que faria a população diminuir. No fim das contas o modelo de crescimento populacional é o que faz a previsão mais razoável, que a população vai continuar crescendo. E é o que a gente acha. Os modelos de regressão são muito bons localmente, ali aonde temos dados, mas eles são péssimos fora da região onde estão os dados. Isso é algo que é bem comum, porque a regressão não assume muita coisa sobre o sistema biológico em si, o preço que se paga é isso.
Claro que se a gente não sabe nada, seria uma informação importante saber o que está acontecendo nesse período de tempo, e como está acontecendo, mas enquanto estamos construindo um conhecimento “bottom-up”, não sabemos nada, e estamos começando. Agora se a gente já tem modelos, a gente pode por eles a prova, algo mais “top-down”. A gente parte do que já sabe e vê se funciona na prática.

Agora quando falamos de chi-quadrado aqui, falamos também de como usar ele para seleção de modelos, mas aqui ferrou, porque para a seleção de modelos, precisamos de modelos aninhados. Por exemplo, modelos de Evolução de DNA temos vários modelos, mas um é uma versão mais simplificada de outro, um caso especial. A mesma coisa temos com os modelos de regressão linear, polinomial ou não. Agora o modelo de crescimento polinomial não aninha com os modelos de regressão, então acho que temos que escolher na argumentação, e não vai ter um valor de p dando sopa para justificar sem discussão essa escolha. Mas é a hora que pensar, apesar de no começo da curva de crescimento populacional (a vermelha) ser ruim no começo, lembramos que no começo a população é pequena, e coisas ao acaso acontecem mais comumente, so lembrar de genetic drift. Depois a população deve estabilizar mais e a curva fica melhor. Muito melhor porque claramente a população ta começando a crescer, e vai crescer bastante, e o modelo de crescimento populacional é o único que está dizendo isso.

Bem é isso ai, uma semana sem posts por conta de provas, mas agora voltando a atividades :), os scripts sempre no repositório recologia e é isso, até o próximo post.

tempo<-seq(0,12,by=2)
N<-c(230,1205,1150,2091,3010,6106,7323)
 
#figura 1
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,pch=19)
 
#modelo de regressão linear
modelo1<-lm(N~tempo)
modelo1
 
#figura 2
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,8000),pch=19)
curve(coef(modelo1)[1]+coef(modelo1)[2]*x,0,12,add=T,lwd=2,lty=1,col="green")
segments(x0=tempo, y0=N, x1 =tempo , y1 = predict(modelo1))
 
#modelo de regressão polinomial
modelo2<-lm(N~poly(tempo, degree=6, raw=TRUE))
modelo2
 
#figura 3
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,8000),pch=19)
curve(coef(modelo2)[1]+coef(modelo2)[2]*x+coef(modelo2)[3]*x^2+coef(modelo2)[4]*x^3+coef(modelo2)[5]*x^4+
      coef(modelo2)[6]*x^5+coef(modelo2)[7]*x^6,0,12,add=T,lwd=2,lty=1,col="blue")
segments(x0=tempo, y0=N, x1 =tempo , y1 = predict(modelo2))
 
#modelo não linear
modelo3<-nls(N~N0+exp(r*tempo),start=list(N0=1,r=1))
modelo3
 
#figura 4
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,8000),pch=19)
curve(coef(modelo3)[1]+exp(coef(modelo3)[2]*x),0,12,add=T,lwd=2,col="red")
segments(x0=tempo, y0=N, x1 =tempo , y1 = predict(modelo3))
 
#figura 5
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,10000),xlim=c(0,20))
curve(coef(modelo1)[1]+coef(modelo1)[2]*x,0,20,add=T,lwd=2,lty=1,col="green")
curve(coef(modelo2)[1]+coef(modelo2)[2]*x+coef(modelo2)[3]*x^2+coef(modelo2)[4]*x^3+coef(modelo2)[5]*x^4+
      coef(modelo2)[6]*x^5+coef(modelo2)[7]*x^6,0,20,add=T,lwd=2,lty=1,col="blue")
curve(coef(modelo3)[1]+exp(coef(modelo3)[2]*x),0,20,add=T,lwd=2,col="red")
points(N~tempo,pch=19,cex=1.2)
legend("topleft",col=c("green","blue","red"),lwd=2,bty="n",title ="Modelo",
       legend=c("Regressão Linear","Polinomial de grau 6","Crescimento populacional"))
 
logLik(modelo1)
logLik(modelo2)
logLik(modelo3)

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)

Snake and Ladders, usando Markov Chain para entender um jogo.

Snake and Ladders, ou Cobras e Escadas em português, é um jogo de tabuleiro originário da índia. Acho que é parecido com a maioria dos jogos de tabuleiros de hoje em dia que so se anda para frente, só mudando a temática, enquanto no Snake and Ladders, você usa uma escada para avançar casas e uma cobra te engole e, bem, você sai pela outra extremidade casa atras, no jogo da vida você pega carona para frente, ou volta no tabuleiro como punição. No wikipedia tem um tabuleiro dele nas antigas.

02

Mas qual o interesse nisso? Olha como o mundo científico é um lugar pequeno. Lembra do senhor John Conway, do jogo da vida que discutimos aqui? Esse mesmo cara escreveu um livro chamado Winning Ways for your Mathematical Plays, e nesse livro ele usa um tipo de cadeia de Markov chamado Absorbing Markov chain para modelar esse joguinho, e aqui vamos olhar esse exemplo.

Bem o jogo é simples, temos esse tabuleiro, mas na nossa simulação vamos usar um tabuleiro um pouco modificado, um com 100 casinhas. Bem para facilitar uma representação no R, vamos usar setinhas ao invés de Escadas e cobras, ou seja se você cair uma no início da seta(lado contrário a ponta) da seta, você vai para onde ela ta apontando, se você cair na ponta da seta, continue sua vida.

Então aqui está nossa representação.

01

Então todo mundo começa com o peão fora do tabuleiro, então joga um dado de 6 faces e vê onde vai parar. Acontece que podemos ver isso como um processo markoviano, em que cada posição no tabuleiro é um estado do sistema, e temos, a cada jogada, as respectivas chances de mudar de estado. So que a gente não vai para qualquer estado, a gente tem uma chance em seis de mover de estado do quadradinho atual para o próximo, uma chance em seis de andar 2 casas, 1 chance em seis de andar 3 casas, até seis.

Isso pode ser representado por uma matriz de estados, com 100 linhas e 100 colunas, uma para cada posição do jogo. 101 na verdade, poque vamos supor que o estado 0 é quando o jogo não começou, o peão não esta ainda no tabuleiro, mas ele tem que começar de algum lugar, que é a posição zero. Assim, a primeira linha vai ser a seguinte.

> M[1,] 0 1 2 3 4 5 6 7 8 9 0.0000000 0.1666667 0.1666667 0.1666667 0.0000000 0.1666667 0.1666667 0.0000000 0.0000000 0.0000000 . . . 90 91 92 93 94 95 96 97 98 99 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 100 0.0000000

Vou suprimir um pouco dos resultado para ficar mais fácil de entender.
Mas veja que essa é a primeira linha da matriz, e ela ta dizendo que uma vez que estamos na posição zero, temos uma chance em seis, 1 divido por 6 que da 0.1666667 arredondando, podemos ir para na casa 1, 2 ,3, 5 ,6. Opa, porque o 4 está com zero de chance? Porque, se você olhar no tabuleiro ali em cima, temos que quando a gente cai no 4, a gente pega a escadinha e vai parar na casa 14.

Podemos colocar esse números na nossa representação gráfica do jogo.

02

Certo, mas a primeira jogada é fácil, mas daqui pra frente que vem a parte bolada.
Se você está na casa 3, e jogar o dado de 6 faces denovo, existem 6 possibilidades ainda.

02

Mas na verdade se a gente jogar o dado duas vezes, podemos saber qual a chance de parar em qualquer lugar ainda, simplesmente multiplicando as chances de cada valor de cada um dos dois dados jogados.

01

Ja vimos isso aqui, em probabilidade condicional. So que a matriz de estados aqui vai respeitar essa propriedade também.

Então podemos fazer uma função para multiplicar linhas dessa matriz.

powermat=function(P,h) {
    Ph=P
    if(h>1) {
        for(k in 2:h) {
            Ph=Ph%*%P
        }
    }
    return(Ph)
}

Essa função multiplica linhas da nossa matriz de transição, e multiplicando as duas primeiras linhas vemos o seguinte.

03

Veja que estamos representando a intensidade das probabilidades com cores, quanto mais vermelho, maior a chance, quanto mais clarinho, menor a chance. Aqui a chance maior é de estarmos na posição sete, já que podemos chegar no sete tirando 1 e 6, 2 e 5, 3 e 4, 4 e 3, opaaaaaa, lembre-se que quando a gente cai no 4, a gente pega a escada para o 14, então com 4 e 3, vamos para o catorze, e depois andamos 3, só que no 17 a gente é engolido pela cobra e volta pro 7, eu diria que isso é uma ironia, já que era para estarmos no 7 mesmo. Mas veja que como temos a matriz, estamos levando em conta as escadas e cobras, o que ignoraríamos, ou teríamos bem mais dificuldade se formos no braços, simulando cada acontecido.

É fácil se perder, se a gente andar 10 vezes, podemos andar de 6 em 6 casas, com muita sorte, você esperaria que ninguém tenha ganho o jogo não? Afinal andando direto, isso da 60 casas, se você começa no zero, so tirando 6 estaria na casa 60, mas lembre-se que quem tem sorte, não vai tirar só 6, vai é pegar as escadas certas, e andar muito mais rápido que isso, veja a distribuição depois de 10 jogadas como fica.

04

As vezes da a impressão que tem algo de especial nessa matriz, mas confira, a gente tem que cada casa só tem chance de ir para 6 outras casas, só voltando em casos de quando encontramos a boca de uma cobra.

> sum(M[2,]>0) [1] 6 > sum(M[10,]>0) [1] 6 > sum(M[50,]>0) [1] 6

Agora vamos supor que depois dessas 10 jogadas, um cara que não sabe perder da um chute na mesa e derruba todos os peões, agora a gente não lembra onde estava, mas podemos pegar a figura acima e podemos então tentar olhar onde você estava, primeiro que sabemos exatamente onde você não pode estar, só ver os lugares com probabilidade zero.

Mais do que isso, se você afirma que estava em algum lugar entre 59, 60 ou 61, podemos ver a probabilidade de estar em cada um desses lugares.

> h=10 > posição<-(initial%*%powermat(M,h))[59:61]/sum((initial%*%powermat(M,h))[59:61]) > posição [1] 0.1597003 0.5168209 0.3234788 > h=10 > posição<-(initial%*%powermat(M,h))[59:61]/sum((initial%*%powermat(M,h))[59:61]) > sum(posição) [1] 1 > posição [1] 0.1597003 0.5168209 0.3234788

Veja que a gente multiplica onde podemos estar depois de 10 jogadas, olhas a probabilidade de cada uma dessas 3 posições, divido pela soma dessas 3 posições, e temos nossa distribuição para a sua afirmação de estar nas casas 59, 60 ou 61.

Veja que podemos seguir a mesma linha, e perguntar quando ainda estaremos jogando, se a gente olhar a distribuição para cada jogada, é só ver que a chance da gente estar jogando ainda, é a chance de não estarmos na posição 100. Veja que uma vez que caímos no estado “posição 100”, ele é estacionário, e todo o jogo tende pra ele, se a gente chegar la, a gente não sai mais, isso é o jogo acabou, então podemos olhar a chance de não estarmos no 100 jogada por jogada e construir uma curva de probabilidade para isso.

05

Após 50 jogadas, é praticamente certeza que o jogo acabou. Veja a distribuição de chances depois de 50 jogadas.

06

Tem que ser muito pé frio para não ter terminado o jogo depois de 50 jogadas, mas veja que isso é uma situação pouco provável, mas não impossível.

Mas em média, vão levar 32 jogadas para finalizar o jogo.

> sum(1-game) [1] 32.16499

Veja que 1 é 100% de chance do jogo ter terminado, quando estamos diminuindo 100% de chance, o que está sobrando é a chance de ainda estarmos jogando. Quando o jogo acabou com certeza, temos 1, ou 100% de chance do jogo ter acabado, então 100% – 100% da zero, a chance da zerada, então estamos somando ai em cima todas as chances de estarmos jogando, para a distribuição ali em cima. Assim talvez seja mais difícil pensar.

Para ficar mais simples, podemos perguntar, quando na distribuição, temos menos de 50% de chance de estar jogando, 50% de chance do jogo ter acabado.

> max(which(1-game>.5)) [1] 29

São vinte e nove jogadas, e já temos 50% de chance do jogo ter acabado.

07

Agora tudo isso que estamos falando, é tendo em mente apenas um jogador. E com dois jogadores, mudaria muito?

Bem, é intuitivo pensar que, quanto mais gente jogando, mais rápido o jogo deve acabar. Se temos dois jogadores, é mais difícil exatamente termos dois jogadores azarados.

08

É só fazer a potência da distribuição anterior, elevar ela pelo número de jogadores, no caso dois, que temos a nova distribuição de quando o jogo deve acabar, podemos calcular também qualquer uma das estatísticas anteriores.

Mais do que isso, podemos usar isso para qualquer quantidade de jogadores.

09

Bem, para quem não pegou ainda a ideia das cadeias de markov, esse foi mais um exemplinho bem legal.
Esse não é o primeiro exemplo de Markov Chain que a gente olha, ja vimos aqui como prever a chuva usando Markov Chain, ou ainda como podemos prever a chuva, sem ver a chuva, apenas pelo comportamento de alguém, nesse exemplo está aqui. Isso pode parecer todo, mas o crescimento populacional pode ser pensado dessa forma, o que a gente ve pode ser apenas uma parte do processo que realmente esta acontecendo, como vimos no crescimento de passarinhos aqui.

Outro processo que temos muito interesse e que é assim é a evolução, a gente viu como simular cadeias de DNA com cadeias de markok aqui. E a evolução é mais ou menos isso, é uma passagem quase sempre, so de ida, não tem volta, a gente dificilmente volta para um estado anterior.

Bem, acho que é isso ai por hoje. Ficamos por aqui. Os scripts sempre no repositório recologia e é isso, até o próximo post. A, e se você curtir esse post, de uma olhada num blog chamado Freakonometrics, foi la que eu vi pela primeira vez essa ideia do Snake and Ladders

#Codigo original
#http://freakonometrics.blog.free.fr/index.php?post/2011/12/20/Basic-on-Markov-Chain-(for-parents)
 
#Matrix de transição
n=100
M=matrix(0,n+1,n+1+6)
rownames(M)=0:n
colnames(M)=0:(n+6)
 
for(i in 1:6) {
    diag(M[,(i+1):(i+1+n)])=1/6
}
 
M[,n+1]=apply(M[,(n+1):(n+1+6)],1,sum)
M=M[,1:(n+1)]
starting=c(4,9,17,20,28,40,51,54,62,64,63,71,93,95,92)
ending  =c(14,31,7,38,84,59,67,34,19,60,81,91,73,75,78)
 
for(i in 1:length(starting)) {
    v=M[,starting[i]+1]
    ind=which(v>0)
    M[ind,starting[i]+1]=0
    M[ind,ending[i]+1]=M[ind,ending[i]+1]+v[ind]
}
 
#Olhando o que criamos
nrow(M)
M[1,]
M[2,]
 
sum(M[2,]>0)
sum(M[10,]>0)
sum(M[50,]>0)
 
#Função para multiplicar a matriz de transição
powermat=function(P,h) {
    Ph=P
    if(h>1) {
        for(k in 2:h) {
            Ph=Ph%*%P
        }
    }
    return(Ph)
}
 
#Função para plotar o calculo das posições
initial=c(1,rep(0,n))
COLOR=rev(heat.colors(101))
u=1:sqrt(n)
boxes=data.frame(index=1:n,ord=rep(u,each=sqrt(n)),abs=rep(c(u,rev(u)),sqrt(n)/2))
 
#
position=function(h){
    D=initial%*%powermat(M,h)
    plot(0:10,0:10,col="white",axes=FALSE,xlab="",ylab="",main=paste("Posição após",h,"jogadas"))
 
    for(i in 1:n) {
        polygon(boxes$abs[i]-c(0,0,1,1),boxes$ord[i]-c(0,1,1,0),col=COLOR[min(1+trunc(500*D[i+1]),101)],border=NA)
    }
 
    segments(c(0,10),rep(0,2),c(0,10),rep(10,2))
    segments(rep(0,2),c(0,10),rep(10,2),c(0,10))
    segments(0:10,rep(0,11),0:10,rep(10,11))
    segments(rep(0,11),0:10,rep(10,11),0:10)
 
    for(i in 1:length(starting)) {
        arrows(x0=boxes$abs[starting[i]]-0.5,
               y0=boxes$ord[starting[i]]-0.5,
               x1=boxes$abs[ending[i]]-0.5,
               y1 =boxes$ord[ending[i]]-0.5,
               lty=3,length=0.10,col="darkgray",lwd=2)
    }
    text(boxes$abs-.5,boxes$ord-.5,boxes$index,cex=.7)
}
 
#figuras
 
position(1)
position(2)
position(10)
 
#Se quiser ver uma animação olhe esse codigo, no Sys.sleep da para configurar
#uma pausa entre os plots, deixei em 1/4 de segundo por jogada.
#for (i in 1:100) {
#                  position(i)
#                  Sys.sleep(0.25)
#                  }
 
 
#Probabilidade parcial
h=10
posição<-(initial%*%powermat(M,h))[59:61]/sum((initial%*%powermat(M,h))[59:61])
sum(posição)
posição
 
#Distribuição de jogadas
distrib<-initial%*%M
game<-rep(NA,1000)
 
for(h in 1:length(game)){
    game[h]<-distrib[n+1]
    distrib<-distrib%*%M
}
 
#Figuras de distribuição
plot(1-game[1:200],type="l",lwd=2,col="red",ylab="Probabilidade de ainda estar jogando",xlab="Jogadas")
position(50)

Forrageamento ótimo: Polinizadores!

Ao falar de polinizadores, agora eu lembro dessa imagem.

av0dbmd_700b

Isso porque polinização é talvez o melhor exemplo de serviço ambiental prestado.
Serviço ambiental é quando você ta recebendo um benefício do ecossistema em volta de você, e podemos até quantificar isso em unidade de dinheiro, reais ou dólar se preferir.

É isso ai, é só pensar de forma simples, sem abelhas, muitas culturas produzem menos sementes, ou até nenhuma semente. E ninguém quer pagar alguém pra ficar pegando pólen em uma plantinha e levando para outra. Nem temos muitas maneiras artificiais eficientes de fazer isso, pensando numa escala de produção de sementes em massa. Então polinização é de maior interesse para todos, mesmo que você não acha de primeira.

Bem todos os organismos forrageiam de forma a manter o seu sucesso maximizado, a gente vem vendo isso aqui e aqui, porque se você não faz isso, você tem menos energia para reprodução ou vai morrer no meio da sua vida e como vimos aqui, quem segue a mesma estratégia falha que você vai rodar. Agora podemos pensar em forrageamento ótimo não somente da perspectiva de comida, por exemplo flores forrageiam polinizadores para conseguir produzir sementes, e otimizam estratégias para isso, já que elas querem atrair polinizadores, mas não predadores.

Do outro lado, a gente pode pensar também numa colonia de abelhas, que de modo geral consiste de uma rainha e trabalhadoras em diferente castas com serviços diversos na colônia. Uma abelha forrageadora é uma trabalhadora com uma tarefa simples, coletar pólen e néctar. A eficiência dessa nossa amiguinha e a quantidade de comida que ela coleta é decisiva para o sucesso reprodutivo da colônia. O pólen e néctar é trazido de volta para a colônia, e guardado como reserva energética para a colônia. Alias, mel é um dos alimentos mais fáceis de estocar, não precisa de geladeira, não funga e até em piramides, encontraram mel intacto depois de mais de mil anos estocado, se a discovery não mentiu pra mim. Mas todos os processos metabólicos, incluindo crescimento, sobrevivência e taxa de reprodução são dependentes da quantidade de néctar e pólen estocados. Assim a seleção natural deve favorecer estratégias que maximizam da taxa líquida de ganho de energia enquanto forrageando atrás de néctar.

Vamos pensar num grupo de abelhas que tem pequenas colonias, pra facilitar o entendimento aqui.
As mamangaba, ou bumblebee em inglês. O nome do robozinho do transformes que era um fusquinha no desenho mas virou camaro nos filmes, já que a Alemanha não mandou dindin pro filme, mas whatever é uma abelha que estabelece colonias em para algumas de suas espécies. Vamos ver um modelinho tentando ver como funciona a interação dessa abelhas mamangaba forrageando no tremoceiro.

No tremoceiro, a flor la de baixo da inflorescência é mais velha, abre primeiro e contém normalmente mais néctar. Conforme o tempo passa, a gente corre na linha do tempo da fenologia da espécie, a flor fica velha demais e “morre”, vai virar fruto e a flor em cima dela passa a ser a última, até acabar as flores. Mas em um dia durante a fase da floração, há uma relação linear entre a posição da flor e o conteúdo calórico dessa flor, para o nosso tremoceiro aqui, a gente pode descrever essa relação da seguinte forma.

{C_i}={Calorias\ por\ flor} - {Inclinacao_i}

Aqui o índice i é referente a posição da flor na inflorescência, 1 é a flor mais de baixo, e se por exemplo a gente tem 20 flores, 20 é a la de cima, veja que a inclinação é negativa, isso acaba por impor um limite no número de flores possíveis, isso é, num teremos uma flor com quantidade negativa de néctar, uma flor não toma néctar da abelhinha, mas é só ser sensato ao usar essa relação, sempre lembrar dos limites dela.

Podemos representar essa relação assim.

01

A nossa mamangaba tem o comportamento de primeiro ir na flor de baixo, que contem mais néctar ou mair retorno em calorias, e vai indo na de cima, depois outra, até o topo da haste da inflorescência. Mas normalmente ela deixa a planta antes de visitar todas as flores daquela inflorescência.

Agora suponha que podemos coletar informação sobre a o consumo de energia de uma espécie de mamamgaba em um dia em particular, em um determinado local.

Matemagicamente, nos podemos modelar a eficiência energética baseado nas informações coletadas. Assim podemos estudar o movimento das mamangabas de uma planta para a próxima no ato de coletar néctar e dai analisar se ela está forrageando de uma maneira ótima ou não. O modelo matemático para calcular a taxa liquida de ganho energético para cada flor vai ser dada pela seguinte formula.

{Net_i} = {{P \sum \limits_i^n (\alpha+\beta \cdot i) - eb - P \cdot ew \cdot ( i - 1) - P \cdot i \cdot ef - ee \cdot (1 - P)} \over{tb + P \cdot tw \cdot (i - 1) + P \cdot i \cdot tf + te \cdot (1-P)}}

Onde

Net_i é o total de energia ganho visitando até a flor i

P é a probabilidade da inflorescência não ter sido drenada de néctar naquele momento, estar vazia.

eb é a energia gasta voando de uma inflorescência para outra

tb é o tempo gasto voando de uma inflorescência para outra

ew é a energia gasta voando de uma flor para outra na mesma inflorescência

tw é o tempo gasto voando de uma flor para outra na mesma inflorescência

ef é a energia gasta coletando todo o néctar da flor.

tf é o tempo gasto coletando todo o néctar da flor

ee é a energia gasta experimentando se a flor está vazia

te é o tempo gasto experimentando se a flor está vazia

n é a posição da ultima flor visitada pela abelha

O numerador é o total de energia ganha pela visita até a flor i e o denominador da formula descreve o tempo gasto forrageando na planta e movendo para a próxima planta.

Quando a razão é maximizada, nos podemos dizer que a abelha no seu forrageamento ótimo.

Vamos implementar essa formula como uma função no R

bumb<-function(P,eb,ew,ef,ee,tb,tw,tf,te,n,alpha,beta) {
    netE<-vector()
    sumi=0
    for(i in 1:n) {
        sumi=sumi+(alpha+beta*i)
        net=(P*sumi-eb-P*ew*(i-1)-P*i*ef-ee*(1-P))/(tb+P*tw*(i-1)+P*i*tf+te*(1-P))
        netE[i]<-net
    }
    return(netE)
}

Nessa função, entramos com os parâmetros das espécie vegetal e a saída é o ganho energético por visitar até a i^th flor, ou seja, da for 1 até a i.

Ai para ver o forrageamento ótimo, é só ver quando a abelha deveria “parar”.

02

Temos ali o ponto critico, se a nossa amiguinha mamangava abandonar a inflorescência antes, ela poderia ter saído daquela inflorescência com um melhor retorno energético, se ela sair depois daquele ponto, ela poderia ter tido um ganho maior, mais ao ficar, o gasto com deslocamento para a próxima flor, para inspecionar a próxima flor foi maior que o retorno, então começou a dar prejuízo energético continuar naquela flor. O ponto indicado é exatamente a hora de parar para ter o maior lucro com o menor investimento, energético.

Mas estamos olhando o forrageamento para um set fechado de parâmetros, e se esses parâmetros fossem um pouco diferentes?

Vamos pensar primeiro na probabilidade da inflorescência não ter sido drenada de néctar naquele momento, estar vazia.
Se essa probabilidade diminui, ou seja, a chance é maior de a flor estar vazia, isso vai diminuir nosso retorno energético como um todo não? Só imaginar o quanto isso vai resultar em encontrar plantas que não tem néctar, mas o gasto energético para ir olhar já foi, essa energia não volta mais.

03

Note ai que o forrageamento ótimo se move para a direita, ou seja, começa a valer mais a pena pegar mais flores, mesmo elas tendo cada vez menos néctar, porque é melhor garantir do que remediar, já que temos grande incerteza quando a próxima inflorescência.

Agora vamos se a gente manter a probabilidade da inflorescência não ter sido drenada de néctar naquele momento constante e alterar a tempo gasto voando de uma inflorescência para outra. Imagine que a vegetação é mais densa, existem mais obstáculos o que gera a necessidade de mais manobras no voo.

04

O ponto de forrageamento ótimo não se altera, mas o ganho energético diminui conforme as plantas estão mais espalhadas no ambiente. Então aglomeração de planta rendem mais energia, o que vai atrair mais polinizadores, ja que esses buscam o melhor ganho. E se a gente pensar, mesmo em muitas espécies, mesmo quando elas não estão exatamente num sistema como esse, tende a estar assim, aglomeradas numa mancha. É manchas podem ser bem mais rentáveis energeticamente do que espalhar uniformemente no campo as flores.

Agora vamos voltar o tempo gasto voando entre inflorescências para o normal e vamos varia o quanto uma flor a mais de néctar do que a próxima, mais especificamente, vamos homogeneizar um pouco mais as flores, diminuir aquela relação linear que falamos la em cima, alterando o coeficiente de inclinação beta.

Bem se a próxima flor tem o mesmo tanto de néctar que a anterior, vamos continuar um pouquinho mais nessa planta, já que da para encher os bolsos, ou perninhas nessa caso, com mais néctar.

05

Veja que o beta é um valor negativo, mas quanto mais reto, mais próximo de zero, onde zero seria todas as flores terem a mesma quantidade de néctar, mais vale a pena passar mais tempo numa inflorescência.

Ja o contrario.

06

É melhor sair antes, e quanto mais flores eu insisto em visitar, mais energia eu desperdiço trabalhando e menos eu levo para casa.

Agora a gente ainda pode fazer uma ultima perguntinha. E se a quantidade de néctar é imprevisível? Se é tudo na loucura e cada hora cada flor tem uma quantidade de néctar. Mesmo a planta fazendo o mesmo investimento energético que faria se fosse previsível a quantidade de néctar inflorescência.
A gente pode simular isso, pegando a média e desvio de conteúdo de néctar para a inflorescência com alta previsibilidade e fazer flores com conteúdo de néctar ao acaso, e jogar na nossa equaçãozinha e ver o que ocorre.

07

Cada linha cinza é uma simulação, uma determinada quantidade de néctar, e a linha vermelha é a média de 250 simulações, o caso médio para todas as simulações.

Olha ai que bolado, primeiro, vemos que o forrageamento ótimo quando é previsível a quantidade de néctar pela abelha, é bem menor para a abelha. Mas agora veja que no finalzinho, na décima flor, a curva de forrageamento ótimo estava subindo para flores com quantidade de néctar ao acaso, e passando as flores com quantidade de néctar previsível. La na frente nesse gráfico, o forrageamento em flores com quantidade de néctar ao acaso superaria até o retorno energético ótimo do forrageamento em flores previsíveis, pense que a próxima flor, se tiver próxima, sempre pode estar cheia até a boca do delicioso néctar, existe essa possibilidade, mas e ai, você teve que passar muito mais tempo forrageando, o que pode ser perigoso, você, uma abelhinha mamangaba pode ser comida no caminho, pode sei la ter um infarto de tanto trabalhar. Prevendo a quantidade de néctar de uma flor, a abelha pode trabalhar bem menos e já conseguir um bom retorno, sendo esse retorno seguro. Agora eu pergunto, como pode a abelha prever a quantidade de néctar?
Ai entra a comunicação entre espécies. Não é a toa que abelhas muitas vezes evoluem para serem especialista, se pensarmos em características das flores como sinais, coloração das pétalas, posição dessas, além de outras características podem funcionar como sinais de que ali tem néctar ou não, pensando por esse ponto de vista claro, existem outros. Além do que se a abelha tem “garantias”, ela pode preferir sempre a mesma espécie de planta, o que garante que o pólen sera transferido de forma mais eficiente, o que é interessante para a planta.

Veja que esse modelinho de polinização é bem legal, e podemos fazer ainda muitas outras modificações pequenas para se adequar a outras situações, por exemplo, trocar inflorescências por flores únicas, adicionar uma distância média entre flores, sorteando valores ao invés de fixo, se as plantas estão ao acaso no campo, ou modificar um pouco a ideia de inflorescência para manchas, e assim vai.

Bem ficamos por aqui, o script está no repositório recologia, polinização é bem legal, mas temos mais um post ainda sobre forrageamento ótimo, os predadores senta-espera maximizando seu forrageamento.

Referência:
Bernstein, R. 2003 – Population Ecology – An Introduction to Computer Simulations. John Wiley & Sons, Ltd 159pp.

Heinrich, B. 1983 – Do bumblebees forage optimally, and does it matter? American Zoologist 23: 273-281.

bumb<-function(P,eb,ew,ef,ee,tb,tw,tf,te,n,alpha,beta) {
    netE<-vector()
    sumi=0
    for(i in 1:n) {
        sumi=sumi+(alpha+beta*i)
        net=(P*sumi-eb-P*ew*(i-1)-P*i*ef-ee*(1-P))/(tb+P*tw*(i-1)+P*i*tf+te*(1-P))
        netE[i]<-net
    }
    return(netE)
}
 
#Calorias
i<-seq(0,10,1)
calories<-20-1.7*i
 
#figura 1
plot(i,calories,type="l",frame=F,ylab="Número de calorias da flor",xlab="Sequência de flores em um ramo da planta",
     main="Relação linear entre posição da flor e conteúdo calorico da flor",xlim=c(0,10),ylim=c(0,20))
 
#figura 2
saida<-bumb(P=0.6,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
plot(1:10,saida,type="l",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, tb=5, tw=0.05, tf=15, te=9, n=10, alpha=20, beta=-1.7, P=0.4",font.main=1,cex.main=0.9)
abline(h=max(saida),v=which(saida==max(saida)),lty=3,cex=2)
text(which(saida==max(saida))+0.5,max(saida)+0.025,"Forrageamento ótimo")
 
#Aumentando a chance da flor ter sido exaurida antes da visita
saida<-bumb(P=0.6,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
plot(1:10,saida,type="l",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, tb=5, tw=0.05, tf=15, te=9, n=10, alpha=20, beta=-1.7",
      font.main=1,cex.main=0.9)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=2,lty=2,lwd=2)
saida<-bumb(P=0.2,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=3,lty=2,lwd=2)
saida<-bumb(P=0.1,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=4,lty=2,lwd=2)
legend("topright",legend=c("P=0.6","P=0.4","P=0.2","P=0.1"),lty=2,col=1:4,lwd=2,bty="n")
 
#Aumentando distancia entre flores
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
plot(1:10,saida,type="l",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, P=0.4, tw=0.05, tf=15, te=9, n=10, alpha=20, beta=-1.7",
      font.main=1,cex.main=0.9)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=7,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=2,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=9,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=3,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=11,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=4,lty=2,lwd=2)
legend("topright",legend=c("tb=5","tb=7","tb=9","tb=11"),lty=2,col=1:4,lwd=2,bty="n")
 
#Flores maturam num intervalo de tempo menor
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
plot(1:10,saida,type="l",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, P=0.4, tw=0.05, tf=15, te=9, n=10, alpha=20, tb=5",
      font.main=1,cex.main=0.9)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.5)
points(1:10,saida,type="l",col=2,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.3)
points(1:10,saida,type="l",col=3,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.0)
points(1:10,saida,type="l",col=4,lty=2,lwd=2)
legend("topright",legend=c("beta=-1.7","beta=-1.5","beta=-1.3","beta=-1.0"),lty=2,col=1:4,lwd=2,bty="n")
 
#Flores maturam num intervalo de tempo menor
jpeg("06.jpg")
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
plot(1:10,saida,type="l",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, P=0.4, tw=0.05, tf=15, te=9, n=10, alpha=20, tb=5",
      font.main=1,cex.main=0.9)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.9)
points(1:10,saida,type="l",col=2,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-2.1)
points(1:10,saida,type="l",col=3,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-2.3)
points(1:10,saida,type="l",col=4,lty=2,lwd=2)
legend("topright",legend=c("beta=-1.7","beta=-1.9","beta=-2.1","beta=-2.3"),lty=2,col=1:4,lwd=2,bty="n")
 
#Experimento
saidaoriginal<-bumb(P=0.6,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
 
 
bumb<-function(P,eb,ew,ef,ee,tb,tw,tf,te,n,media,desvio) {
    netE<-vector()
    sumi=0
    for(i in 1:n) {
        sumi=sumi+rnorm(1,media,desvio)
        net=(P*sumi-eb-P*ew*(i-1)-P*i*ef-ee*(1-P))/(tb+P*tw*(i-1)+P*i*tf+te*(1-P))
        netE[i]<-net
    }
    return(netE)
}
 
#figura 7
plot(1:10,saida,type="n",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, P=0.4, tw=0.05, tf=15, te=9, n=10, tb=5",font.main=1,cex.main=0.9)
 
matriz<-matrix(NA,ncol=10,nrow=250)
for(i in 1:250) {
    saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,media=mean(calories),desvio=sd(calories))
    matriz[i,]<-saida
    points(1:10,saida,type="l",col="gray",lty=3,lwd=0.7)
}
 
points(1:10,colSums(matriz)/250,type="l",col="red",lty=3,lwd=3)
points(1:10,saidaoriginal,type="l",col="black",lty=3,lwd=3)
legend("bottomright",legend=c("Conteúdo de néctar conhecido","Conteúdo de néctar ao acaso"),lwd=3,lty=3,
       col=c("black","red"))

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)