Visualizando a distribuição normal multivariada com elipses

Uma distribuição interessante é a distribuição normal multivariada, que é um grupo de distribuições normais.

A parte interessante vem do fato que muito das coisas que a gente olha, são multivariadas, vem de várias distribuições, então entender como ela funciona pode ser útil, para captar melhor os modelos que fazemos para descrever o mundo.

Mas vamos ser mais específico, alguém já ficou pensando o que diabos é essa matriz que sai junto com alguns modelos?

Relembrando esse post aqui

Correlation of Fixed Effects: levsl1 levsl2 levsl3 levsl4 levsl2 0.990 levsl3 0.990 0.990 levsl4 0.990 0.990 0.990 levsl5 0.990 0.990 0.990 0.990

Veja que no final do modelo misto, sai essa matriz ai, que a gente, eu pelo menos, ignorava por um bom tempo, por não saber como diabos funciona uma matriz, para um monte de distribuição normal que são os parâmetros calculados ali em cima no modelo, opa, olha algo repetido aqui, monte de distribuição normal, e a distribuição normal multivariada é um monte de distribuição normal junta. Bem acho que agora já não parece que a gente ta so olhando mais uma distribuição maluca a esmo, e cedo ou tarde, isso vai aparecer.

Uma distribuição normal multivariada, é então como a distribuição normal, possui média e variância, com a diferença que ela é multi, então são várias médias e várias variâncias, so que não somente isso, ela introduz o conceito de correlação entre as amostras. Sem isso, a gente poderia trata-las de forma independente, mas ela é uma distribuição multivariada porque não se trata somente das médias e das variâncias, mas da correlação entre as amostras.

Na página da wikipedia sobre a distribuição, que tem aquele quadrinho legal do lado que sumariza várias coisas legais, como quais são os parâmetros da distribuição, ele fala que entre os parâmetros temos as médias, que são a localização no gráfico, e a matriz de correlação.

Vamos pensar num exemplo que é simples de representar graficamente, que é o de duas distribuições normais, um exemplo 2d.

Para gerar dados de uma distribuição, a gente pode usar a função rnorm, que vai perguntar quantas amostras queremos, e alguns parâmetros, que vão ser a média e o desvio-padrão (que ao quadrado é a variância)

Description Density, distribution function, quantile function and random generation for the normal distribution with mean equal to mean and standard deviation equal to sd. Usage dnorm(x, mean = 0, sd = 1, log = FALSE) pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) rnorm(n, mean = 0, sd = 1)

Mas existe, no pacote MASS, uma função equivalente ao rnorm para a distribuição normal multivariada que é o mvrnorm.

Ele vai precisar do número de amostras que você quer e os parâmetros, que aqui vai ser um vetor de médias, que é a localização das amostras, e o vetor de covariância, que é como se fosse a variância, mas tem a correlação entre as duas distribuições.

Vamos fazer aqui o seguinte exemplo, teremos primeiro como os parâmetros

   Parametros = \left\{       \begin{array}{lr}         \mu = {0,0}\\         \Sigma = \left(\begin{array}{cc} 1&0 \\ 0&1 \end{array} \right)       \end{array}     \right.

Veja então que é um vetor de médias e uma matriz de variância

mu<-c(0,0)
mu
Sigma <- matrix(c(1,0,0,1),2,2)
Sigma

E podemos gerar uma amostra

amostras<-mvrnorm(n=1000, mu, Sigma)

E ver como fica

03

Veja que basicamente é um círculo, ou seja a gente tem o centro do círculo mais ou menos no ponto [0,0] e o raio vai ser mais ou menos a variância, mais ou menos porque é uma distribuição normal, então temos variação.

Agora so como exercício, olha que legal. Olhando o circulo no wikipedia, mais especificamente a equação paramétrica, a gente consegue reproduzir ela no R e desenha um círculo.

plot(amostras)
t<-seq(0, 2*pi, length = 100)
r<-1
a<-0
b<-0
points(a+r*cos(t),b+r*sin(t),type="l",lwd=5,lty=3,col="red")

04

Muito legal né, inútil mas legal. So que com esse círculo ali, a gente pode ver quantos porcentos desses monte de ponto estão dentro do círculo. Como?
Por isso tem esse monte de formula de círculo, várias formas de equação, que cada uma serve pra uma coisa, a forma paramétrica a gente usa para desenhar, como ali em cima. A equação do círculo, que tem ali em cima da equação paramétrica no wikipedia, a gente faz uma função para verificar quem está dentro, quem está fora.

> ifelse(amostras[,1]^2 + amostras[,2]^2 <= r^2, TRUE, FALSE) [1] TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE TRUE . . . [991] TRUE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE

Agora a gente pode contar, quantos pontos estão dentro desse círculo.

> sum(ifelse(amostras[,1]^2 + amostras[,2]^2 <= r^2, TRUE, FALSE))/nrow(amostras) [1] 0.387

So contar com sum, lembrando que TRUE vale 1 e FALSE 0, então a gente soma esse vetor e divide pelo total de amostras e temos a nossa contagem.

Mas a gente pode tentar achar o valor de raio que englobe 95% dos pontos, que desenhe um circulo que englobe 95% dos pontos. Como? Repita isso para vários valores de raio.

> for(r in seq(0,4,0.5)) { + cat("Raio ",r," cobre ",(sum(ifelse(amostras[,1]^2 + amostras[,2]^2 <= r^2, TRUE, FALSE))/nrow(amostras))*100,"% dos pontos \n",sep="") + } Raio 0 cobre 0% dos pontos Raio 0.5 cobre 10.8% dos pontos Raio 1 cobre 38.7% dos pontos Raio 1.5 cobre 68.7% dos pontos Raio 2 cobre 86.6% dos pontos Raio 2.5 cobre 94.8% dos pontos Raio 3 cobre 98.7% dos pontos Raio 3.5 cobre 99.8% dos pontos Raio 4 cobre 100% dos pontos

Ha, olha ai, um raio de 2.5 cobriria praticamente 95% dos pontos. E como usar essa informação de um jeito legal, a gente pode fazer tipo um kernel pros pontos, ir colocando uns círculos pra gente ter uma noção melhor da distribuição dos pontos nesse plano.

05

Massa né, mas agora esse esquema usando o círculo vai pifar se tivermos covarianças. Vamos alterar a matriz de covariância.

Vamos usar os seguintes parâmetros agora para a distribuição

   Parametros = \left\{       \begin{array}{lr}         \mu = {0,0}\\         \Sigma = \left(\begin{array}{cc} 1 & 1 \\ 1 & 1 \end{array} \right)       \end{array}     \right.

E olha o que vira.

06

Agora ja da pra compreender melhor essa matriz de covariância.

Veja que as variâncias tanto n eixo x como no eixo y são próximos de 1. Que é o que ta na diagonal central da matriz.
E o 0.8 ali é o quanto um eixo esta correlacionado um com o outro.

> var(amostras[,1]) [1] 0.9343743 > var(amostras[,2]) [1] 0.9513073 > cor(amostras) [,1] [,2] [1,] 1.0000000 0.7826158 [2,] 0.7826158 1.0000000

E a gente pode ir mexendo nesses parametros para ver o como vai ficando a figura.

07

E o circulo acabou de ir pro bebeleu, que ele não serve mais para tentar representar isso ai, mas a gente tem a elipse, que pode servir.

Olha no wikipedia la, a formula da elipse, é diferente, mas podemos simplesmente transcrever ela pro R a plotar para ver.

t<-seq(0, 2*pi, length = 100)
Xc<-0
Yc<-0
alpha<-pi/3
a<-1
b<-1
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")

E agora ao plotar no R vemos a elipse.

08

So que a elipse tem bem mais parâmetros que o circulo, que só tem o raio, e o ponto de origem. Aqui a gente tem o ponto de origem, um a e b e angulo alpha que vão controlar a forma da elipse. E o negócio é sempre sair jogando valor e ver no que da para pegar o feeling, mudando somente o alpha, veja o que da para fazer

08

Agora não vai dar mais para sair usando força bruta, como fizemos com o círculo, num vai dar muito certo aqui, porque são muitas combinações possíveis de parâmetros. Mas no pacote MASS do R, existe uma função chamada cov.mve que é um estimador pronto dos parâmetros da distribuição multivariada, que é o que precisamos para desenhar a elipse. Então podemos usar ela, e assim como no caso do circulo, podemos traçar uma elipse que cobre uma % desejada dos pontos.

09

Agora é uma questão de cobertura, assim como fizemos com o círculo, podemos desejar outros tamanhos de coberturas, ou %.

10

Mas a elipse so vai servir para a distribuição multivariada normal, as vezes quando a gente ta olhando esse tipo de dados, em 2d lembra muito dados espaciais, como quando olhando mapas, como no caso do John Snow, aqui a gente precisa de outras abordagens, como a distribuição posterior, que tem uma função pronta para o calculo no pacote de nosso querido Ben Bolker (O que ele faria?).

12

Outras soluções, talvez para uma primeira visualização, é usar imagens, onde cada pixel tem a densidade local expressa na intensidade de cor, por exemplo usando o hex do ggplot2.

13

E é isso ai, a gente so deu uma olhada numa distribuição, pra matar a curiosidade, e viu algumas formas de representação dela, bem como outras formas de representação de dados em 2d, para sei la, facilitar a visualização. Qualquer mapinha você vai precisar dessas coisas. A, outra coisa interessante.

01

02

Veja como a média, que é a localização, pode ser mantida, mesmo mudando os parâmetros de covariância. Veja que em cima e do lado do gráfico, a gente vê a distribuição de uma dimensão, mas no meio a gente ve a distribuição em duas dimensões.

Fim de ano ta difícil ter tempo para mexer com esses posts, mas já já acaba. Por enquanto é isso ai, apenas uma olhada em no mundo 2d. Os scripts sempre no repositório recologia e até o próximo post.

#Muito do codigo veio desse post
#http://www.sumsar.net/blog/2014/11/how-to-summarize-a-2d-posterior-using-a-highest-density-ellipse/
 
library(MASS)
library(cluster)
library(plotrix)
library(emdbook)
library(ggplot2)
library(hexbin)
#install.packages("devtools")
#library(devtools)
#install_github("rasmusab/bayesian_first_aid")
library(BayesianFirstAid)
 
Sigma <- matrix(c(1,0,0,1),2,2)
Sigma
 
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
 
#
plot(amostras)
t<-seq(0, 2*pi, length = 100)
r<-1
a<-0
b<-0
points(a+r*cos(t),b+r*sin(t),type="l",lwd=5,lty=3,col="red")
 
#Testando quem está dentro do circulo de raio r
sum(ifelse(amostras[,1]^2 + amostras[,2]^2 <= r^2, TRUE, FALSE))/nrow(amostras)
 
for(r in seq(0,4,0.2)) {
    cat("Raio ",r," cobre ",(sum(ifelse(amostras[,1]^2 + amostras[,2]^2 <= r^2, TRUE, FALSE))/nrow(amostras))*100,"% dos pontos \n",sep="")
}
 
#
plot(amostras,frame=F)
i=2
for(r in c(3.1,2.5,1.2)) {
    points(a+r*cos(t),b+r*sin(t),type="l",lwd=5,lty=3,col=i)
    i=i+1
}
legend("topright",title="cobertura",legend=c("99%","95%","50%"),bty="n",lwd=4,lty=3,col=c(2,3,4))
 
#Distribuição não redondinha
Sigma <- matrix(c(1,0.8,0.8,1),2,2)
Sigma
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
#
plot(amostras)
 
var(amostras[,1])
var(amostras[,2])
cor(amostras)
 
 
#
par(mfrow=c(2,2))
Sigma <- matrix(c(1,0.8,0.8,1),2,2)
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
plot(amostras)
addtable2plot (x=-2,y=2,table=Sigma,display.colnames=F,bty="o",hlines=T,vlines=T)
 
Sigma <- matrix(c(1,-0.8,-0.8,1),2,2)
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
plot(amostras)
addtable2plot (x=2,y=2,table=Sigma,display.colnames=F,bty="o",hlines=T,vlines=T)
 
Sigma <- matrix(c(1,0.8,0.8,2),2,2)
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
plot(amostras)
addtable2plot (x=-3,y=2,table=Sigma,display.colnames=F,bty="o",hlines=T,vlines=T)
 
Sigma <- matrix(c(2,0.8,0.8,2),2,2)
amostras<-mvrnorm(n=1000, rep(0, 2), Sigma)
plot(amostras)
addtable2plot (x=3,y=-3,table=Sigma,display.colnames=F,bty="o",hlines=T,vlines=T)
 
 
 
#Desenhando Elipse
t<-seq(0, 2*pi, length = 100)
Xc<-0
Yc<-0
alpha<-pi/3
a<-1
b<-1
 
#
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")
 
 
#
par(mfrow=c(2,2))
alpha<-pi/3.5
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")
legend("topleft",legend=paste("alpha=",round(alpha,2)),bty="n")
alpha<-pi/3
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")
legend("topleft",legend=paste("alpha=",round(alpha,2)),bty="n")
alpha<-pi/2
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")
legend("center",legend=paste("alpha=",round(alpha,2)),bty="n")
alpha<-pi/1.5
plot(Xc+a*cos(t)*cos(alpha)-b*sin(t)*sin(alpha),
     Yc+a*cos(t)*sin(alpha)-b*sin(t)*cos(alpha)
     ,type="l",lwd=3,lty=1,col="black")
legend("topright",legend=paste("alpha=",round(alpha,2)),bty="n")
 
 
#Ajuste de elipse
ajuste <- cov.mve(amostras, quantile.used = nrow(amostras) * 0.50)
pontos_na_ellipse <- amostras[ajuste$best, ]
limite_ellipse <- predict(ellipsoidhull(pontos_na_ellipse))
plot(amostras)
points(limite_ellipse[,1],limite_ellipse[,2],type="l",lwd=5,lty=3,col="red")
legend("topleft", "50%",lwd=5,lty=3,col="red")
 
 
 
#Varios ajustes
plot(amostras)
i<-2
for(cobertura in c(0.99, 0.95, 0.5)) {
    ajuste <- cov.mve(amostras, quantile.used = nrow(amostras) * cobertura)
    pontos_na_ellipse <- amostras[ajuste$best, ]
    limite_ellipse <- predict(ellipsoidhull(pontos_na_ellipse))
    points(limite_ellipse[,1],limite_ellipse[,2],type="l",lwd=5,lty=3,col=i)
    i<-i+1
}
legend("topleft",legend=c("99%","95%","50%"),lwd=5,lty=3,col=2:4)
 
 
#Plot com poligonos
plot(amostras,type="n")
i<-2
for(cobertura in c(0.99, 0.95, 0.5)) {
    ajuste <- cov.mve(amostras, quantile.used = nrow(amostras) * cobertura)
    pontos_na_ellipse <- amostras[ajuste$best, ]
    limite_ellipse <- predict(ellipsoidhull(pontos_na_ellipse))
    polygon(limite_ellipse, col = i, border = NA)
    i<-i+1
}
points(amostras[,1],amostras[,2],pch=19)
legend("topleft",legend=c("99%","95%","50%"),pch=1,col=2:4)
 
#Distribuição posterior
plot(amostras)
HPDregionplot(amostras, prob = c(0.95, 0.75, 0.5), col=2:4, lwd=3,lty=3, add=TRUE)
legend("topleft", legend = c("95%", "75%", "50%"), col = 2:4, lty=3, lwd=3)
 
 
 
#ggplot2 hexagonos
qplot(amostras[,1], amostras[,2], geom=c("hex"))
 
#Bayes plot
amostras<-mvrnorm(n=50, c(0,0), Sigma)
ajuste <- bayes.cor.test(amostras[,1], amostras[,2])
plot(ajuste)

Simulando sequências de DNA com Cadeias de Markov (Markov Chain)

03

Opa, uma forma excelente de entender melhor analise estatística, é simular dados. Eu sempre digo isso por aqui, mas aprendi nos livros do Marc Kéry.
Sempre que estamos simulando dados, estamos fazendo o processo inverso da análise de dados. Na análise a gente olha amostras e tenta entender como elas são geradas, quais as regras, para biologia as regras impostas pela natureza, que produzem essas amostras.

Agora para aprender melhor a construir árvores filogenéticas, usar elas para corrigir filogenia em análises como aqui, a gente precisa da capacidade de simular sequâncias de DNA, para testar nossas análises, se elas funcionam.Mas como sempre, nem tudo é tão simples como parece.

O DNA é composto de um alfabeto de quatro letras, A de adenina, C de citosina, T de timina e G de guanina.

O jeito mais simples de fazer isso, é simular o sorteio repetido dentro desse alfabeto.
Vamos supor um alfabeto mais simples, uma moeda, a moeda só cara ou coroa, então uma moeda honesta a gente simularia assim no R

> sample(c("Cara","Coroa"),30,rep=TRUE,prob=c(0.5,0.5)) [1] "Coroa" "Coroa" "Coroa" "Cara" "Coroa" "Cara" "Cara" "Cara" "Coroa" "Coroa" "Cara" "Cara" [13] "Cara" "Coroa" "Cara" "Coroa" "Cara" "Coroa" "Cara" "Cara" "Coroa" "Cara" "Coroa" "Coroa" [25] "Cara" "Coroa" "Coroa" "Coroa" "Cara" "Cara"

Usamos a função sample, que o primeiro argumento é um vetor com o alfabeto que queremos retirar as amostras, depois informamos que o sorteio é com repetição, o default é sem repetição. Temos que informar também quantas amostras queremos, aqui 30, e por ultimo um vetor de probabilidades, que para uma moeda honesta é 50% para cara e 50% para coroa.

Agora a gente pode seguir a mesma linha e simular dessa forma uma sequencia de DNA. Basta trocar o alfabeto usado no sorteio e o vetor de probabilidades, que tem que ser do tamanho do alfabeto, claro.

> sample(c("A","G","C","T"),30,rep=TRUE,prob=c(0.25,0.25,0.25,0.25)) [1] "G" "C" "G" "C" "G" "A" "T" "T" "A" "G" "G" "A" "T" "G" "C" "T" "C" "C" "A" "A" "C" "T" "A" "A" [25] "T" "G" "T" "A" "G" "C"

Mas talvez sequências geradas assim sejam simplesmente um monte de letras, e sequências de DNA não são tão aleatórias assim. Como assim?
Bem por exemplo, podemos citar o conteúdo de GC na sequência. A se liga com T e C com G, so que quando A se liga com T, essa ligação possui duas pontes de hidrogênio, enquanto C com G formam três pontes de hidrogênio, assim quanto mais o conteúdo de CG na cadeia, mais estável ela fica, do jeito que estamos gerando sequências, a tendência é elas sempre terem quantidades iguais de AT e CG, já que estamos sorteando cada base de forma independente.
Poderíamos gerar sequências com maior conteúdo de CG simplesmente aumentando a probabilidade dessas letras do nosso alfabeto de DNA, mas podemos também condicionar a simulação a quando tiver acumulando muito AT, termos uma chance maior de tender a um estado com mais CG usando cadeias de markov para simular o DNA.

Vamos tentar fazer um exemplo pra ficar mais fácil.

A gente tem uma sequência com variáveis \{X_1, X_2, X_3, X_4 , \dots, X_n \} e cada variável tem um estado E, estado pra gente é qual letra do do alfabeto um X_n qualquer está, sendo o alfabeto ATCG, bem E = \{A,\ C,\ ,T,\ G\}X_n.

Mas para pegar o jeito primeiro, vamos assumir um alfabeto mais simples, pense que E = \{1,\ 2\} e temos a sequência \{X_1, X_2, X_3, X_4 , \dots, X_n \}.
Quando olhamos para X_n = i qualquer, o nosso processo está no estado i no tempo n.

A expressão P(X_1=i) denota a probabilidade que o processo esteja no estado i no tempo 1.

O evento no qual o processo muda seu estado de i para j (transição) entre os tempos um para o dois corresponde ao evento (X_2 = j \vert X_1 = i), sendo que essa barrinha significa “dado que”, lembre-se de probabilidade condicional que a gente comentou aqui.
Essa probabilidade é denotada como P(X_2 = j \vert X_1 = i), mas de modo geral, a probabilidade de transição de i para j do ponto n ao n+1 pode ser descrita de forma mais geral como P(X_{n+1} = j \vert X_n = i). Dessa forma, podemos pensar ja nas probabilidades relacionadas a toda a sequencia que vamos gerar e probabilidades ligadas a todo o alfabeto, todos os estados.
Todas essas probabilidades pode ser visualizadas de uma forma mais simples numa matriz, a matriz de probabilidades de transição, onde cada elemento da matriz é P_{ij}(X_{n+1} = j \vert X_n = i).

A matrix de probabilidade de transição contem a probabilidade discreta (condicional) em cada uma de suas linhas (que somam 1 sempre). Para um processo de Markov, o estado no ponto n+1 depende do estado do ponto n, mas não nos estados dos tempos anteriores.

Então no nosso processo com o alfabeto E = \{1,\ 2\}, temos o seguinte.

  \begin{array}{ccc}   & 1 & 2 \\  1 & P_{11} & P_{12} \\  2 & P_{21} & P_{22} \end{array}

E como a gente viu la em cima, isso implica no seguinte.

  \begin{array}{ccc}   & 1 & 2 \\  1 & (X_{n+1} = 1 \vert X_n = 1) & (X_{n+1} = 2 \vert X_n = 1) \\  2 & (X_{n+1} = 1 \vert X_n = 2) & (X_{n+1} = 2 \vert X_n = 2) \end{array}

Certo, agora a gente pode colocar isso na pratica, a gente pode fazer uma função para simular sequencias de dna usando cadeias de Markov

markov1 <- function(x,P,n) {
    seq <- x
    for(k in 1:(n-1)){
        seq[k+1] <- sample(x, 1, replace=TRUE, P[seq[k],])
    }
    return(seq)
}

Lembrando que nosso exemplo ainda não é o dna, os estados são E = \{1,\ 2\}.
Na função a gente recebe x, que é o alfabeto, dai podemos usar a função sample para sortear o próximo estado, e usar as probabilidades da linha da matriz que do estado atual.
Bem temos que iniciar também a cadeia, aqui vamos simplesmente começar com x.

Temos que estabelecer uma matriz de transição, vamos inventar alguns valores aqui.

  \begin{array}{ccc}   & 1 & 2 \\  1 & {5\over{6}} & {1\over{6}} \\  2 & {1\over{2}} & {1\over{2}} \end{array}

No R é so a gente fazer uma matriz.

> P <- matrix(c(1/6,5/6,0.5,0.5), 2, 2, byrow=TRUE) > P [,1] [,2] [1,] 0.1666667 0.8333333 [2,] 0.5000000 0.5000000

Ai temos que arrumar os nomes de linhas e colunas da matriz, para ela funcionar corretamente no samples, estabelecemos os estado atual e é só simular nossa sequencia.

> rownames(P)<-c(1,2) > colnames(P)<-c(1,2) > x<-c(1,2) > markov1(x,P,30) [1] 1 2 2 2 1 2 1 2 2 1 2 2 2 2 2 2 2 1 2 1 1 2 2 2 2 1 2 2 1 2

Certo, agora como isso funciona para sequências de DNA?
Vamos tentar adicionar a propriedade que falamos la em cima, vamos colocar uma maior quantidade de C e G na sequência.

P <- matrix(c(
1/6,5/6,0,0,
1/8,2/4,1/4,1/8,
0,2/6,3/6,1/6,
0,1/6,3/6,2/6),4,4,byrow=TRUE)
rownames(P) <- colnames(P) <- StateSpace <- c("a","c","g","t")

O que a gente precisa é pensar na matriz de transição.
Olha la, depois de um A, a gente poe uma grande chance de ir para C
Se a gente tiver em C, a gente tem a maior chance de ficar em C, segundo maior em G.
E assim vai.

> P a c g t a 0.1666667 0.8333333 0.00 0.0000000 c 0.1250000 0.5000000 0.25 0.1250000 g 0.0000000 0.3333333 0.50 0.1666667 t 0.0000000 0.1666667 0.50 0.3333333

A gente ainda tem que fazer uma pequena modificação, vamos mandar a parte uma chance de cada uma das bases para iniciar a cadeia, chamando ela de pi0. No mais temos o mesmo sistema anterior, e vamos guardar os resultado num vetor de character.

markov2 <- function(StateSpace,P,pi0,n) {
    seq <- character(n)
    seq[1] <- sample(StateSpace, 1, replace=TRUE, pi0)
    for(k in 1:(n-1)) {
        seq[k+1] <- sample(StateSpace, 1, replace=TRUE, P[seq[k],])
    }
    return(seq)
}

Então geramos um vetor para iniciar a sequência.

pi0 <- c(1/4,1/4,1/4,1/4)

E bang, usamos nossa função, podemos ver o começo da cadeia com a função head, e podemos usar a função table e ver quanto temos de cada base.

> x <- markov2(StateSpace,P,pi0,1000) > head(x) [1] "c" "t" "t" "t" "g" "g" > table(x) x a c g t 54 399 364 183

E olha ae, temos bastante C e G, pouquinho A a um pouco mais de T.

Mudando um pouco a matriz de transição, podemos ainda aumentar a chance de alguns aminoácidos, lembre-se que três bases vão ser codificados em um aminoácido, As sequencia que produzem phenylalanine são TTT e TTC. No começo da tabelinha de codon a gente ve esses dois grupinhos.

Geramos uma nova matriz de transição

P <- matrix(c(.01,.01,.01,.97,
.01,.01,.01,.97,
.01,.01,.01,.97,
.01,.28,.01,0.70),4,4,byrow=T)
rownames(P) <- colnames(P) <- StateSpace <- c("a","c","g","t")
x <- markov2(StateSpace,P,pi0,30000)

E geramos nossa sequência de DNA.

> library(seqinr) > table(getTrans(x)) * A C D F H I L M N P Q R S T V W Y 2 2 62 1 5157 20 105 2255 2 1 15 1 23 2192 1 93 1 67

Podemos usar a função getTrans para traduzir a sequencia em aminoácidos, mas a gente ja fez essa função aqui.
E pimba, temos 5157 phenylalanina, enquanto o aminoácido que tem a segunda maior quantidade é a Leucina que tem 2255.

Agora você pode pensar, poxa, agora como diabos eu vou ficar inventando essas matrizes de transição? Eu não sou criativo, como eu vou tirar esse monte de número da minha cabeça?
Mas você não precisa ficar inventado, a na verdade essa matriz de transição já existe, nas sequencias que a gente já sequenciou. Você pode pegar uma sequência qualquer no genbank e estimar a matriz de transição da qual ela veio.

A gente pode contar o número de transições e dividir pelo total da linha. Podemos assim estimar as probabilidades de transição de maior probabilidade, mais verossimilhantes. Alias qualquer probabilidade zero será exato nesse caso, porque serão casos que nunca existem.

Primeiro a gente conta as transições.

> nr <- count(x,2) > nr aa ac ag at ca cc cg ct ga gc gg gt ta tc tg tt 4 1 2 295 66 56 61 6419 1 4 5 267 231 6542 208 15837

Temos que fazer uma matriz para receber os resultados e colocar as contagens nele

A <- matrix(NA,4,4)
A[1,1]<-nr["aa"]; A[1,2]<-nr["ag"]; A[1,3]<-nr["ac"]; A[1,4]<-nr["at"]
A[2,1]<-nr["ga"]; A[2,2]<-nr["gg"]; A[2,3]<-nr["gc"]; A[2,4]<-nr["gt"]
A[3,1]<-nr["ca"]; A[3,2]<-nr["cg"]; A[3,3]<-nr["cc"]; A[3,4]<-nr["ct"]
A[4,1]<-nr["ta"]; A[4,2]<-nr["tg"]; A[4,3]<-nr["tc"]; A[4,4]<-nr["tt"]

E vemos o seguinte

> A a c g t a 4 2 1 295 c 1 5 4 267 g 66 61 56 6419 t 231 208 6542 15837

Igual a matriz de probabilidade de transições, mas temos contagens aqui.
ai a gente transforma essas contagens em probabilidade, dividindo pela soma da linha.

rowsumA <- apply(A, 1, sum)
Phat <- sweep(A, 1, rowsumA, FUN="/")
rownames(Phat) <- colnames(Phat) <- c("a","g","c","t")

E podemos agora comparar os valores originais

> P a c g t a 0.01 0.01 0.01 0.97 c 0.01 0.01 0.01 0.97 g 0.01 0.01 0.01 0.97 t 0.01 0.28 0.01 0.70

Contra os valores estimados

> round(Phat,3) a g c t a 0.013 0.007 0.003 0.977 g 0.004 0.018 0.014 0.964 c 0.010 0.009 0.008 0.972 t 0.010 0.009 0.287 0.694

Fica bem próximo não, claro que esperamos melhores resultado com sequências maiores, mas é um resultado expressivo. Pensando desse jeito, a gente poderia pegar qualquer sequência conhecida, estimar a matriz de transição dela e simular outras sequências para testar analises. Como a gente disse, a gente faz analise de dados para estimar um modelo, mas com o modelo a gente simula dados.

Ficamos por aqui, por agora, porque ainda vamos melhorar isso ae. Os scripts sempre no repositório recologia e é isso, até o próximo post.

Referências:
Wim P. Krijnen – Applied Statistics for Bioinformatics using R

sample(c("Cara","Coroa"),30,rep=TRUE,prob=c(0.5,0.5))
 
sample(c("A","G","C","T"),30,rep=TRUE,prob=c(0.25,0.25,0.25,0.25))
 
markov1 <- function(x,P,n) {
    seq <- x
    for(k in 1:(n-1)){
        seq[k+1] <- sample(x, 1, replace=TRUE, P[seq[k],])
    }
    return(seq)
}
 
P <- matrix(c(1/6,5/6,0.5,0.5), 2, 2, byrow=TRUE)
 
rownames(P)<-c(1,2)
colnames(P)<-c(1,2)
x<-c(1,2)
 
markov1(x,P,30)
 
markov2 <- function(StateSpace,P,pi0,n) {
    seq <- character(n)
    seq[1] <- sample(StateSpace, 1, replace=TRUE, pi0)
    for(k in 1:(n-1)) {
        seq[k+1] <- sample(StateSpace, 1, replace=TRUE, P[seq[k],])
    }
    return(seq)
}
 
P <- matrix(c(
1/6,5/6,0,0,
1/8,2/4,1/4,1/8,
0,2/6,3/6,1/6,
0,1/6,3/6,2/6),4,4,byrow=TRUE)
 
rownames(P) <- colnames(P) <- StateSpace <- c("a","c","g","t")
 
pi0 <- c(1/4,1/4,1/4,1/4)
 
x <- markov2(StateSpace,P,pi0,1000)
 
head(x)
table(x)
 
P <- matrix(c(.01,.01,.01,.97,
.01,.01,.01,.97,
.01,.01,.01,.97,
.01,.28,.01,0.70),4,4,byrow=T)
rownames(P) <- colnames(P) <- StateSpace <- c("a","c","g","t")
x <- markov2(StateSpace,P,pi0,30000)
 
#install.packages("seqinr")
library(seqinr)
table(getTrans(x))
 
nr <- count(x,2)
nr
 
A <- matrix(NA,4,4)
rownames(A) <- colnames(A) <- c("a","c","g","t")
A[1,1]<-nr["aa"]; A[1,2]<-nr["ag"]; A[1,3]<-nr["ac"]; A[1,4]<-nr["at"]
A[2,1]<-nr["ga"]; A[2,2]<-nr["gg"]; A[2,3]<-nr["gc"]; A[2,4]<-nr["gt"]
A[3,1]<-nr["ca"]; A[3,2]<-nr["cg"]; A[3,3]<-nr["cc"]; A[3,4]<-nr["ct"]
A[4,1]<-nr["ta"]; A[4,2]<-nr["tg"]; A[4,3]<-nr["tc"]; A[4,4]<-nr["tt"]
 
A
 
rowsumA <- apply(A, 1, sum)
Phat <- sweep(A, 1, rowsumA, FUN="/")
rownames(Phat) <- colnames(Phat) <- c("a","g","c","t")
 
P
round(Phat,3)

Encontrando a máxima verossimilhança analiticamente

Em 2012, a gente comentou o que é e como encontrar a máxima verossimilhança de uma distribuição binomial aqui, agora vamos olhar com mais carinho o que fizemos la e tentar entender um pouco melhor como esse negócio funciona de vedadinha.

Para uma única observação da distribuição binomial, a verossimilhança (ou likelihood) de k ocorrências de um total de N tem uma probabilidade p (e.g.k arbustos com aranhas de um total de N arbustos ou k passarinhos parasitados de um total N de passarinhos observados) é:

\mbox{Prob}(k|p,N)=\binom{N}{k} p^k (1-p)^{N-k}

Pra quem não sabe o que \binom{N}{k} é, de uma olhada em combinação.

Mas continuando, agora se a gente tem n observações, cada uma com o mesmo número N de arbustos (no caso das aranhas, ou pode ser n observações com N passarinhos em cada amostra, ou o que diabos for do seu interesse), e o número de arbustos habitados por aranhas da i^{\mbox{th}} observação é k_i, então a verosimilhança (ou likelihood) vai ser:

{\cal L} = {\prod_{i=1}^n} {\binom{N}{k_i} p^{k_i} (1-p)^{N-k_i}}
Quem nunca viu esse simbolo {\prod_{i=1}^n}, ele é a mesma coisa que somatório, mas ao invés de somar a gente multiplica as coisas, ou seja ele é o simbolo do produtório, lembre-se que a chance de 2 coisas acontecerem é uma coisa vezes a outra, lembre-se dos “e” ou “ou” da estatística.

Nos já comentamos o porque as pessoas transformam a verosimilhança em log (log-likelihood), que segundos as regras de log, uma multiplicação vira uma soma, além de ter que lidar com números com mais casas decimais, permitindo maior precisão, mas no caso do log-likelihood, não sei como escreve em português, se log-verossimilhança ou verossimilhança logarítmica, se alguém puder me corrigir seria legal, mas o que temos é:

 L=\sum_{i=1}^n (\log \binom{N}{k_i} + k_i \log p + (N-k_i) \log (1-p))

No R a gente escreve isso ai em cima assim:

sum(dbinom(k,size=N,prob=p,log=TRUE))

Agora a gente vê formulas pra caramba em ciência, faz muitas contas que simplesmente sabemos que é. E dificilmente descobrimos o porque é, pegue por exemplo quando a gente calcula a prevalência, o total de casos ocorrido em todas as observações pelo total de casos avaliado em todo o experimento de algum parasita.

Agora vamos olhar algo legal relativo a equação ali de cima, partindo de uma abordagem analítica (leia calculo ^^)

Nesse caso nos podemos resolver o problema analiticamente, diferenciando e equação a partir do p e igualando a derivada a zero. Existe também esse pequeno post aqui sobre derivadas no R. Então a gente não precisa necessariamente saber fazer essa conta, mas é legal o que se desenrola daqui.

Primeiro lembre-se como se parece o gráfico da verossimilhança de uma distribuição binomial.

.

Veja que ele pode começar em baixo, vai chegar uma hora que sobe, numa barriguinha e depois vai descer.

Se você sair na calculando derivadas e desenhando retas tangente assim.

01

O topo da barriguinha do gráfico da verossimilhança da distribuição binomial vai ter uma reta paralela ao eixo x, ou seja a inclinação vai ser zero, ou seja com derivada zero.

Então suponha que \hat p é a máxima verossimilhança da estimativa do parâmetro p que satisfaça:

 \frac{dL}{dp} =  \frac{d\sum_{i=1}^n \left(  \log \binom{N}{k_i} + k_i \log p + (N-k_i) \log (1-p) \right)}{dp}  =  0.

Como a derivativa das soma é igual a soma das derivativas, temos:

  \sum_{i=1}^n \frac{d \log \binom{N}{k_i}}{dp} +  \sum_{i=1}^n \frac{d k_i \log p}{dp} +  \sum_{i=1}^n \frac{d (N-k_i) \log (1-p)}{dp} = 0

O termo \log \binom{N}{k_i} é uma constante em relação ao p, logo sua derivativa é zero, e assim o primeiro termo desaparece.

Como k_i e (N-k_i) são fatores constantes, eles saem das derivativas e a equação então se torna:

  \sum_{i=1}^n k_i \frac{d \log p}{dp} +  \sum_{i=1}^n (N-k_i) \frac{d  \log (1-p)}{dp} = 0.

A derivativa de \log p é 1/p, então a regra da cadeia diz que a derivada de

\log (1-p) é

d(\log (1-p))/d(1-p) \cdot d(1-p)/dp = -1/(1-p).

Vamos escrever o valor particular de p que queremos como \hat p.

Então…

 \frac{1}{\hat p} \sum_{i=1}^n k_i - \frac{1}{1-\hat p} \sum_{i=1}^n (N-k_i)  =  0

Jogamos aquele pedacinho negativo pro outro lado, ou como é o correto que aprendi no Khan Academy, somamos o que está depois do igual dos dois lados.
 \frac{1}{\hat p} \sum_{i=1}^n k_i  =  \frac{1}{1-\hat p} \sum_{i=1}^n (N-k_i)

Dai a gente troca os 1-\hat p e \hat p que estão embaixo do 1 de lado…
 (1-\hat p) \sum_{i=1}^n k_i  =  \hat p \sum_{i=1}^n (N-k_i)

Depois aplica a distributiva ali no 1-\hat p,
 \sum_{i=1}^n k_i - \hat p \sum_{i=1}^n k_i =  \hat p \sum_{i=1}^n (N-k_i)

Passa pro outro lado
 \sum_{i=1}^n k_i =  \hat p \sum_{i=1}^n (N-k_i) + \hat p \sum_{i=1}^n k_i

E coloca o \hat p em evidência…
 \sum_{i=1}^n k_i  =  \hat p \left( \sum_{i=1}^n k_i + \sum_{i=1}^n (N-k_i) \right)

Agora note que aqui,\sum_{i=1}^n k_i + \sum_{i=1}^n (N-k_i) , a gente ta somando os casos que aconteceram, k, com o total de casos N, menos os casos que acontecerem k, poxa, isso simplismente N.

 \sum_{i=1}^n k_i = \hat p \sum_{i=1}^n N

Sem desistir de entender o que está acontencedo, a somatoria \hat p \sum_{i=1}^n N é o numero de observações com N amostras, ou  \hat p n N , depois vamos ainda passar tudo que multiplica o \hat p pro outro lado dividindo….
 \sum_{i=1}^n k_i  =  \hat p n N

E acabamos com…
 \hat p  =  \frac{\sum_{i=1}^n k_i}{nN}

Então a estimativa da máxima verissimilhança \hat p, é apenas o fração do total ocorrências de algo, \sum k_i, dividido por pela soma de todas as observações, nN.

Por exemplo o total de arbustos com aranhas pelo total de arbustos vistoriados, ou o total de passarinhos parasitados dividido pelo total de passarinhos verificados.

Parece um bocado de esforço para provar o obvio, que a melhor estimativa do parasitismo é meramente a frequência de casos de parasitismos, mas a gente começa a ver a relação daquela continha que sempre fazemos, com as curvas que likelihood que gostamos de desenhar.

Outras distribuições, como a distribuição do poisson se comportam de forma similiar. Se a gente diferenciar a verissimilhança (likelihood) ou a verossimilhança logaritimizada (log-likelihood) e resolver para a estimativa da máxima verossimilhança, nos conseguimos a resposta.
Com um pouco mais de esforço (ja que a distribuição normal tem dois parâmetros, média e desvio padrão), podemos fazer o mesmo para média e desvio padrão da distribuição normal (também conhecida com Gaussiana).

Para a distribuição de poisson, a estimativa do parâmetro de frequência \hat \lambda é igual a média das contagens observadas por amostra.

Para a distribuição normal, com dois parâmetros \mu e \sigma^2, nos temos que calcular as derivativas parciais da verossimilhança para os dois parâmetros, e resolver as duas equações. (\partial L/\partial \mu =\partial L/\partial \sigma^2 = 0).

E como a gente disse, sem conta nenhuma, a resposta obvia daquilo que a gente calculou a vida inteira.
\hat \mu=\bar x (a estimativa é a média observada nas amostras) e a variância é \hat{\sigma^2}=\sum (x_i-\bar x)^2/n (a estimativa da variância é a variância das amostras)

Um pequeno detalhe, é que a máxima verossimilhança é uma estimativa tendenciosa da variância, dividindo a soma dos quadrados \sum (x_i-\bar x)^2 por n-1 ao invés de n é o usado comumente. Mas a explicação fica para a próxima hehe.

Para outras distribuições simples como a binomial negativa, e mais alguns problemas complexos, não há uma solução analítica simples assim, como a que vimos aqui. E temos que aplicar soluções numéricas para encontrar a máxima verossimilhança para parâmetros das distribuições e modelos que temos interesse. Para praticamente tudo, existem programas prontos, mas o ponto legal do calculo visto aqui é criar alguma intuição do que está acontecendo em casos simples, e que as estimativas fazem sentido quando pensamos nas distribuições, ou quando olhamos um histograma.

Bem eu ainda não tenho os culhões para deduzir tudo isso sozinho, essa conta tem no livro Ecological models and Data in R, que é um dos livros mais legais que ja li, apesar de não entender 3/4 do que está escrito, foi o primeiro livro que vi tentando explicar estatística e modelagem, de uma perspectiva mais biológica que bolinhas pretas e brancas em potinhos, alias eu sou um grande fa do autor, o Ben Bolker, uma das imagens do cabeçalho do blog é a foto dele, se prestar atenção nesse nome, direto vai encontrar ele espalhado pelos pacotes do R. Além de tudo acho um link legal para ver onde o calculo entra na bio-estatística, e como a gente usa calculo muito mais rotineiramente que pensamos, gostamos dele ou não. E ainda bem que a distribuição binomial tem somente um parâmetro, para facilitar entender a derivação, que é bem complicada, mas com algum esforço da para entender.

Bem ficamos por aqui hoje. Feliz ano novo a todos que passarem por aqui, e que tenhamos um 2014 melhor que 2013.

Referência:

Bolker, B. 2008 – Ecological Models and Data in R 408pp Princeton University Press

Noções de notação matemática para o teste t

Opa, após alguma semanas parado apanhando pra conseguir voltar o blog, surge um novo post.

Vamos ver aqui um pouquinho de como é a notação desses modelos que a gente faz, como escrever eles com aquelas letrinhas bonitinhas da matemática que normalmente confunde as pessoas.
A parte boa de tentar entender essa notação, é que ela pode ajudar nosso entendimento dos modelos, assim como melhorar nossa intuição do que deve estar acontecendo.

Bem teste-t é uma análise bem básica, que a gente ja viu milhares de vezes, como aqui e aqui, para ver diferença de médias.
Mas vemos ver de novo ele para os seguintes dados.

massa<-c(6,8,5,7,9,11)
regiao<-factor(c("Campo Grande","Campo Grande","Campo Grande","Campo Grande","Dourados","Dourados"))

Podemos fazer um gráfico rápido para facilitar o entendimento.

01

Então temos um fator, de dois níveis, que são as cidades, campo grande e dourados, e temos a massa de indivíduos de alguma espécie qualquer.

No R, fazer o teste t significa digitar

t.test(massa~regiao,var.equal=T)

E temos como resposta…

Two Sample t-test data: massa by regiao t = -3.0551, df = 4, p-value = 0.03784 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -6.6808175 -0.3191825 sample estimates: mean in group Campo Grande mean in group Dourados 6.5 10.0

Certo, mas isso na verdade a gente está fazendo um modelo linear, é algo como um caso especial, não um caso especial, simplesmente um caso de modelo linear com 1 fator de 2 níveis, falar um caso significa que queremos ver o teste-t dentro da teoria que usamos para fazer todos os modelo lineares ou modelos gerais linearizados.

Bem um modelo linear no R é feito simplesmente usando o comando lm.

lm(massa~regiao,var.equal=T)
Call: lm(formula = massa ~ regiao) Coefficients: (Intercept) regiaoDourados 6.5 3.5

Bem se a gente calcular a média das massas para as duas cidades a gente vai ver que:

> aggregate(massa,list(regiao),mean) Group.1 x 1 Campo Grande 6.5 2 Dourados 10.0

Agora esse 6.5, que está indicado como intercept da para saber o que é, é a media de campo grande, agora e esse 3.5 de dourados, e aonde estão os valores p?

Primeiro, temos que perceber que antes de tudo temos que estimar parâmetros, para construir um modelo, mas o que diabo é o modelo? Agora quando a gente tentar escrever a notação matemática, vamos conseguir talvez entender melhor, visualizar o quadro maior.

Quando a gente fala modelo, a gente ta falando em fazer uma função matemática, que representa a situação que a gente observa, qual a melhor função matemática para descrever o que a gente ta vendo?
Uma função exata não da para fazer quase sempre, um função probabilística sempre da
O porque não conseguimos fazer uma função matemática perfeita para os dados? Porque existe quase sempre um componente estocástico, um componente chamado acaso, que a gente pode esperar algo, mas não temos como saber exatamente o que. Mesmo que você tenha todos os dados dos mundo, a posição exata de cada atomo, sua carga, toda a informação, ainda sim não da para prever o tempo perfeitamente. Perfeitamente não, mas razoavelmente bem sim.

Ta ficando confuso, mas muita calma nessa hora.

Matematicamente, o que estamos fazendo aqui é uma função assim:

massa_i = \alpha + \beta \cdot regiao_i + \epsilon_i

Poxa, a resposta então é a massa, ali a gente tem a equação de reta, a+bx e um erro, esse sendo representado por esse simbolo epsilon (\epsilon_i).
Os i nas são referentes aos índices, índices do que? Dos casos, das amostras que temos, vamos olhar a massa, no R a gente tem:

> massa [1] 6 8 5 7 9 11 > massa[1] [1] 6 > massa[3] [1] 5

Então a gente tem 6 amostras, o i é o índice, quando o índice = 1, o valor de massa que observamos é 6, quando o índice=3, o valor observado é 5, isso é o i, índice de posição, o mesmo vale para a região.

Beleza, mas a notação completa do modelo é assim;

massa_i = \alpha + \beta \cdot regiao_i + \epsilon_i
\epsilon_i \sim Normal(0,\sigma^2)

Ou as vezes assim:

massa_i = Normal(\alpha + \beta \cdot regiao_i,\sigma^2)

Então a primeira linha é a equação determinística, mas nela ta somado um número que vem de uma distribuição normal com média zero, a segunda linha, no segundo caso a gente so juntou tudo isso numa linha, mas o primeiro jeito eu acho que é mais simples de entender.

Mas como a gente ve isso em relação aos nossos dados? Assim.

  6  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_1 \\  8  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_2 \\  5  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_3 \\  7  = \alpha \cdot 1 + \beta \cdot 0 + \epsilon_4 \\  9  = \alpha \cdot 1 + \beta \cdot 1 + \epsilon_5 \\  11 = \alpha \cdot 1 + \beta \cdot 1 + \epsilon_6 \\

Essencialmente esse é o meu modelo, eu to tentando ver ele em relação a cada amostra de dados que tenho.
Quando a gente estima parâmetros, o que eu quero saber é qual valor eu posso colocar para alpha(\alpha) e para beta(\beta) para solucionar esse monte de equação ai. So que lembre-se que o epsilon (\epsilon) é um valor que vem de uma distribuição normal, uma distribuição de eventos estocásticos, dados ao acaso, então ele pode ser qualquer valor, eu posso colocar qualquer número pra ele, qualquer número mesmo, ou seja, eu sempre vou conseguir solucionar todas essa equações ae em cima, para qualquer valores de alpha e beta. Por exemplo , eu posso colocar o valor de alpha como 100, beta como um milhão e epsilon como -94 e ta resolvido.

 6  = 100 \cdot 1 + 1000000 \cdot 0 + (-94)

Espero não ter errado, mas 100 vezes 1 da 100, dai um milhão vezes zero da zero, e tirando 94 da 6. So que essa solução é ruim, porque esses valores de alpha e beta fazem a gente ter que usar valores muito grandes para epsilon, mais precisamente, muito diferentes para cada equação. Imagina la em baixo quando o um milhão é multiplicado por 1, a gente vai precisar de um epsilon de mais de um milhão la, enquanto usamos 94 na primeira equação, ou seja para conseguir valores assim precisamos de um desvio muitooooo grande para o epsilon, e a pergunta é, quais valores eu ponho em alpha e beta para usar o menor desvio possível no epsilon. Pra variar o mínimo possível os valores que vou ter que colocar no epsilon para satisfazer todas essas equações ai de cima.

Certo, até aqui tudo bem, mas como isso está no R?

Bem o calculo não é feito com aquele monte equação, tem um jeito mais simples de fazer isso no R, que é usando álgebra linear(ou não, dependendo do tamanho do seu ódio para com álgebra linear). Mas a gente usa operações com matrizes.

Fica assim a conta para o nosso exemplo.

 \left(\begin{array}{c} 6 \\ 8 \\ 5 \\ 7 \\ 9 \\ 11 \end{array} \right) = \left(\begin{array}{cc} 1&0 \\ 1&0 \\ 1&0 \\ 1&0 \\ 1&1 \\ 1&1 \end{array} \right) \ast  \left(\begin{array}{c} \alpha \\ \beta \end{array} \right) +  \left(\begin{array}{c} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \\ \epsilon_4 \\ \epsilon_5 \\ \epsilon_6 \end{array} \right)

Se a gente resolver essa multiplicação e a soma de matrizes e vetores, a gente vai resolver exatamente aquele conjunto de equaçãozinha la em cima.

Agora é so pensar, antes do igual, é nosso vetor chamado massa:

> massa [1] 6 8 5 7 9 11

Essa é a primeira matriz, antes do sinal de igual, a matriz de 1 coluna e 6 linhas(o número de linhas é o número de amostras) de resposta.

A segunda matriz, é a matriz do modelo. Pra ver ela é so ver como é guardado os dados da região.
Ele é um fator, que tem 2 componentes, os níveis e como eles estão distribuídos.

Com o comando levels a gente ve os níveis no R

> levels(regiao) [1] "Campo Grande" "Dourados" > as.integer(regiao) [1] 1 1 1 1 2 2

E temos números que representam onde eles estão, é assim porque a gente faz conta com número, para nos os nomes tem significado, mas para a álgebra das matrizes não.

Mas é so unir essas duas informações que temos nossa informação devolta

> levels(regiao)[as.integer(regiao)] [1] "Campo Grande" "Campo Grande" "Campo Grande" "Campo Grande" "Dourados" "Dourados"

Bem com essa informação, o R constrói aquela matriz que de 0 e 1 do modelo

> model.matrix(~regiao) (Intercept) regiaoDourados 1 1 0 2 1 0 3 1 0 4 1 0 5 1 1 6 1 1 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$regiao [1] "contr.treatment"

E assim ja temos a nossa função, agora é só minimizar os quadrados para obter a resposta. Isso é os outros dois componente la de cima.

O comando lm ta achando aquele vetor que tem o alpha e o beta:

> lm(massa~regiao) Call: lm(formula = massa ~ regiao) Coefficients: (Intercept) regiaoDourados 6.5 3.5

Se você não se ligou ainda, essas são as médias, do que todas são dadas a partir do intercepto, ou seja a média de campo grande é 6.5, e a média de dourados é 10, que é o 6.5+3.5, aqui é quando o beta é multiplicado por 1 na matriz la em cima.

E com os menores resíduos possíveis, resido é os números que temos que atribuir a epsilon, esses aqui:

> resid(lm(massa~regiao)) 1 2 3 4 5 6 -0.5 1.5 -1.5 0.5 -1.0 1.0

Se a gente substituir esses valores la, temos a solução de todas as equações, usando o menor desvio possível.

Bem, então o nosso exemplo é assim:

 \left(\begin{array}{c} 6 \\ 8 \\ 5 \\ 7 \\ 9 \\ 11 \end{array} \right) = \left(\begin{array}{cc} 1&0 \\ 1&0 \\ 1&0 \\ 1&0 \\ 1&1 \\ 1&1 \end{array} \right) \ast  \left(\begin{array}{c} 6.5 \\ 3.5 \end{array} \right) +  \left(\begin{array}{c} -0.5 \\ 1.5 \\ -1.5 \\ 0.5 \\ -1 \\ 1 \end{array} \right)

Essa operação ai, é exatamente esse conjunto de equações.

  6  = 6.5 \cdot 1 + 3.5 \cdot 0 + -0.5 \\  8  = 6.5 \cdot 1 + 3.5 \cdot 0 + 1.5 \\  5  = 6.5 \cdot 1 + 3.5 \cdot 0 + -1.5 \\  7  = 6.5 \cdot 1 + 3.5 \cdot 0 + 0.5 \\  9  = 6.5 \cdot 1 + 3.5 \cdot 1 + -1 \\  11 = 6.5 \cdot 1 + 3.5 \cdot 1 + 1 \\

Por exemplo, o primeiro caso, 6.5 vezes 1 da 6.5, mais 3.5 vezes zero da zero, e continuamos com 6.5, tirando -0.5 da 6, que é o nosso valor observado.

E qual diabos é a função ai? Podemos escrever ela assim.

   f(x) = \left\{       \begin{array}{lr}         6.5 + \epsilon & : x = 1\\         6.5 + 3.5 + \epsilon & : x = 2       \end{array}     \right.

Onde x é nada mais que a região. E lembre-se que a gente tem a região representada assim

> regiao [1] Campo Grande Campo Grande Campo Grande Campo Grande Dourados Dourados Levels: Campo Grande Dourados

Mas o computador le assim:

> levels(regiao) [1] "Campo Grande" "Dourados" > as.integer(regiao) [1] 1 1 1 1 2 2

E dai vem o 1 e 2 da equação la em cima, x so pode assumir esses dois valores para essa equação. Agora imagine que para generalizar aquela equação para 3 níveis, como numa anova, é so adicionar uma linha, simples certo?

Agora da para visualizar como o teste-t é um caso especial de anova com apenas 2 níveis.

Bem vamos olhar mais uma vez o resultado de um modelo linear no R

> summary(lm(massa~regiao)) Call: lm(formula = massa ~ regiao) Residuals: 1 2 3 4 5 6 -0.5 1.5 -1.5 0.5 -1.0 1.0 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.5000 0.6614 9.827 0.000601 *** regiaoDourados 3.5000 1.1456 3.055 0.037841 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.323 on 4 degrees of freedom Multiple R-squared: 0.7, Adjusted R-squared: 0.625 F-statistic: 9.333 on 1 and 4 DF, p-value: 0.03784

Bem agora a gente sabe o que são aqueles resíduos, os coeficientes ou parâmetros estimados, e aquele itemzinho chamado “Residual standard error” é o desvio usado no Epsilon.

> sqrt(deviance(lm(massa~regiao))/df.residual(lm(massa~regiao))) [1] 1.322876

A menor quantidade de erro para satisfazer o modelo com esses parâmetros ae.

Assim se não houve-se diferença na massa entre dourados e campo grande, esperaríamos um valor de beta próximo de zero. A nossa equação ficou.

   f(x) = \left\{       \begin{array}{lr}         6.5 + \epsilon & : x = 1\\         6.5 + 3.5 + \epsilon & : x = 2       \end{array}     \right.

Deveria ter a mesma média para x=1 e x=2, o termo de diferença, que é o beta, deveria ser zero.
Como não é, é porque tem essa diferença de 3.5, lembrando que é uma diferença de mais ou menos 3.5, não exatamente, porque o epsilon pode assumir qualquer valor.

Bem, ficamos por aqui, por enquanto. Mas acho que agora talvez esses modelos façam mais sentido para todo mundo. Eu me pergunto porque eu nunca vi uma aula as pessoas mostrando esses modelos desse jeito, acho que eu entenderia melhor, mas ta bom, ta ruim mas ta bom.

Bem apesar do pau no site, o repositório recologia sempre esteve la. E é isso, até o próximo post.

massa<-c(6,8,5,7,9,11)
regiao<-factor(c("Campo Grande","Campo Grande","Campo Grande","Campo Grande","Dourados","Dourados"))
 
plot(massa~regiao,frame=F)
 
regiao
 
levels(regiao)
 
as.integer(regiao)
 
levels(regiao)[as.integer(regiao)]
 
t.test(massa~regiao,var.equal=T)
 
lm(massa~regiao)
 
resid(lm(massa~regiao))
 
model.matrix(~regiao)
 
summary(lm(massa~regiao))
 
model.matrix(~regiao-1)
 
summary(lm(massa~regiao-1))
 
aggregate(massa,list(regiao),mean)
 
sqrt(deviance(lm(massa~regiao))/df.residual(lm(massa~regiao)))

Um exemplo do algoritimo de Metropolis–Hastings aplicado a distribuição binomial

E la vamos nos denovo olhar o algoritmo Metropolis–Hastings

Bem depois de olhar o exemplo do como estimar o valor de pi usando o método de monte carlo, devemos progredir um pouco no entendimento desse negócio de mcmc.

Ja tivemos posts aqui e aqui, onde olhamos como o algoritmo funcionava para fazer uma regressão linear.

Mas vamos tentar um exemplo mais simples para estimar a probabilidade encontrar ou não aranhas num arbusto.

Certo então existe uma probabilidade \theta (letra grega theta) na qual eu encontro ou não as aranhas nos arbustos.

Então eu fui la e olhei 14 arbustinhos e vi o seguinte:

> meusdados [1] 0 1 1 0 1 0 0 1 0 1 1 1 1 1

Ok 1 significam que nesse arbusto eu encontrei uma aranha, e 0 significa que não tinha aranha.

Eu posso estimar a verosimilhança para várias probabilidades, para vários valores de \theta, ou mais precisamente, o p(D|\theta) onde D são meus dados

likelihood <- function( theta , data ) {
  z <- sum( data == 1 )
  N <- length( data )
  pDataGivenTheta <- theta^z * (1-theta)^(N-z)
  pDataGivenTheta[ theta > 1 | theta < 0 ] <- 0
  return( pDataGivenTheta )
}

Ou seja qual a chance de determinado \theta dado meus dados (p(D|\theta)). Aqui tem uma revisão rápida de probabilidade condicional.

Bom mais fácil que ficar nisso é fazer um gráfico para entender o que essa função ta fazendo.

01

Então temos a como são as chances dado um valor de \theta, ou seja olhando os dados la tivemos 9 arbustos com aranhas, dos 14 que olhamos, sabemos que o maximum likelihood nesse caso é 9/14, ou seja o pico do likelihood é quanto \theta esta num valor igual a 0.64, ou seja um \theta de 0.64 é o mais provável de ter produzido esses resultado, dado os dados que temos, veja que se você pegar um valor de vamos supor 0.2, quase não tem como achar 9 arbustos com aranhas de 14 arbustos verificados, teria que ser uma cagada muito grande.

Mas lembramos que em inferência bayesiana, ponderamos nossos dados, ou seja multiplicamos esse likelihood por um prior.

Então fazemos uma função para gerar esse prior.

prior <- function( theta ) {
  prior <- dbeta(theta,1,1)
  #exemplo de um prior bimodal
  #prior < dbeta( pmin(2*theta,2*(1-theta)),2,2 )
  prior[ theta > 1 | theta < 0 ] <- 0
  return( prior )
}

Aqui, o prior ta definido ali como uma distribuição beta com parâmetros a=1 e b=1. Isso vai dar uma gráfico assim:

02

Veja que todas as possibilidades estão em 1, ou seja eu estou atribuindo uma chance igual para todos os valores possíveis de \theta.

Veja que podemos mudar nossa crença inicial mudando os parâmetros a e b para outros valores, se você não sabe o que acontece com diferentes valores, é mais ou menos isso, veja alguns exemplos:

01

Certo, mas para facilidade do exemplo vamos assumir aquele prior e não ficar discutindo muito isso agora.

Agora se lembrarmos o teorema de bayes usado na inferência bayesiana:

p(\theta|D) = {p(D|\theta) * p(\theta) \over {\int p(D|\theta)p(\theta)}}

Então o que falta é sabe qual o valor desse divisor. Na verdade para esse exemplo da pra resolver analiticamente como feito aqui ou aproximar como aqui, mas vamos usar o algorítimo de Metropolis–Hastings para construir a probabilidade posterior.

A ideia é a mesma de quando estimamos o valor de pi. Relembrando os passos:

1. Defina um domínio de todas as possibilidades 2. Gere aleatoriamente entradas de uma distribuição de probabilidade que cobre o domínio. 3. Faça uma computação determinística das entradas 4. Reuna os resultados

Então o domínio são as probabilidades de achar aranha em um arbusto, que são valores entre 0 e 1.

Eu tenho que gerar valores ao acaso de probabilidade, que cubram todas as possibilidades e fazer uma computação determinística para levar ele em conta ou não, no caso do pi o que a gente gerava um ponto ao acaso dentro de um quadradinho e tínhamos um círculo dentro desse quadrado, ai olhávamos se ele estava ou não dentro do círculo, para no final fazer uma razão.

Aqui o que temos são valores possíveis de \theta e eu quero saber qual gerou os dados, bom isso eu nunca vou saber, mas eu quero ter uma estimativa, uma distribuição que se o valor de \theta seja bem razoável, ele tem que ser alto, senão ele é baixo.
Talvez eu esteja falando bobeira aqui, mas o que acontece é o seguinte, no caso do Metropolis–Hastings, ficamos escolhendo um valor de \theta ao acaso, e colocamos ele ou não na nossa cadeia numa probabilidade que depende de o quão bom o novo \theta é dividido pelo \theta anterior.

É difícil explicar isso, acho que entender também, mas fazemos uma função assim:

targetRelProb <- function( theta , data ) {
  targetRelProb <-  likelihood( theta , data ) * prior( theta )
  return( targetRelProb )
}

Essa função multiplica o likelihood pelo prior, dai geramos uma valor de theta vai ter essa multiplicação do valor novo de theta pela multiplicação com o antigo.

Vamos olhar alguns exemplos, bem eu usei 0.7 como chance para gerar os dados. Então vamos supor que o valor anterior na cadeia é 0.7 e fosse sugerido fosse de 0.4, qual a chance de aceitar esse valor e colocar na cadeia?

> targetRelProb(0.4,meusdados)/targetRelProb(0.7,meusdados) [1] 0.2078775

20%, baixinho certo, ou seja eu quero “visitar” esse valor pouco, quero poucas ocorrências dele na cadeia, assim no histograma da distribuição posterior ele vai ter uma barrinha baixinha, mas impossível ver um valor desse? Não uai, mas é pouco provável.

Agora vamos mais pra perto, se propusermos um valor de 0.5.

> targetRelProb(0.5,meusdados)/targetRelProb(0.7,meusdados) [1] 0.6224313

Olha ai, temos 60% de chance de enfiar esse valor na cadeia, porque ta mais pertinho do valor real. E quanto mais próximo vamos do valor real, agora com 0.55.

> targetRelProb(0.55,meusdados)/targetRelProb(0.7,meusdados) [1] 0.8666388

Temos mais e mais chance de enfiar esses valores pra dentro da cadeia.
Agora até aqui eu sempre tava dividindo com um theta de 0.7, mas e se sei la estamos no começo da cadeia, e estamos em valores bem distantes desses. Por exemplo estamos em 0.2 e temos uma proposta de ir para 0.1

> targetRelProb(0.1,meusdados)/targetRelProb(0.2,meusdados) [1] 0.003519595

Meu, 0.2 ta longe pra caramba das probabilidade real, mas 0.1 é muito mais longe, então a chance de aceitar esse valor é muito muito baixa. Agora num caso ao contrário, estamos em 0.9, que é longe e recebemos uma proposta de 0.8, que é bem mais próximo, qual a chance?

> targetRelProb(0.8,meusdados)/targetRelProb(0.9,meusdados) [1] 11.08606

Puxa, um valor muito maior que 1, veja que isso pode acontecer, na verdade aqui o que fazemos e com certeza aceitar o valor proposto, como se estamos em 0.5 e recebemos a proposta de 0.6

> targetRelProb(0.6,meusdados)/targetRelProb(0.5,meusdados) [1] 1.690757

Mas para fugir desses valores e ficar com valores sempre entre 0 e 1 na cadeia temos que fazer umas correções nos valores, mas basicamente é isso o algoritimo, gera um theta novo e aceita ele ou não, um bocado de vezes, mas um bocado são muitas e muitas vezes mesmo.

Então o algoritimo vai ser isso aqui:

for ( t in 1:(Ncadeia-1) ) {
    valoratual = cadeia[t]
    valorproposto = rnorm( 1 , mean = 0 , sd = 0.1 )
    chanceaceitar = min( 1,
        targetRelProb( valoratual + valorproposto, meusdados ) /
        targetRelProb( valoratual , meusdados )
        )
    if ( runif(1) < chanceaceitar ) {
        cadeia[ t+1 ] <- valoratual + valorproposto
        } else {
            cadeia[ t+1 ] <- valoratual
	}
}

A gente salva numa variável o valor atual, na prática a gente gera uma quantidade, que é o valor proposto, que é uma quantidade para mudar o valor de theta atual, dai calculamos como comentamos acima a chance de aceitar esse novo valor, e se ele for aceito colocamos ele na cadeia, senão repetimos o último valor. E milagrosamente vamos visitar cada valor possível de theta um número de vezes proporcional ao valor de theta que gerou os dados. Na verdade não é milagre, essa estabilidade é prevista para cadeias de markov dado a chance de transição entre estados, mas pra mim é mágico como funciona 🙂

Então ao usarmos nosso algoritimozinho, o que vemos como resultado?

03

Uma boa estimativa de valor real de theta, o intervalo de 95% da distribuição cobre o valor real de theta.

Mas temos uma distribuição continua de thetas, não quero ficar olhando as coisas na forma de histograma, podemos olhar a mesma distribuição na forma de densidade contínua.

04

Lembrando que aqui a gente ja arrancou o início da cadeia, que como começamos num lugar qualquer, traz estimativas ruim. Vamos dar uma olhadinha na cadeia, para o início e um pedaço la na frente.

05

Como começamos num valor arbitrário qualquer, começamos num valor de theta que não tem nada a ver com o valor real e vamos caminhando ao valor próximo do real, ai chegando la ficamos rodeando em volta dele. Veja que o tamanho dos passos aqui depende do valor proposto, se propomos valores pouco diferentes do anterior, vamos andando devagarinho na cadeia, agora poderíamos propor passos mais largos, mais ai também, assim como poderíamos chegar num valor próximo do real mais rápido, poderíamos derrepente sair dele pra muito longe. Esse tipo de ajustamento, como o tamanho do passo, como propor uma esse passo que é a parte difícil, aqui tudo funciona mas em algum momento alguém deve ter ficado testando possibilidades e ficou frustado em como as coisas não funcionavam. Outra coisa é que pode ser proposto valores maiores que 1 ou negativos, menor que zero, veja que ficamos mudando esses valores para os extremos, 0 e 1, mas eu não sei qual o impacto de fazer isso se o valor de theta fosse 1 por exemplo, mas todas essas coisas influenciam no funcionamento e eficácia do uso de MCMC. Nisso pacotes como OpenBugs ou Stan salvam a vida, ao possibilitar a gente abstrair dessa coisas e focar em entender melhor nossas observações da natureza. Mas ainda sim é importante entender pelo menos um pouco do de tudo que estamos fazendo.

Bem acho que ja esta bom por aqui, depois olhamos mais sobre MCMC. Esse script como outros do blog estão sempre disponíveis no repositório do github aqui.

Referência:
John K. Kruschke 2011 Doing Bayesian Data Analysis: A Tutorial with R and BUGS. Academic Press / Elsevier

#Adaptado de uma função do livro Bayesian Data Analysis: A Tutorial with R and BUGS
#Semente para replicar os resultado do post.
set.seed(51)
 
#Gerando uma amostra de dados
probabilidade<-0.70
meusdados<-rbinom(n=14, size=1, prob=probabilidade)
 
#Definindo a função de likelihood binomial.
likelihood <- function( theta , data ) {
  z <- sum( data == 1 )
  N <- length( data )
  pDataGivenTheta <- theta^z * (1-theta)^(N-z)
  pDataGivenTheta[ theta > 1 | theta < 0 ] <- 0
  return( pDataGivenTheta )
}
 
plot(seq(0.01,0.99,0.01),likelihood(seq(0.01,0.99,0.01),meusdados),frame=F,
     xlab="Probabilidade",ylab="Likelihood",type="l")
 
#Função para definir o prior p(D),
prior <- function( theta ) {
  prior <- dbeta(theta,1,1)
  #exemplo de um prior bimodal
  #prior < dbeta( pmin(2*theta,2*(1-theta)),2,2 )
  prior[ theta > 1 | theta < 0 ] <- 0
  return( prior )
}
 
plot(seq(0.01,0.99,0.01),prior(seq(0.01,0.99,0.01)),frame=F,
     xlab="Probabilidade",ylab="Likelihood",type="l")
 
#Função auxiliar para ajudar a calcular a probabilidade relativa
targetRelProb <- function( theta , data ) {
  targetRelProb <-  likelihood( theta , data ) * prior( theta )
  return( targetRelProb )
}
 
# Tamanho da cadeia
Ncadeia = 100000
 
#Inicializando um vetor para salvar os resultados do MCMC
cadeia = rep( 0 , Ncadeia )
 
# Iniciando o MCMC com um valor arbritario
cadeia[1] = runif(1,min=0,max=1)
 
# Especificando um valor de burnning, ou seja o quanto do inicio da cadeia vamos jogar fora
# Aqui vamo jogar fora 10%
burnIn = ceiling( .1 * Ncadeia )
 
# Vamos contar quantas vezes aceitamos ou não os valores propostos
nAceito = 0
nRejeitado = 0
 
# Agora iniciando o "Random Walk"
for ( t in 1:(Ncadeia-1) ) {
    #Valor Atual
    valoratual = cadeia[t]
    #Propomos uma quantidade para alterar esse valor
    valorproposto = rnorm( 1 , mean = 0 , sd = 0.1 )
    #Computamos a chance de aceitar a mudança
    #Essa chance é igual a Likelihood*prior com o novo valor de theta
    #dividido pele Likeliood*prior antigo
    chanceaceitar = min( 1,
        targetRelProb( valoratual + valorproposto, meusdados ) /
        targetRelProb( valoratual , meusdados )
        )
    #Agora aqui a gente ve na sorte se aceitamos ou não o novo valor.
    #Veja que temos uma chance de fazer a mudança ou manter o ultimo valor
    #Se a gente sortei um número qualquer de 0 a 1 e compara se ele é
    #menor que a chance de aceitar, dessa forma mudamos numa probabilidade
    #da chance de mudar
    if ( runif(1) < chanceaceitar ) {
        #Se aceitarmos mudamos ele na cadeia
        cadeia[ t+1 ] <- valoratual + valorproposto
        #So vamos começar a contar depois do burnIn que escolhemos
        #Afinal, com o começo arbitrario, esperamos aceitar bastante
        #a não ser que escolhemos um valor proximo do final no começo
        if ( t > burnIn ) { nAceito <- nAceito + 1 }
 
        } else {
            # Quando não aceitamos o valor proposto mantemos o anterior
            cadeia[ t+1 ] <- valoratual
            # E incrementamos aqui apos o BurnIn, como acima
            if ( t > burnIn ) { nRejeitado <- nRejeitado + 1 }
	}
}
 
#Agora vamos olhar como ficou a cadeia final, ou seja, tirando os valores iniciais
cadeiafinal <- cadeia[ (burnIn+1) : length(cadeia) ]
 
#Alguns Gráficos
 
#-----------------------------------------------------------------------
HDIofMCMC = function( sampleVec , credMass=0.95 ) {
    sortedPts = sort( sampleVec )
    ciIdxInc = floor( credMass * length( sortedPts ) )
    nCIs = length( sortedPts ) - ciIdxInc
    ciWidth = rep( 0 , nCIs )
    for ( i in 1:nCIs ) {
        ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
    }
    HDImin = sortedPts[ which.min( ciWidth ) ]
    HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
    HDIlim = c( HDImin , HDImax )
    return( HDIlim )
}
 
 
#------------------------------------------------------------------------
intconf<-HDIofMCMC(cadeiafinal)
 
hist(cadeiafinal,ylab="Frequência na cadeia final",xlab="Probabilidades",xlim=c(0,1))
lines(c(intconf[1],intconf[2]),c(500,500),lwd=4,lty=3,col="red")
abline(v=mean(cadeiafinal),lwd=4,col="red")
abline(v=probabilidade,lwd=2,col="blue")
 
legend("topleft",lty=c(3,1,1),lwd=c(4,4,2),bty="n",cex=0.8,col=c("red","red","blue"),
       legend=c("Intervalo de confiança","Média das estimativas","Probabilidade usada nas estimativas"))
 
plot(density(cadeiafinal),ylab="Frequência",xlab="Probabilidades",xlim=c(0,1),frame=F)
lines(c(intconf[1],intconf[2]),c(0.25,0.25),lwd=3,lty=2,col="red")
abline(v=mean(cadeiafinal),lwd=1,col="red")
abline(v=probabilidade,lwd=1,col="blue")
 
legend("topleft",lty=c(3,1,1),lwd=c(2,1,1),bty="n",cex=0.8,col=c("red","red","blue"),
       legend=c("Intervalo de confiança","Média das estimativas","Probabilidade usada nas estimativas"))
 
layout(matrix(c(1,2,3,4),ncol=2))
hist(cadeia[1:100],ylab="Frequência na cadeia",xlab="Probabilidades",main="100 primeiros valores",xlim=c(0,1))
plot(0,0,xlim=c(0,100),ylim=c(0,1),type="n",frame=F,main="Cadeia",
     xlab="Posição na cadeia",ylab="Probabilidade")
points(cadeia[1:100],type="l")
 
hist(cadeia[10000:10500],ylab="Frequência na cadeia",xlab="Probabilidades",main="Do 10000 ao 10500",xlim=c(0,1))
plot(0,0,xlim=c(0,500),ylim=c(0,1),type="n",frame=F,main="Cadeia",
     xlab="Posição na cadeia",ylab="Probabilidade",xaxt="n")
axis(1,at=seq(0,500,by=50),labels=seq(10000,10500,by=50))
points(cadeia[10000:10500],type="l")
 
# Evidencia para o modelo, p(D).
 
media<-mean(cadeiafinal)
desvio<-sd(cadeiafinal)
 
a =   media   * ( (media*(1-media)/desvio^2) - 1 )
b = (1-media) * ( (media*(1-media)/desvio^2) - 1 )
 
wtdEvid = dbeta( cadeiafinal , a , b ) /
    ( likelihood( cadeiafinal , meusdados ) * prior( cadeiafinal ) )
pData = 1 / mean( wtdEvid )
 
format(pData,scientific = F)