Diferença entre estatística Bayesiana e Frequentista

Se você ainda não sacou a diferença entre os Paradigmas, segue um comentário que achei legal.

A Bayesian is one who, vaguely expecting a horse, and catching a glimpse of a donkey, strongly believes he has seen a mule.

 

Em português eu acho que ficaria algo como:

Um adepto da estatística Bayesiana é aquele que, vagamente esperando um cavalo, e vislumbrando um burro, acredita fortemente que viu uma mula.

 

Eu vi essa frase enquanto olhava o twitter do John Cook, mais especificamente a conta Stat Fact.

Mas não demorou muito e ele recebeu uma resposta:

A Frequentist, hoping to tell between a horse and a donkey, and glimpsing a mule, concludes he may as well have seen a horse.

Que continuando com minha tradução se utilizando de meu inglês nórdico é:

Um adepto da estatística frequentista é aquele que, querendo diferenciar entre um cavalo e um burro, ao vislumbrar uma mula, pode também concluir que viu um cavalo.

Convergência da cadeia em estatística Bayesiana.

Bem a gente fica falando de estatística Bayesiana, mas ainda não comentamos nenhum jeito de olhar se tudo esta funcionando.

Quando a gente olhou o que são cadeias de markov aqui e aqui, nos vimos que elas tem uma propriedade, que é estabilizar, mas como temos que começar de algum ponto, precisamos continuar acompanhando ela tempo suficiente para ela estabilizar. Agora existem formar de olhar isso, se a cadeia estabilizou?

Bem primeiro vamos olhar denovo aquele modelinho de “Anova Bayesiana” que fizemos a alguns post atrás, mas simplificando ele e deixando somente as estimativas de médias de cada grupo.

Primeiro geramos alguns dados e vamos relembrar como eles são com o seguinte gráfico.

01

Daqui para frente, tudo é bem igual o último post da anova. Rodamos o modelo no Openbugs, verificamos as estimativas, tiramos conclusões e a vida segue.

Porém, a gente pode ser questionado, ou deveríamos nos questionar na verdade, será que ta dando certo? Será que a cadeia ta convergindo ou eu parei num estimativa ruim, e minhas conclusões não vão gerar problemas?

Então quando nos olhamos nossos resultados, nesse caso basicamente olhamos isso:

> print(out, dig = 3) Inference for Bugs model at "anova.txt", Current: 3 chains, each with 1200 iterations (first 200 discarded), n.thin = 2 Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha[1] 2.181 0.236 1.708 2.018 2.178 2.342 2.653 1.001 2700 alpha[2] 2.198 0.239 1.730 2.035 2.200 2.358 2.672 1.003 770 alpha[3] 3.282 0.242 2.823 3.119 3.279 3.448 3.756 1.002 1900 alpha[4] 4.462 0.237 3.989 4.303 4.465 4.623 4.913 1.001 3000 sigma 0.927 0.091 0.768 0.862 0.921 0.985 1.121 1.001 2900 deviance 159.585 3.250 155.300 157.200 158.900 161.300 167.800 1.001 3000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 4.9 and DIC = 164.5 DIC is an estimate of expected predictive error (lower deviance is better). >

Além de olhar as estimativas, a gente pode dar uma olhada naquele Rhat, “hat” de chapeuzinho em inglês, seria algo assim  \bar{r} , mas aqui as coisas são texto simples, então a sai escrito e não o nome, alias por isso ta alpha ali e não  \alpha , mas isso é só um detalhe que pode facilitar a ver formulas de artigos nessas saídas, porque as vezes fica tão distante o que a gente le em livros do que ve, principalmente em relação a esses simbolos.
Mas chega de divagar, então vimos que temos esse número ai, como principio básico, regra de ouro, rule of thumb, ele tem que ser menor que 1.1, quanto maior, quer dizer que a cadeia em questão não esta convergindo, quanto mais próximo de 1, quer dizer que a cadeia ta convergindo legal, note que ali em baixo ele explica o que é Rhat

For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

Que quer dizer, para cada parâmetro, n.eff é uma medida crua do número efetivo do tamanho da amostra (amostra da cadeia de markov) e Rhat é uma escala potencial de redução do fator.

Então primeiro passo é olhar esses valores, ver se todos são menor que 1.1:

> out$summary[,"Rhat"] alpha[1] alpha[2] alpha[3] alpha[4] sigma deviance 1.001234 1.003087 1.001563 1.000584 1.001182 1.000701 > all(out$summary[,"Rhat"] < 1.1) [1] TRUE

Mas lembre-se, é apenas um principio básico, não quer dizer que 1.100000001 é uma cadeia que não convergiu, não presta para nada e 1.099999999 é uma cadeia que convergiu.

Podemos usar um argumento do comando bugs, o “codaPkg=TRUE” e salvar as cadeias de markov inteiras para investigar mais.

Depois usamos o pacote coda do R que traz mais algumas coisas legais.

Agora a gente roda denovo o comando e ao olhar a saida:

> out<-bugs(data=bugs.data,inits=inits, parameters.to.save=params, + model.file="anova.txt",n.thin=nt,n.chains=nc,n.burnin=nb,n.iter=ni, + codaPkg=TRUE) out > [1] "/tmp/Rtmpb1FBh5/CODAchain1.txt" "/tmp/Rtmpb1FBh5/CODAchain2.txt" "/tmp/Rtmpb1FBh5/CODAchain3.txt"

A gente ve que não tem o output igual la em cima, mas ele salvou cadeia por cadeia em um arquivo, e nos da o local onde esses estão salvos, mas com o comando read.bugs() ta preparado para ler essa informação.

Ai depois de ler isso, a gente pode ver as cadeias inteiras:

02

Aqui a gente pode olhar se realmente a cadeia convergiu, que é ela não estar com uma tendencia, subindo descendo, dando volta ou sei la, coisas estranhas assim. Outra coisa é que a gente repetiu 3 cadeias, então todas as 3 deveriam ir ou tender para a mesma estimativa, se cada uma vai para um lado, seria um problema, mas aqui parece que esta tudo ok.

Outro coisa que podemos ver, é como são as densidades das estimativas.

La em cima a gente ta olhando:

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha[1] 2.181 0.236 1.708 2.018 2.178 2.342 2.653 1.001 2700

Média, desvio e os quantiles, a gente devia olhar se realmente essas estimativas tem uma distribuição normal, como ela é e tal, ai fazemos mais um gráfico e vemos:

03

Que parece tudo ok, notem que é importante olhar se as cadeias estão indo na mesma direção, convergindo para a mesma coisa, como aqui ta praticamente uma em cima da outra, isso é um bom sinal.

Agora a gente olhar se existe alguma tendencia temporal

04

E não é temporal de dias, temos uma cadeia que tem uma sequência, valor 1, valor 2, valor 3. Aqui a gente olha a auto-correlação, ou o fato do valor 1 ser muito parecido com o 2, o 2 com o 3, se isso ocorrer, logicamente isso vai influenciar na nossas medidas finais, pense que o isso deixaria por exemplo o gráfico acima de densidades mais gordinho, se eu não estou imaginando errado. Então verificamos aqui se existe uma tendência, de subir, descer, fazer curva. So lembrar que isso é parecido com o que fizemos com os dados do John Snow aqui, so que ao invés de autocorrelação espacial, aqui é temporal, temporal no sentido da ordem das amostragens. Se existe problemas aqui, uma coisa é alterar o Thinning la em cima, que é ir pulando de n em n valores e dai pegar uma amostra, para evitar essa auto-correlação. Agora da para entender porque existe aquilo la em cima.

E para quem gosta de teste, a gente pode fazer o “Gelman-Rubin convergence diagnostic” com o comando gelma.diag. O shrink factor tem que ser menor que 1.05 (Regra de ouro, não guie sua vida por esses números).

> gelman.diag(out.coda) Potential scale reduction factors: Point est. Upper C.I. alpha[1] 1 1.00 alpha[2] 1 1.02 alpha[3] 1 1.01 alpha[4] 1 1.00 deviance 1 1.00 sigma 1 1.02 Multivariate psrf 1.01

Certo a gente tem uma estimativa por cadeia, e uma multivariada para o modelo todo. Podemos olhar esse shink factor num grafico também aqui.

05

Certo, eu não sei vocês, mas eu não entendi nada desse treco de shink factor. Mas deve ser importante porque tem até uma pagina no wikipedia sobre shink factor . Mas se olharmos a documentação da função do Gelman da uma luz. Basicamente, como a gente viu nos post de markov chain e hidden markov chain, a gente começa em um ponto, aqui até a gente gera valores aleatórios para iniciar a cadeia, e ai vamos seguindo na cadeia, logo as estimativas no começo pode ser ruim, e vão variar bastante, mas mais para frente elas vão estabilizando. Shink factor é exatamente isso, shrink é encolher em português, então se a variação do valor ta encolhendo cadeia a frente ou não, encolhendo a variação, encolher no sentido que no começo de cada cadeia era um valor para cada cadeia, mas no final todas as cadeias tão concluindo o mesmo valor. Um dos jeitos de calcular isso para ter uma ideia é exatamente algo similar a Razão F, como a gente viu aqui na estatística F.

Assim temos:

{\bar{\sigma}^2 }={(n-1) \cdot {W \over n} + {B \over n}}

Onde n é o número de iterações, o tamanho da cadeia, W é a variação de todas as cadeias combinadas e B é a variação de cada cadeia. Igual a gente fazia com a distribuição F.
Se você olhar o help da função você vai ver que ele da uma melhorada nessa idéia, já que ela não funciona para alguns casos, mas por enquanto vamos ficar com ela, apenas para ter um feeling que não é nada do outro mundo esse Rhat ou da onde ele sai.

O gráfico é essa idéia ao longo do tempo (tempo = iterações aqui), logo o alpha[3] ficou bem estranho não? Deveríamos verificar ele.

E por último, com a convergência garantida, podemos usar esses resultados para recuperar as estimativas e intervalos de confiança para cada cadeia.

> out.summary$stat Mean SD Naive SE Time-series SE alpha[1] 2.1807257 0.23620474 0.004312489 0.004229698 alpha[2] 2.1981760 0.23868415 0.004357756 0.004353213 alpha[3] 3.2820577 0.24153206 0.004409752 0.004408873 alpha[4] 4.4619603 0.23705762 0.004328060 0.004329306 deviance 159.5846000 3.24957681 0.059328884 0.062254119 sigma 0.9270018 0.09066346 0.001655281 0.001714132 > out.summary$q 2.5% 97.5% alpha[1] 1.707950 2.653050 alpha[2] 1.729975 2.672025 alpha[3] 2.823000 3.756025 alpha[4] 3.988900 4.913025 deviance 155.300000 167.800000 sigma 0.768390 1.121025

Que nada mais é que o output quando o “codaPkg=TRUE” não esta ligado, ou o “codaPkg=FALSE”, que é o default do comando bugs, salvo o arredondamento.

Bem é isso ai, o script esta aqui e no repositório recologia do github também.
Se eu escrive alguma bobeira é so avisar pelos comentários 🙂

#Criando dados, abrindo pacotes, etc
library(R2OpenBUGS)
set.seed(15)
dose<-factor(rep(c("Controle","Dose1","Dose2","Dose3"),each=15))
medidas<-rnorm(60,rep(c(2,2,3,4),each=15))
 
#Grafico dos dados
boxplot(medidas~dose,frame.plot=F,ylim=c(0,6),ylab="Medidas",xlab="Dose")
 
# Escrevendo o modelo para o OpenBugs
sink("anova.txt")
cat("
model {
  # Priors
  for (i in 1:4){
    alpha[i] ~ dnorm(0, 0.001)
    }
  sigma ~ dunif(0,100)
  # Likelihood
  for (i in 1:60) {
    medidas[i] ~ dnorm(mean[i], tau)
    mean[i] <- alpha[dose[i]]
    }
  # Quantidades Derivadas
  tau <- 1 / ( sigma * sigma)
  }
  "
,fill=TRUE)
sink()
 
# Juntando os dados.
bugs.data <- list(medidas=medidas,dose=as.numeric(dose))
# Gerador alearorio de valores iniciais
inits <- function(){ list(alpha = rnorm(4, mean = mean(medidas)), sigma = rlnorm(1) )}
# Parametros a serem estimados
params <- c("alpha","sigma")
# Configurações do MCMC
ni <- 1200
nb <- 200
nt <- 2
nc <- 3
 
# Começando o Gibbs MCMC
out<-bugs(data=bugs.data,inits=inits, parameters.to.save=params,
model.file="anova.txt",n.thin=nt,n.chains=nc,n.burnin=nb,n.iter=ni)
# Olhando o resultado, restringindo a 3 digitos
print(out, dig = 3)
 
#Estatistica Rhat
out$summary[,"Rhat"]
 
#Verificando se todos ficaram menor que 1.1
all(out$summary[,"Rhat"] < 1.1)
 
#Salvando as cadeias inteiras: codaPkg=TRUE
out<-bugs(data=bugs.data,inits=inits, parameters.to.save=params,
model.file="anova.txt",n.thin=nt,n.chains=nc,n.burnin=nb,n.iter=ni,
          codaPkg=TRUE)
 
out
out.coda <- read.bugs(out)
 
#Abrindo o pacote coda
library(coda)
 
#Olhando as cadeias inteiras
xyplot(out.coda)
 
#Densidade dos parametros estimados
densityplot(out.coda)
 
#Procurando autocorrelação
acfplot(out.coda)
 
#Shrinkage factor
gelman.diag(out.coda)
gelman.plot(out.coda)
 
#Retirando estimativas dos parametros
out.summary <- summary(out.coda, q=c(0.025, 0.975))
out.summary$stat
out.summary$q

Crescimento logístico e estabilidade vs instabilidade

Algumas vezes eu vi noticias sobre infestações de caramujos por aqui no meu estado, noticias como essa aqui e essa aqui.

Ai eu pensava, poxa uma população que tem hora que aparece, depois some, talvez sumir porque todo mundo começa a matar os caramujos, mas como ela chega a um pico assim, como não extingue?

Bem essa espécie dessas notícias na verdade é uma espécie introduzida, o caramujo africano.

Pensando na equação de crescimento logística, é claro que a taxa de crescimento pode depender da onde a espécie esta. Como esta veio da africa, chegou num lugar onde não tem predadores nativos, especializados em devorar pobres caramujinhos indefesos, podemos imaginar que ela consiga ter uma taxa de crescimento maior que do lugar da onde ela veio originalmente.

Certo, então lembrando a formulinha que estamos usando para crescimento logístico, se não entender vale a pena dar uma olhada nesse post aqui, aqui e aqui.

  N_{t+1} = r\cdot N_{t} \cdot (1- \alpha * N_{t} )

Certo então o que estamos falando acima é o que acontece conforme o r (taxa de crescimento) fica maior e maior e maior? Sera que tudo tende a estabilidade como aqui?

Então criamos varias populações com valores de taxa de crescimento (r) diferentes e colocamos num gráfico e olhamos, tentativa e erro 🙂

01

Olha ai, populações com valor do r de 1.3 e 1.6 são bem o que a gente ve em todo santo grafiquinho de crescimento populacional correto?
Mas quando a taxa de crescimento (r) ficou 1.9 ou maior, ela parece que ficou oscilando, ja não ficou aquela linha reta e estável na capacidade suporte do ambiente (k, ou nesse caso 1/alpha que é a mesma coisa que k) como normalmente ocorre.

Vamos fazer o seguinte, vamos desmembrar essas várias populações cada uma em um gráfico diferente, usando a função xyplot do pacote lattice do R a gente consegue fazer facilmente esse gráfico por população, ai vamos olhar ele denovo.

02

O que ta acontecendo é que quando o r esta chegando a um certo valor, não esta mais tendo uma “estabilidade”, um único ponto de estabilidade quando o tamanho populacional esta chegando perto da capacidade suporte, ela parece que fica oscilando entre vários pontos, vários tamanhos populacionais.

Mas como isso é possível?

Vamos olhar a formulinha do crescimento logístico:

  N_{t+1} = r\cdot N_{t} \cdot (1- \alpha * N_{t} )

Veja so,   N_{t+1} depende do tamanho populacional no tempo anterior. Certo, mas a taxa de crescimento (r) só é afetada pela população que existe (1- \alpha * N_{t}), tudo o que ta dentro do parenteses é o que diminui o r, ja discutimos isso. Acontece que se a taxa de crescimento é alta o suficiente, vamos supor uns números tolos, a capacidade suporte é 100, e a população esta em 99 indivíduos, se a taxa de crescimento é muito alta, a população pula de 99 para 119, muito acima da capacidade suporte, ai vai rolar uma alta superpopulação, e ela vai ter uma queda brusca, o r vai cair para um número negativo, a população vai cair para menos de 100, mas depois disso ela tem uma subida acima da capacidade suporte denovo porque ela cresce muito rápido, e nisso ela vai fica oscilando e não vai simplesmente estabilizar e ficar constante no tamanho populacional da capacidade suporte do ambiente. Esse é o caos, é essa oscilação que a gente ta vendo.
Aqui é um caso onde ao invés de pensar em capacidade suporte (k), é mais fácil pensar no impacto que um indivíduo exerce no outro (alpha ou \alpha), ou seja, o ambiente enche rápido demais ai o impacto da superpopulação não simplesmente baixa o tamanho populacional para a capacidade suporte do ambiente, baixa abaixo disso, o quanto ela abaixa depende do tamanho populacional atual e do alpha(\alpha), mas uma vez abaixo da capacidade suporte, ta aberta a possibilidade de subir denovo o tamanho da população, inclusive acima da capacidade suporte denovo.

Certo, agora ja ta fazendo sentido biológico essa oscilação ai na função, beleza. Mas o quanto ela oscila?

Novamente partimos para a força bruta e fazemos muitas e muitas populações, com vários e vários valores de taxa de crescimento (r) diferentes, e vemos o que acontece, vemos que valores, depois da capacidade suporte, o tamanho populacional pode assumir, basicamente entre quais pontos aquela oscilação que vamos ali em cima vai estar rolando.

03

Certo, agora fizemos esse gráfico para todas as populações e vamos olhar 3 pontos que estão marcados com uma linha colorida para entender melhor.
Primeiro o risco vermelho, quando a taxa de crescimento é 1.8, simplesmente a população cresce e para na capacidade suporte, e la vai estacionar até ser perturbada.
Aqui conforme vamos da esquerda para a direita do eixo x, estamos olhando taxas de crescimento maiores, agora no risco verde temos uma taxa de crescimento de 2.4. Aqui a população ja fica oscilando entre depois tamanhos populacionais, que não são a capacidade suporte, um acima e um abaixo.
Agora vamos mais para frente ainda e olhar uma taxa de crescimento 2.8, aqui ja virou uma zona, igual quando olharmos no gráfico la em cima, a população fica oscilando la em entre vários tamanhos, que nem uma louca, podemos dizer que ela é bem instável.

Outro ponto ainda, é que se começarmos a populações com tamanhos diferentes, veja o que acontece no gráfico abaixo.

04

Agora um último ponto importante, é que fica tão difícil prever, que dependendo do tamanho da população que começamos a observar, se nossas observações começam com tamanhos de população diferentes, ela oscila de forma diferente, ela não tende a oscilar de um jeito e fica assim para sempre. Ou seja dependendo de qual o tamanho populacional, ela vai ficar variando de um jeito.

Isso torna bem difícil qualquer tipo de precisão, e impõem difíceis limites para o uso de métodos como esse aqui, que discutimos para estimar o crescimento populacional, se usarmos o crescimento logístico ao invés do crescimento populacional ai, temos que ter em mente que é mais complicado se a população tiver uma taxa de crescimento maior que dois mais ou menos.

E agora o que isso tem haver com os caramujos, desde a primeira vez que vi aquelas notícias, e ja tinha visto essa idéia discutido nesse post, o que acontece quando a taxa de crescimento é muito elevada, eu ficava pensando se a população desses caramujos simplesmente não tem uma taxa de crescimento assim, gigantesca, só que quando há um pico populacional muito alto, ai ele vira uma “praga”, incomoda a população. Mas devido a população ficar grande demais, ele mesmo faz com que a própria população diminua, sei la, deve faltar buraco, comida, de tudo um pouco para os caramujinhos, ai a população vai la para baixo, enquanto a população esta la em baixo, ou esta oscilando entre números baixos, ninguém mais liga para ele, ele não vira mais “notícia”, ai derrepente a população chega num momento que da um surto novamente, e pimba, ele vira “notícia” novamente.

É claro que não tenho certeza que o caso dos caramujos entra aqui, nos casos de taxa de crescimento muito alta, mas imagino que essa é uma boa hipótese de porque a população desse caramujo vira e mexe vira uma “praga” e depois some do mapa por um tempo.

Bem galera é isso. So mais uma coisa legal, depois de tanto ver as pessoas falando de github no R-Bloggers, eu resolvi criar o meu para tentar aprender como funciona.
A partir de agora eu pretendo deixar os scripts la no repositorio chamado recologia.
Eu tenho mais alguns repositórios também no qual vou ir pondo tranqueiras, mas de qualquer forma, esse script abaixo vai estar la disponível também se alguém quiser.

O endereço do reposito é: https://github.com/Squiercg/

Agora sim é isso 🙂

Referencias:

Stevens, M. H. – 2011 – A Primer of Ecology with R. Srpinger

###################################################
### Estabilidade vs Instabilidade
###################################################
rm(list=ls())
 
#Crescimento Logístico
dlogistic <- function(alpha=0.01, rd=1, N0=2, t=15) {
  N <- c(N0, numeric(t))
  for(i in 1:t) {
    N[i+1] <- N[i] + rd * N[i] * (1 - alpha * N[i])
  }
  return(N)
}
 
#Criamos um vetor com varia possibilidade de taxa de crescimento
#entre 1.3 e 2.8
rd.v <- seq(1.3, 2.8, by=.3)
t <- 15
 
#Usamos nossa função para ver como o tamanho populacional ao
#longo do tempo
Ns <- data.frame(sapply(rd.v, function(r) dlogistic(rd=r, t=t) ))
 
#Então fazemos um grafico:
matplot(0:t,Ns,type="l",lty=1,xlab="Tempo",ylab="Tamanho populacional",frame=F)
legend("bottomright",legend=rd.v,lty=1,col=1:length(rd.v),bty="n")
 
#Reorganizando os dados
tmp <- data.frame(rd=as.factor(rd.v), t(Ns))
Ns2 <- reshape(tmp, varying=list(2:ncol(tmp)), idvar="rd", v.names="N",direction="long")
 
#Refazendo o gráfico de forma diferente
library(lattice)
xyplot(N ~ time|rd, data=Ns2, type="l", layout=c(3,2,1), col=1,
       ylab="Tamanho Populacional",xlab="Tempo" )
 
#Pegando mais possibilidades de taxa de crescimento
#entre 1 e 3
num.rd <- 201
rd.s <- seq(1,3, length=num.rd)
t <- 400
 
#Usando a função logistica nessa varias possibilidades
tmp <- sapply(rd.s, function(r) dlogistic(rd=r, N0=99, t=t) ) 
 
#Reorganizando os dados
tmp.s <- stack( as.data.frame(tmp) )
names(tmp.s) <- c("N", "Old.Column.ID")
tmp.s$rd <- rep(rd.s, each=t+1)
tmp.s$time <- rep(0:t, num.rd)
N.bif <- subset(tmp.s, time > 0.5 * t )
 
#Graficos dos pontos de parada do tamanho populacional
plot(N ~ rd, data = N.bif, pch=".", xlab=quote("r"["d"]),frame=F)
abline(v=1.8,col="red",lty=2,lwd=2)
abline(v=2.4,col="green",lty=2,lwd=2)
abline(v=2.8,col="blue",lty=2,lwd=2)
legend("bottomleft",legend=c(1.8,2.4,2.8),lty=2,col=c("red","green","blue"),
       bty="n")
 
#Selecionando diferentes pontos iniciais da população
N.init <- c(97,98,99)
t <- 30
Ns <- sapply(N.init, function(n0) dlogistic(rd=2.7, N0=n0, t=t) )
 
matplot(0:t,Ns,type="l",lty=1,xlab="Tempo",ylab="Tamanho populacional",frame=F)
legend("bottomleft",legend=N.init,lty=1,col=1:length(N.init),bty="n",
       title="Tamanho Inicial")

Encontrando a máxima verosimilhança de um conjunto de dados usando um grid para a busca.

A boa e velha verosimilhança, ou likelihood em inglês, para quem preferir. Eu prefiro likelihood porque tem menos letras e não tem acento, eu sempre erro o acento da palavra verosimilhança, aliás eu erro a acentuação de um monte de palavras, talvez erre até um pouco na vida, infelizmente.

Mas a gente já falou a muito tempo atras de verossimilhança aqui. Aliás o que vamos fazer nesse post a gente já fez aqui também, quando falamos de aproximação por grid em estatística bayesiana, só que para a distribuição binomial.

Eu sempre imagino se sou só eu que demoro entender as coisas, ou tem mais gente lenta como eu.

Por exemplo, uma coisa que me perturbou muito tempo foi a diferenças entre  \mu (letra grega mu minúscula), que representa a média e porque as vezes tinha um  \bar{\mu} assim, com esse chapéuzinho. Se era a mesma coisa, se eram coisas diferentes.

Ai tinha a explicação que com chapéuzinho, com a barrinha em cima era a média das amostras, e sem chapeuzinho era a média da distribuição, e isso era estranho porque amostra vem dos dados, não deveria ser a mesma coisa?

Depois que eu vi nos livros do Mark Kéry a idéia de simular dados, para depois realizar uma análise de dados para ver se recuperávamos o padrão simulado que isso começou a fazer mais sentido, quem deu uma olhada nesse post aqui e nesse outro aqui viu que eu fiz um gráfico tentando resumir tudo, vou colocar ele denovo aqui em baixo.

02
03

A diferença é que quando eu simulo dados, eu tenho como colocar no gráfico as medidas “reais”, agora no segundo caso, quando usamos a análise em dados reais, não tem linha “real” nem nunca vai ter, podemos pensar que a natureza fez os dados, mas ela não mostra os parâmetros que usou, a gente estima e acredita que funcionou e para por ai. A gente nunca vai saber como diabos são os dados reais, e essa é a diferença entre  \mu e  \bar{\mu} , ou qualquer coisa com e sem chapeuzinho, qualquer parâmetro, a gente estima eles a partir de amostras, então a gente tem a medida com chapeuzinho, a medida sem chapeuzinho nunca saberemos ao lidar com dados reais, so que como as vezes a gente simula dados, simula amostras a partir de um  \mu e consegue ter um  \bar{\mu} parecido, a gente acredita que esse procedimento funciona e usa ele nos dados reais, mas não temos como ter certeza, nada garante. No caso da estatística bayesiana a gente ainda combina esse  \bar{\mu} com um prior, com algo que a gente acha para gerar uma distribuição posterior, que é um intervalo de valores que achamos que o  \mu ta la dentro, mas nada é exato.

Certo, um bocado de blablabla, mas vamos fazer algo legal, vamos estimar o likelihood usando um grid de valores de média ( \bar{\mu} ) e desvio ( \bar{\sigma} ).

Então a gente gera uma amostra a partir de uma distribuição usando a função rnorm do R que gera números ao acaso seguindo os parâmetros que estipulamos.

media <- 15
desvio <- 2.5
n <- 100
x<-rnorm(n,media,desvio)

Certo, agora imagine que esse passo foi feito por outro pessoa, que a gente não tem contato, a única coisa que a gente ve são os valores em x, ou seja os 100 valores que estão dentro de x.

Como sempre e todo o sempre, a primeira coisa que a gente faz é um gráfico para ver como são os dados que temos em mãos

01

Certo, so de bater o olho nos parece uma distribuição normal, mas porque?
Oras, porque tem esse jeitão em forma de sino, com bordinhas, ai um levantando no meio, uma barriguinha, acredito todo mundo deve ter uma imagem de sino na cabeça.

Mas então o primeiro passo é imaginar uma distribuição adequada. Uma idéia aqui também seria fazer uma tabelinha com as forma das distribuições mais conhecidas, e ir no olhômetro.

Mas supondo que estamos satisfeitos com a idéia de usar uma distribuição normal, precisamos ajustar os parâmetros da distribuição a usar, a distribuição normal é baseada em dois parâmetros, a média e o desvio (também conhecido como variância, precisão, depende do caso e do uso do freguês). Esses parâmetros, como falamos ali em cima, podem ser descritos comumente com as letras gregas  \mu e  \sigma (isso é um sigma minúsculo). Mas quando eu falo ajustar parâmetros eu falo que temos que dar valores para eles.

Mas vejamos, temos uma forma, a distribuição normal, agora precisamos encaixar ela naquele histograma la em cima.
Vamos ver o que acontece para diferentes valores desses parâmetros:

02

Olha ai, vejam que alguns casos tudo fica muito bonitinho, um encaixe bom, todas as barrinhas bem perto da linha, mas alguns casos, por exemplo média 18 e desvio 5, o encaixe é horrível, esses valores são muito ruins, mas com média 14 e desvio 3 até que a curvinha, a distribuição normal encaixa bonitinho nos dados. Mas aqui usamos 4 valores de média, 12, 14, 16, 18 e usamos 4 valores de desvio, 2, 3, 4 e 5. Ai fomos no olhômetro e ficamos tentando perceber qual a melhor combinação que ficava melhor.

Basicamente é isso que a gente quer, no entanto a gente poderia aumentar a precisão aumentando a quantidade de valores possíveis, ao invés de 12 a 18 pulando de 2 em 2, a gente podia pular de 1 em 1, de 0.5 em 0.5, ou menos, mas ai o problema é que conforme aumentamos esse grid de busca, aumentamos a quantidade de gráficos que temos que olhar, ou seja usando a técnica do olhômetro nos estamos limitados a quantos gráficos podemos olhar se o encaixe ficou bom ou não. E se aumentarmos a precisão então ferrou, porque as possibilidades crescem exponencialmente conforme aumentamos o número de desvios e médias que queremos olhar, ja que o número de casos é o número de médias vezes o número de desvios usado.

Agora a gente chegou ao ponto de o que diabos o likelihood ou verosimilhança faz, ela substitui o nosso “olhômetro”, ele ve se o encaixe ficou legal, a gente calcula um número que quanto maior, melhor o encaixe.

Algebrismos a parte, agora já deu para imaginar o procedimento, pegar o feeling da coisa, a gente vai escolher valores arbitrários, bastantes valores, claro que sendo razoável nessa escolha, não vou testar o caso de média um zilhão, porque pelo gráfico da para ver que não faria sentido um valor assim, então nem precisa desperdiçar esforço num número assim, mas escolhendo números razoáveis, numa precisão que de tempo para o computador processar, a gente vai fazer a curva, ver se ela encaixa legal, e registrar o valor do encaixe, ou o likelihood, depois é so ver qual é melhor.

Como a gente viu na figura anterior, como eu testo todas as combinações, eu posso registrar numa matriz todos os valores de todas as combinações possíveis de likelihood, e depois vejo o maior.

Notem que os valores sempre ficam negativos, mas o menos negativo (mais próximo de zero) é maior que um número muito muito grande mas negativo.

Depois de tudo isso eu posso olhar o resultado do maior valor na matrix, ou procurar direto o maior valor.

> max_lh<-max(log_lh,na.rm =T) > max_lh [1] -233.3122

E temos ai nosso melhor encaixe, é so olhar onde ele esta na matriz que temos o valor de média e desvio que deveríamos estar mais inclinados a acreditar que deva ser o valor real dos dados.

Mas eu sou sempre horrível para ver valores e matrizes, então eu faço gráficos.

Nesse caso a gente pode até usar uma analogia, a média e desvios reais são uma flecha que acertou o alvo em cheio, so que a gente não ve a flecha, a gente so sabe que ela acerta em cheio, o que a gente ta tentando construir é o alvo, a gente tenta desenha o alvo e espera que a flecha tenha acertado em cheio. Mas o alvo que a gente faz nunca é perfeito. As vezes pior ou melhor, dependendo das nossas escolhas para desenhar ele, mas olhando esse caso aqui específico ficou assim.

03

Aqui, olha ai como a gente sabe onde a flecha está, a gente pinta um pontinho cinza para ver se ta dando certo esse esquema, ai pintamos de vermelho a area mais provavel, e conforme a gente muda as cores do vermelho para o azul, a chance é menor (os likelihood são valores menores).

E a flechada acertou mesmo no alvo que a gente pintou, então a gente fez um bom trabalho, esse procedimento parece funcionar bem.

Se a gente compara o valor de média que ficou no melhor encaixe, o maior valor de likelihood, ou maximum likelihood, ou máxima verosimilhança, vemos o seguinte, em relação as possibilidades que testamos:

> mu_mle <- muhat[maxrow] > mu_mle [1] 14.37779 > mean(x) [1] 14.67122

Esse é nosso melhor chute de valor para média, ou o  \bar{\mu} , se fazemos um gráfico com todos os valores que usamos no grid para média, e colocamos esse como o separador de águas, menor que ele pintando de verde e maior pintando de azul vemos que:

04

Olha ai a média real ta quase junto do “divisor de águas”. Muito legal né.

Agora qual a relação disso com os modelos que fazemos todos os dias, a gente faz esse processo toda hora que chama o comando lm, ou glm para outras distribuições, ou muitos outros comandos do R.

Veja os valores que estimamos para média e desvio com nosso método, e compare com a boa e velha formulinha da média e do desvio padrão.

> mu_mle [1] 14.37779 > sd_mle [1] 2.58246

Agora note que são quase a mesma coisa que o valor de média e desvio usando as funções mean e sd.

> mean(x) [1] 14.67122 > sd(x) [1] 2.507243

Média e desvio são o corte de caminho para todo esse procedimento? Ai você pensa, então tudo isso é inútil?

Não, lembra la atras quando usamos um procedimento de grid para estatística bayesina? Isso é o procedimento geral, da para usar com qualquer distribuição, no mundo não existe somente a distribuição normal, existe muitas delas, na página do John Cook tem uma diagrama só para termos uma idéia das possibilidades de distribuições, ache a normal nesse bolo ai.

Agora a relação com quando a gente usa o comando lm do R é o seguinte. Como só temos um nível, o x, a gente pode usar ele para estimar a média falando na formula que so existe um nível (usando x~1).

Então usamos o comando:

summary(lm(x~1))

E recebemos a resposta:

Call: lm(formula = x ~ 1) Residuals: Min 1Q Median 3Q Max -6.2405 -1.7020 -0.5832 1.7584 5.8910 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 14.6712 0.2507 58.52 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.507 on 99 degrees of freedom

Agora olha que fascinante, notem que aqui:

Estimate Std. Error t value Pr(>|t|) (Intercept) 14.6712 0.2507 58.52 <2e-16 ***

Estimate é a média, não ta idêntico ao nosso valor porque aqui é usado cálculo infinitesimal, que é o mesmo que um grid infinitamente preciso, com infinitos números, mas como cálculo é complicado, começamos de um jeito mais fácil, com um grid, e temos um jeito agora de conferir se quando tentarmos isso com calculo estaremos indo no caminho certo, ou não, viu é muito útil seja para um fim como para outro.
Mas então, o estimate é o meio da bola vermelha do gráfico la em cima, e o Std. Error, ou erro, é a precisão com que podemos falar que esse valor ta certo. Lembre-se que na figura la em cima, o vermelho tem uma considerável área, a gente ta falando que qualquer valor ali ta valendo, ou seja, esse valor é relativo ao tamanho da bolota vermelha la naquele gráfico, se ele fosse menor, a bola seria menor, se ele for grande a bola é grande do local onde assumimos que a flecha acertou.

Ta certo, agora olha esse número aqui

Residual standard error: 2.507 on 99 degrees of freedom

Olha ai é o nosso desvio do modelo. Esses são os parâmetros do modelo que dão a forma aquela curva para encaixar nos dados.

Por último a gente ve se ficou centralizadinho, porque se não ficou centralizado, não encaixou legal, certinho, o que pode ocorrer se la em cima escolhemos o modelo errado, a distribuição errada.
Então diminuímos todos os valores de média dos dados e fazemos um histograma.

05

E se quiser fazer isso com o ajuste que o comando lm fez, é so usar o comando resid no modelo gerado pelo lm e fazer um histograma, dessa forma:

hist(resid(lm(x~1)))

Aqui eu to fazendo um histograma do resíduo do modelo produzido pelo lm, se achou complicado é so fazer parte a parte.

modelo<-lm(x~1)
residuos<-resid(modelo)
hist(residuos)

Ambos os jeitos vão fazer exatamente a mesma coisa, mas as vezes usar um comando dentro do outro pode ser confuso logo de primeira, até aprender as coisas.

Que legal né, agora a gente não precisa olhar os valores juntos, média e desvio juntos, a gente pode olhar de um em um, ou olhar o que nos interessa mais. Basta a gente pegar o Maximum likelihood para os outros valores (nesse caso o desvio) e olhar como fica a distribuição dos valores para o que nos interessa, nesse caso a média.

06

Isso que as pessoas chamam de profile, é olhar so o que interessa, você mantem todos os outros parâmetros no melhor encaixe e olha o que te interessa e tira suas conclusões.

Bem, eu acho legal esse tópico, primeiro que ajuda a fazer uma relação com qual a linha de raciocínio que as análises usam para produzir conclusões, além de ajudar a interpretar mais valores que somente o p quando fazemos alguma análise.

E outra coisa, é que, quando começamos a perceber esse esquema do likelihood, da para perceber melhor a relação entre diferentes análises, e ver como regressão linear, anova, regressão de poisson, regressão binomial ou logística, é tudo o mesmo esquema, variando alguns passos apenas.

Bem até o próximo post para quem teve saco de ler até aqui ^^

#################################################################
#Máxima verossimilhança da média e desvio utilizando
#grid para a busca
#veja o script original em:
#http://www.sortie-nd.org/lme/Course_Schedule_2012/Day_1/Lab%201%20-%20Section%201.R
##################################################################
set.seed(51)
 
#Simulando uma amostra com distribuição normal
media <- 15
desvio <- 2.5
n <- 100
x<-rnorm(n,media,desvio)
#A partir daqui a gente tem que fingir que não sabe os valores
#reais de média e desvio
 
#Histograma
hist(x,main="Histograma",ylab="Frequência",ylim=c(0,50))
 
 
#Grid testando alguns valores de parametros
par(mfrow=c(4,4))
for(i in seq(12,18,by=2)) {
  for(j in seq(2,5,by=1)) {
    hist(x,prob=T,main=paste("Média=",i,"Desvio=",j),ylim=c(0,0.3),xlim=c(6,24))
    curve(dnorm(x,i,j),add=T,col="red",lty=3,lwd=4)
  }
}
 
#Agora vamos fazer um grid para calcular a verosimilhança de diferentes
#combinações de estimativas
 
#Defina as dimensões do grid
range <- 1
increment <- 0.01
 
#Crie um vetor com a mudança dos valores do grid
muhat <- seq(mean(x)-(range*mean(x)),mean(x)+(range*mean(x)),by=mean(x)*increment)
sdhat <- seq(sd(x)-(range*sd(x)),sd(x)+(range*sd(x)),by=sd(x)*increment)
 
#Iniciando uma matrix para salvar os resultados
log_lh <- array(0,dim=c(length(sdhat),length(muhat)))
 
#Cuidado, na função original acho que o autor invertou o loop do j e do i
#Como a matrix original era quadrada não tinha problema, mas se vc alterar o
#grid vc tera um erro ao tentar olhar fora do limite da matrix.
 
for (i in 1:length(sdhat)) {
  for (j in 1:length(muhat)) {
    log_lh[i,j] <-  sum(dnorm(x,muhat[i],sdhat[j],log=T))
  }
}
 
#Essa linha:
#log_lh[i,j] <-  sum(dnorm(x,muhat[i],sdhat[j],log=T))
#É o mesmo que pensar nesse loop
#for (k in 1:length(x)) {
#  log_lh[i,j] <- log_lh[i,j] + dnorm(x[k],muhat[i],sdhat[j],log=T)
#}
#So para ir mais rapido :)
#Mas vale a pena olhar como fica o loop, pelo menos eu
#entendo melhor esses loops com for que operação com vetor
 
#Qual o Maximum likelihood
max_lh<-max(log_lh,na.rm =T)
max_lh
 
# Agora o gráfico do alvo de arco e flecha
xlimit <-range(muhat)
ylimit <-range(sdhat)
clevels <- seq(max_lh - 10,max_lh,by=0.5)
cores<-colorRampPalette(colors=c("blue","red"))
 
#Grafico Alvo
par(mfrow=c(1,1))
filled.contour(x=muhat,y=sdhat,z=log_lh,xlim=c(13,17),ylim=c(1,4),
               levels=clevels,color.palette=cores,frame=F,
               plot.axes = { axis(1); axis(2);
                             contour(x = muhat, y = sdhat, z = log_lh,
                                     xlim = c(13,17), ylim =c(1,4),
                                     levels=seq(min(clevels),max(clevels),
                                       length=11)[-11],drawlabels=F,add=T,lwd=2,
                                     lty=3)
                             points(media,desvio,pch=19,cex=1.5,col="darkgray")
                           }
               )
 
 
#Mas ainda não esta pronto, agora extraimos os indices dos maximos valores da
#maxima verosemelhança
mx<-max.col(log_lh)
 
#Agora vamos linha por linha e achamos qual linha tem o maior valor
#iniciando um valor para ser substituido, so colocar um valor gigante
max <- -1000000000
for (i in 1:nrow(log_lh)) {
  if (log_lh[i,mx[i]] > max)
    {  maxrow <- i
       maxcol <- mx[i]
       max <- log_lh[i,mx[i]]
     }
}
 
#Estimativa da maxima verosimilhança da média
mu_mle <- muhat[maxrow]
mu_mle
mean(x)
 
#Figura 4
plot(muhat,rep(1,length(muhat)),col=ifelse(muhat<=mu_mle,"green","blue"),
     yaxt="n",frame=F,xlim=c(12,16),ylab="",pch=19)
abline(v=mean(x),lwd=2,lty=2,col="red")
legend("topleft",bty="n",col=c("green","blue"),pch=19,
       legend=c("Valores de muhat<=mu_mle",
         "Valores de muhat>mu_mle"))
legend("topright",bty="n",lwd=2,lty=2,col="red",legend="Média Real")
 
 
#Estimativa da maxima verosimilhança do desvio padrão
sd_mle <- sdhat[maxcol]
sd_mle
sd(x)
 
#Usando o comando lm para estimar os valores
summary(lm(x~1))
 
#Agora calculamos os residuos (erros) e vemos se ele são
#normalmente distribuidos com média zero
res <- x - muhat[maxcol]
 
#Histograma com curva estimada
hist(res,freq=F,xlab="Residual",main="Residual standard error")
 
#Nos podemos sobrepor uma PDF normal nesse grafico, dados os MLEs
pdfx <- seq(from = -2*desvio, to = 2*desvio, by = 0.5)
pdfy <- dnorm(pdfx,0,sd_mle)
points(pdfx,pdfy,type="p",pch=19,col="Blue",cex=2)
lines(pdfx,pdfy,type="l",col="Blue",lwd=2)
 
 
#Profile likelihoods são graficos da mudança do likelihood de uma parametro
#em especial enquanto o outro parametro é mantido no maximo.
#Por exemplo o Profile da média
plot(muhat,log_lh[,maxcol],xlim=c(0,30),frame=F)
abline(max_lh-2,0)
abline(v=media,lwd=2,col="red",lty=2)
legend(0,-300,legend="Média Real",lty=2,lwd=2,col="red",bty="n")

Problema – Rosalind – Calculating Protein Mass – PRTM

E aqui vamos nos com o problemas do projeto Rosalind. Bem, em um problema anterior, esse aqui mais especificamente, a gente viu como podemos a partir de uma sequencia, ver qual proteína será produzida. Agora uma coisa legal é que se a gente sabe o que vai ser produzido e algumas características dos tijolinhos que produzem, podemos prever varias características do que ta no final da linha de produção desse sistema todo ae.

O problema te fornece uma tabela com as massas para cada aminoácido, essa tabela tem la na wikipedia, assim como alguns pacotes que englobam funções voltadas para bioinformática tem essas e outras informações para previsões úteis.

Então o python para resolver o serviço, a idéia é bem simples, criamos um dicionario no python com o nome do aminoácido e sua respectiva massa. Agora eu entendi porque diabos as pessoas usam um dicionario e não uma lista para fazer isso. Se essa lista fosse muito longa, sempre para procurar um peso, na pior das hipóteses teríamos que rodar toda a lista procurando, em dados como um dicionario (usando {} chaves), o dicionário tem uma função ja feita nele de hash-table, isso basicamente é um índice remissivo, que torna a busca mais rápida, lembre-se que se queremos achar uma palavra, se você olhar o livro todo você ta ferrado, mas se você pega a primeira letra e vai procurar no índice remissivo no final do livro, serão menos páginas para olhar, logo a busca é mais rápida porque temos que olhar menos coisas, particularmente aqui isso não deve ter efeito nenhum, mas se a lista fosse grande teria 🙂

Certo, mas sabendo disso, nos optamos pelo dicionario ao invés de uma lista e dai é so somar letra a letra os pesos e temos as resposta.
Como so temos uma sequência, a leitura do arquivo é simples, eu so não sei se desse jeito o arquivo fica na memória, deve ficar la tomando espaço caso não usemos o close()

Mas então a função final fica assim:

mono_mass_table = {
    'A': 71.03711,
    'C': 103.00919,
    'D': 115.02694,
    'E': 129.04259,
    'F': 147.06841,
    'G': 57.02146,
    'H': 137.05891,
    'I': 113.08406,
    'K': 128.09496,
    'L': 113.08406,
    'M': 131.04049,
    'N': 114.04293,
    'P': 97.05276,
    'Q': 128.05858,
    'R': 156.10111,
    'S': 87.03203,
    'T': 101.04768,
    'V': 99.06841,
    'W': 186.07931,
    'Y': 163.06333,
}
 
def peso(seq_prot):
    n = 0
    for i in seq_prot:
        n += mono_mass_table[i]
    return n
 
dados = "SKADYEK"
#dados = open("/home/augusto/Downloads/rosalind_prtm.txt").read().strip()
 
print peso(dados)

Que vai nos dar o resultado de:

>>> 821.39192

E todos vivemos felizes para sempre.