Teste t de Poisson

O titulo de post é um nome bem esquisito, e errado também, teste t usa a estatística t, então não da para ter um teste t de Poisson exatamente. Mas a gente acaba associando alguns nomes a situações.
Teste t para avaliar um fator de dois níveis, anova para um fator de vários níveis, ancova para variáveis continuas e variáveis discretas e assim vai.
A questão é que todos esses métodos estão na verdade dentro de uma categoria, existe uma teoria que unifica todos esses métodos, esses e outros, como aqueles que já comentamos usando distribuição binomial, ou mesmo com distribuição de poisson como vemos aqui são todos casos da mesma “estratégia” em estatística aplicada, se meus termos aqui não são ruins. Mas em 1972, John Nelder e Robert Wedderburn unificaram mais ou menos tudo que falamos ali em cima na ideia geral de Modelos gerais linearizados.. Eles mostraram que todas essas técnicas, representadas como tipos diferentes de analises, incluindo chi-quadrados até que a gente ve muito em génetica, podiam todos serem representados como casos especiais de modelos lineares. Quando eu fiz estatística na faculdade, na graduação, eu nunca fui apresentado a essa ideia que tudo é “modelo linear”, que tudo esta interligado e por isso muitas ideias são comuns, teste t e anova eram coisas distintas, técnicas que se você tem a situação a você aplica isso, se tem a situação b, aplica aquilo e se por ventura tem a situação c você se ferrou, porque não tem um teste pronto pra você. E isso gera o problema, que as vezes as pessoas abandonam a ideia, pressupostos da biologia por exemplo, para forçar os dados serem aplicáveis a um teste e temos que lidar com um monte de problemas decorrentes disso. Além é claro de isso muitas vezes diminuir o avanço cientifico, se é que alguém ainda se importa com isso.

Mas então, um modelo geral linear vai ter mais ou menos três componentes.

1. Uma distribuição estatística, usada para descrever a parte aleatória da resposta do modelo, o eixo y dos gráficos que a gente normalmente faz. Essa é a parte estocástica do sistema que estamos tentado descrever.
2. Uma função de ligação, que é aplicada ao valor esperado, um E(y).
3. E um preditor linear, os a+bx, intercepto e inclinação, as variáveis preditoras que na maioria das vezes estamos interessados. se a função de ligação é a função g, então um g(E(y)). Essa é a parte determinística do modelo, veja que aqui entra qualquer função matemática na verdade, mas essa ai é bem comum.

Nesse sentido a regressão linear que a gente faz é o seguinte.

1. Distribuição: y_i \sim Normal(\mu_i , \sigma^2)
2. Função de ligação: Identidade, \mu_i = E(y_i) que é o preditor linear.
3. Preditor linear: \alpha + \beta \cdot x_i

Se você pensa que nunca fez isso, veja que todas as simulações que fazemos aqui a gente meramente vai do item 3 pro item 1, a gente gera os valores determinísticos, transforma pela função de ligação e depois gera a parte aleatória dos dados, e depois roda uma analise para ver se recuperamos os valores do preditor linear.

Mas e um modelo usando Poisson?

Ai teremos o seguinte:

1. Distribuição: C_i \sim Poisson(\lambda_i)
2. Função de ligação: \log{(\lambda_i)} ou \log({E(C_i)})
3. Preditor linear: \alpha + \beta \cdot x_i

Viu, a mesma ideia, só mudamos a função de ligação, que antes era a identidade, que nada mais é que o mesmo valor do número, agora usamos o log, log na base natural, que é o número do euler la, que a gente ja esbarrou por ele aqui.

Bem agora a analise, primeiro vamos entender os dados.

Vamos supor que queremos estudar o número de aranhas que encontramos em arbustos. Então a gente pega um arbusto inteiro, poe ele num sacão e leva para o laboratório, la vasculhamos ele e contamos quantas aranhas conseguimos achar. Particularmente os arbustos são de uma espécie de planta que ocorre tanto na zona rural como na zona urbana, em terrenos baldios. Mas nos achamos que aranhas não conseguem dispersar para dentro da zona urbana, logo esperamos que elas sejam restritas a áreas rurais, certo? Claro que uma ou outra aranha vai acabar chegando na cidade, ainda mais que estamos falando “aranhas” sem considerar as espécies delas. Mas de modo geral esperamos menos aranhas na cidade, ou mais claramente falando, esperamos que em média, arbustos na cidade tenham menos aranhas, ou seja em média as contagens de aranhas serão menor na cidade quando comparado a zona rural.

Então temos a pergunta, existe diferença no número de aranhas entre a área urbana e a área rural?

Bem sabendo o que queremos saber, fica fácil, é so pegar arbustos na cidade e na área rural e contar a aranhas e ver se existe diferença nas contagens de aranhas em cada área. Mas porque não usar um teste t normal? Porque primeiro aranhas so vem inteiras, não existe meia aranha nem um quarto de aranha, como a distribuição normal é continua, ela vai prever uma variável continua, o que é ruim. Mas além disso, a distribuição normal, se eu ver um número muito baixo de aranhas, algo que deixe a média perto de zero, eu vou prever números negativos de aranhas em arbusto, o que não pode existir, se num local existem menos aranhas, eu tenho que ver mais zeros e não números negativos. Esses dois problemas são solucionados usando a distribuição de Poisson, existem mais coisas a se considerar, mas por agora vamos em frente, outras considerações ficam como cenas dos próximos capítulos.

Então fomos la no mato, coletamos os arbustos e contamos a aranhas, e vamos ver o que conseguimos:

02

Olha ai, pegamos 20 arbustos no total, 10 na zona rural e 10 na zona urbana. Como ja suspeitávamos, na cidade não vimos mais que 4 aranhas por arbusto, enquanto no campo olha la, a gente tinha muito mais, até 8 aranhas. Podemos fazer um box-plot.

01

Certo, são os mesmo dados, mas agora acho que visualmente da para ver outra questão, aqui temos uma variação muito maior dos dados na área rural e isso quebra uma das premissas de modelos lineares em geral, como o teste t, que é homogeniedade de variâncias. Mas não é difícil pensar que se a gente sempre vê em média uma aranha, não tem como ter uma variação igual quando a gente tem muitas, pense que não tem como contar algo para menos que zero então um número pequeno implica de certo modo que só podemos ter uma variação pequena, agora se a gente contar em média 100 aranhas nos arbustos, seria fácil pensar por exemplo que erramos e contamos 90, ou 110, e isso faz um desvio muito grande em relação a uma contagem pequena, o ponto é que o desvio tem uma relação com o quanto a gente conta, ou seja com o \lambda, mas na distribuição de poisson é exatamente assim que funciona, o \lambda determina a média e a variação dos dados, no caso quanto maior a contagem média, maior a variação. Temos como mexer um pouco nessa relação, mas isso é outro post, por agora já da para ter uma noção de outra característica do modelo que estamos fazendo, que é em modelos usando distribuição de poisson, a média e o desvio estão relacionados entre si.

Mas podemos ver as médias das áreas também, para ver como ficaram.

> aggregate(C, by = list(x), FUN = mean) Group.1 x 1 Urbana 2.4 2 Rural 5.3

E a média da área rural é quase o dobro da média na área urbana. Agora podemos tentar definir isso como parâmetros de um modelo linear generalizado usando a distribuição de poisson.

Primeiro temos que definir o modelo na linguagem bugs.

# Definindo o modelo na linguagem Bugs
sink("Poisson.t.txt")
cat("
model {
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i]
# Ajuste
    residuo[i] <- (C[i] - lambda[i]) / sqrt(lambda[i])
}
}
",fill=TRUE)
sink()

Além do modelo em si, vamos salvar os resíduos.

Agora se pararmos um momento para olhar o likelihood, podemos ver os itens que definimos la em cima.
O modelo é:

# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i]
}

E olha aqui, temos um loop para processar todos as amostras, e dentro dele temos:

1. A distribuição de poisson

C[i] ~ dpois(lambda[i])

2. A função de ligação

log(lambda[i])

3. A parte determinística do modelo.

alpha + beta *x[i]

Esta numa linguagem de computador, mas é exatamente o que definimos la em cima.
Mas agora com o modelo pronto, definimos as configurações do MCMC e rodamos ele.
Com os resultados em mão, antes de fazer qualquer inferência, primeiro de tudo a gente olha se as cadeias convergiram. Existem muitas formas de ver isso, mas primeiro podemos olhar para os valores de Rhat, eles são simples de interpretar, porque é so buscar se estes estão inferiores a 1.1, o valor limite, acima disso temos que pensar bem como prosseguir, mas nesse caso…

03

Todos menores que 1.1.

Podemos fazer um teste rapido, e simplesmente testar se alguém ficou acima de 1.1

which(out$summary[,8] > 1.1)

Certo, outra coisa, podemos olhar agora os resíduos, estes tem que ter uma distribuição normal, se eles não tiverem uma distribuição normal, significaria que alguma coisa muito importante que não estamos levando em conta esta acontecendo, que tem um efeito claramente mudando muito tudo. Se os resíduos estiveram normal, com média 0, quer dizer que estamos errando igualmente tanto para mais quanto para menos.

04

Olha ai, passamos uma linha no zero, e temos mais ou menos a mesma quantidade de pontos pra cima e para baixo, e nenhuma tendência aparente. Podemos olhar também esses resíduos na forma de um histograma.

05

E apesar de umas barras talvez esquisitas, podemos ver também a densidade em vermelho, que é um “histograma” infinito, mais suave, e vemos bem uma distribuição normal, com média perto do zero, marcado em azul na figura.

Existem também testes para testar a normalidade, como o teste de shapiro-wilks.

> shapiro.test(out$mean$residuo) Shapiro-Wilk normality test data: out$mean$residuo W = 0.9374, p-value = 0.2144

Mas a gente tem que tomar cuidado para não entrar num paradoxo, se a gente começar a testar o teste do teste do teste, a gente não vai parar nunca, e nunca voltaremos para a questão que realmente nos interessa.
Mas pelo menos eu ja estou mais que convencido que o modelo teve um ajuste bom, não tem problemas, então poderíamos olhar as conclusões que ele nos traz, quais inferências podemos fazer.

06

E exatamente o \beta é o que nos interessa, ele é a diferença entre a área rural e a área urbana.
Voltando a nossa pergunta inicial, que era se existia diferença no número de aranhas entre a área urbana e a área rural, um valor de \beta=0 significaria que não há diferença. Como em todas as amostras do MCMC, não temos nenhum 0 nesse caso, nem nada menor que zero estamos fortemente propensos a acreditar que existe uma diferença. Podemos calcular o intervalo de confiança de 95%, ou seja, qual intervalo mais comum de valores e falar que este é maior que zero. Poderíamos sei la, calcular um intervalo de 99%, mas aqui o que interessa é que podemos fazer uma afirmação mais ou menos assim.

“Temos 95% de certeza que deve haver uma diferença no número de aranhas presente em arbustos na zona rural e zona urbana.”

Podemos apresentar então quantas aranhas esperamos encontrar em cada local, e a partir daqui é tentar entender melhor porque existe essa diferença.

07

E assim concluímos nosso primeiro glm usando distribuição de Poisson, já que esse post está ficando grande mais. Mas da para ver como os glms abrem uma grande gama de perguntas que podemos abordar. Espero que não tenha falado muita bobeira no caminho, mas modelos lineares usando distribuição normal, binomial e de poisson cobrem a grande maioria das coisas que lemos comumente na literatura científica relacionada a biologia.

Bem o script está disponível no repositório recologia no github, além de aqui em baixo, e acho que é isso, até o próximo post.

Referencias:
Marc Kéry (2010) Introduction to WinBUGS for Ecologists. Academic Press, Burlington.

set.seed(123)
### Gerando os dados
#Gerando o preditor linear, o item 3 do post.
n.arbustos <- 10
x <- gl(n = 2, k = n.arbustos, labels = c("Urbana", "Rural"))
n <- 2*n.arbustos
 
#Aplicando a função de ligação, exponencial é o contrario de logaritimo,
#estamos no caminho contrario aqui, esse é o item 2 do post
lambda <- exp(0.69 + 0.92*(as.numeric(x)-1))
 
#Adicionando a parte estocastica, o item 1 do post
C <- rpois(n = n, lambda = lambda)
 
#Vamos olhar o que criamos
aggregate(C, by = list(x), FUN = mean)
 
boxplot(C ~ x, col = "grey", xlab = "Tipo de ambiente", ylab = "Contagem de aranhas", las = 1,frame=F)
 
stripchart(C ~ x,xlab = "Tipo de ambiente", ylab = "Contagem de aranhas", las = 1,frame=F,
           method="stack",vertical=T,pch=19,col=1,at=c(1.4,2.1),offset = 0.8)
legend("topleft",pch=19,legend="Um arbusto")
 
 
### Analise Bayesiana
# Definindo o modelo na linguagem Bugs
sink("Poisson.t.txt")
cat("
model {
 
# Priors
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- alpha + beta *x[i]
 
# Ajuste
    residuo[i] <- (C[i] - lambda[i]) / sqrt(lambda[i])
 
}
 
}
",fill=TRUE)
sink()
 
# Juntando os dados para mandar para o Openbugs
dados.bugs<- list(C = C, x = as.numeric(x)-1, n = length(x))
 
# Função geradora de valores iniciais
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))}
inits()
 
# Parametros a estimar
params <- c("alpha", "beta","residuo")
 
# Configurações do MCMC
nc <- 3
ni <- 3000
nb <- 1000
nt <- 2
 
# Iniciando o Gibss Sampler
library(R2OpenBUGS)
out <- bugs(data=dados.bugs, inits=inits, parameters.to.save=params,model.file="Poisson.t.txt",
            n.thin=nt, n.chains=nc, n.burnin=nb,n.iter=ni)
 
 
### Checando convergencia e ajuste do modelo
 
hist(out$summary[,8],col="grey",main="Valores de Rhat")
which(out$summary[,8] > 1.1)
 
plot(out$mean$residuo, las = 1,frame=F,ylab="Resíduos",xlab="Amostras",xlim=c(0,20))
abline(h = 0)
 
hist(out$mean$residuo,xlab="Resíduos",freq = F)
abline(v = 0,lwd=2,col="blue",lty=3)
points(density(out$mean$residuo),type="l",lty=2,lwd=2,col="red")
 
shapiro.test(out$mean$residuo)
 
### Inferencia baseada no modelo que criamos
print(out, dig = 3)
 
#Função do Livro Doing Bayesian Analysis
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 )
}
 
int95<-HDIofMCMC(out$sims.list$beta)
hist(out$sims.list$beta, col = "grey", las = 1, xlab = "Coeficiente para area rural", main = "",freq=F)
points(density(out$sims.list$beta),type="l",lty=2,lwd=3,col="red")
arrows(int95[1],0.25,int95[2],0.25,col="blue",angle =90,lwd=3,lty=3,code=3)
legend("topleft",legend="Intervalo de confiança de 95%",lwd=3,lty=3,col="blue",bty="n")
 
par(mfrow = c(2,1))
int95<-HDIofMCMC(exp(out$sims.list$alpha))
hist(exp(out$sims.list$alpha), main = "Área urbana",col = "grey", xlab = "", xlim = c(0,10), breaks = 20,freq=F)
points(density(exp(out$sims.list$alpha)),type="l",lty=2,lwd=3,col="red")
arrows(int95[1],0.1,int95[2],0.1,col="blue",angle =90,lwd=3,lty=3,code=3)
legend("topright",legend="Intervalo de confiança de 95%",lwd=3,lty=3,col="blue",bty="n")
 
int95<-HDIofMCMC(exp(out$sims.list$alpha + out$sims.list$beta))
hist(exp(out$sims.list$alpha + out$sims.list$beta), main = "Área rural",freq=F,
     col = "grey", xlab = "Número de aranhas esperado", xlim = c(0,10), breaks = 20)
points(density(exp(out$sims.list$alpha + out$sims.list$beta)),type="l",lty=2,lwd=3,col="red")
arrows(int95[1],0.1,int95[2],0.1,col="blue",angle =90,lwd=3,lty=3,code=3)
 
### Analise Frequentista
poisson.t <- glm(C ~ x, family = poisson)
summary(poisson.t)
anova(poisson.t, test = "Chisq")

Vetores no R, representando o espaço vetorial.

Espaço vetorial é algo bem complicado, mas ajuda muito a conseguir visualizar o que está acontecendo por exemplo quando a gente tem um modelo de regressão com duas variáveis que são significativas, como é um gráfico disso, mais complexo que isso é quanto temos três variáveis preditoras, ai ja nem gráfico resolve muito bem, ja que precisa necessariamente um pouco de abstração para pensar, três variáveis preditoras mais uma variável resposta resulta em 4 dimensões, e não da para ver a quarta dimensão, claro que de modo geral, tem gente que deve dizer que vê quarta dimensão, quinta, tudo que é dimensão, depende da cabeça.

Um conceito que eu achei meio complicado, é o de vetores serem linearmente independentes, bem vetores são essas flechinhas, antes que haja duvida porque vale a pensa pelo menos pensar um momento sobre isso.

03

Veja que de qualquer lugar, elas são iguais, as de cima e as de baixo da linha, entre elas. Outra coisa legal é que da para fazer operações com elas, como vimos aqui:

04

Mas com elas estamos conseguindo entender uma tendencia no espaço, sendo que esse espaço é formado por duas dimensões, dois vetores independentes que são as abundancias das duas especies, essas figuras são desses post aqui, se alguém quiser olhar.

Como eu so consigo entender bem as coisas fazendo gráficos, eu fui tentar ver um jeito legal de representar isso no R, e um pacote legal para fazer isso é o scatterplot3d.

Esse pacote, a primeira vez que você olha, ele so tem uma função, o que faz ele parecer simples, mas o autor fez milagre com ele, quando você ve ele com calma.

Bem ele basicamente faz gráficos em três dimensões assim:

01

Para termos três dimensões, precisamos de três vetores linearmente independentes, porque se não teríamos um plano ou menos dimensões. Mas o que diabos quer dizer linearmente independente? Para esse caso, de três dimensões, quer dizer que com operações em três vetores, podemos construir qualquer um quarto vetor, dentro desse espaço.

Salvo os tropeços que devo estar dando aqui, vamos pensar que dois vetores LI (linearmente independentes) fazem um plano, se temos ai por exemplo x, y e z, os vetores x e y fazem um plano, os vetores x e z fazem outro plano, e os vetores y e z outro plano ainda distinto.

02

Legal né, essas três paredinhas. Agora veja que ai dentro a gente pode colocar uma setinha qualquer…

03

Veja que ela esta iniciando no ponto x=0, y=0 e z=0, e terminando no ponto x=1, y=1 e z=1, a pontinha da flecha.

É meio difícil enxergar, mas a gente pode subir o plano formando por x e y para facilitar.

04

Agora o que quer dizer, falar que podemos decompor ele em três vetores linearmente independentes?
Veja que esse cara ai pode ser construído assim:

Podemos começar somente pensando no plano formado pelos vetores x e y, e somar os dois vetores la, os vetores x e y. Bem eu estou pecando um pouco aqui em não ficar definindo as formulas, nem mostrando os vetores, mas a ideia é so visualizar o que esta acontecendo. E a soma é pela regra do triangulo. Eu não vou dar conta de explicar aqui agora, mas da para pegar o “feeling”, é só a gente começar um vetor no final do outro e fechar o triângulo.

05

Somando os vetores x e y, azuizinhos, temos o vetor resultante vermelho. Mas até aqui estamos olhando somente o plano formado por x e y, e agora?
Bem, podemos ver que esse vetor resultante da soma de x e y acaba em um ponto, agora podemos passar uma reta ali, furando o plano xy.

06

E essa reta ai é paralela a um outro vetor, o vetor z.

08

Lembra que o vetor, pode ser representado em começando de qualquer ponto, o vetor é o conjunto de todas as setinhas possíveis no espaço, mas então podemos representar z ali naquela reta, começando naquele ponto que é o fim de vetor somado de x e y.

09

Agora, o primeiro vetor que a gente olhou, é na verdade o x+y somado ao vetor z, ou seja descrevemos um vetor pela soma de x+y+z. E podemos construir qualquer vetor nesse espaço desse jeito, somando esses três carinhas. Ou seja com esses três vetores de base, fazemos qualquer vetor nesse espaço.

10

Bem agora da para ficar brincando de apontar vetores em 3d no R.

Eu fiz essas 2 funções para facilitar, e ter que digitar menos.

vga.vetor<-function(vetor,origem=c(0,0,0),...) {
    arrows(x0=R3$xyz.convert(x=origem[1],y=origem[2],z=origem[3])$x,
           y0=R3$xyz.convert(x=origem[1],y=origem[2],z=origem[3])$y,
           x1=R3$xyz.convert(x=origem[1]+vetor[1],y=origem[2]+vetor[2],origem[3]+vetor[3])$x,
           y1=R3$xyz.convert(x=origem[1]+vetor[1],y=origem[2]+vetor[2],origem[3]+vetor[3])$y,
           ...)
    }
 
vga.reta<-function(vetor,origem=c(0,0,0),...) {
    xyz1 <- R3$xyz.convert(origem[1],origem[2],origem[3])
    xyz2 <- R3$xyz.convert(origem[1]+vetor[1],origem[2]+vetor[2],origem[3]+vetor[3])
    segments(xyz1$x, xyz1$y, xyz2$x, xyz2$y,...)
    }

Acho que quase todo curso de engenharia tem uma disciplina chamada vga, vetores e geometria analitica, então é legal ter uma ferramente para ficar desenhando vetores.

Nesse pacote, quando você faz um grafico, ao inves de simplismentes digitar

scatterplot3d(0,0,0)

Para plotar o ponto x=0,y=0 e z=0, salve em um objeto o plot…

R3<-scatterplot3d(0,0,0)

E para esse caso, ficaremos com varias funções, chamadas a partir de R3$função para fazer coisas no gráfico, além de converter pontos ou coordenadas para usar funções básicas do R, como ali onde usamos arrows para flechinhas e segments para plotar retas. É realmente muito engenhoso o jeito que isso é feito, e acho que tem que entender muito de álgebra linear, mas no fim pode ajudar bastante a estudar, e entender um pouco melhor esses conceitos complicados.

Bem é isso, o script esta aqui no fim do post e no repositório recologia do github. Um abraço e até o próximo post.

#install.packages("scatterplot3d")
library(scatterplot3d)
?scatterplot3d
x=0
y=0
z=0
 
R3<-scatterplot3d(x,y,z,xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),zlim=c(-1.5,1.5),type="n",main="",box=T)
 
 
#Plano YZ
x0 <- -1.5
xyz1 <- R3$xyz.convert(rep(x0, 7), rep( -1.5, 7), seq(-1.5, 1.5, by=0.5))
xyz2 <- R3$xyz.convert(rep(x0, 7), rep(1.5, 7), seq(-1.5, 1.5, by=0.5))
segments(xyz1$x, xyz1$y, xyz2$x, xyz2$y, lty="dotted",col="gray")
 
xyz1 <- R3$xyz.convert(rep(x0, 7), seq(-1.5, 1.5, by=0.5), rep(-1.5, 7))
xyz2 <- R3$xyz.convert(rep(x0, 7), seq(-1.5, 1.5, by=0.5), rep(1.5, 7))
segments(xyz1$x, xyz1$y, xyz2$x, xyz2$y, lty="dotted",col="gray")
 
 
#Plano XZ
y0 <- 1.5
xyz1 <- R3$xyz.convert(rep(-1.5, 7),  rep(y0, 7),seq(-1.5, 1.5, by=0.5))
xyz2 <- R3$xyz.convert(rep( 1.5, 7),  rep(y0, 7),seq(-1.5, 1.5, by=0.5))
segments(xyz1$x, xyz1$y, xyz2$x, xyz2$y, lty="dotted",col="gray")
 
xyz1 <- R3$xyz.convert(seq(-1.5, 1.5, by=0.5), rep(y0, 7), rep(-1.5, 7))
xyz2 <- R3$xyz.convert(seq(-1.5, 1.5, by=0.5), rep(y0, 7), rep( 1.5, 7))
segments(xyz1$x, xyz1$y, xyz2$x, xyz2$y, lty="dotted",col="gray")
 
arrows(x0=R3$xyz.convert(x=0,y=0,z=0)$x,
       y0=R3$xyz.convert(x=0,y=0,z=0)$y,
       x1=R3$xyz.convert(x=1,y=1,z=1)$x,
       y1=R3$xyz.convert(x=1,y=1,z=1)$y,
       length=0.1,lwd=2,col="red")
 
R3$plane3d(Intercept=c(0,0,0),lty=2,col="blue",lwd=0.5)
 
vga.vetor<-function(vetor,origem=c(0,0,0),...) {
    arrows(x0=R3$xyz.convert(x=origem[1],y=origem[2],z=origem[3])$x,
           y0=R3$xyz.convert(x=origem[1],y=origem[2],z=origem[3])$y,
           x1=R3$xyz.convert(x=origem[1]+vetor[1],y=origem[2]+vetor[2],origem[3]+vetor[3])$x,
           y1=R3$xyz.convert(x=origem[1]+vetor[1],y=origem[2]+vetor[2],origem[3]+vetor[3])$y,
           ...)
    }
 
vga.reta<-function(vetor,origem=c(0,0,0),...) {
    xyz1 <- R3$xyz.convert(origem[1],origem[2],origem[3])
    xyz2 <- R3$xyz.convert(origem[1]+vetor[1],origem[2]+vetor[2],origem[3]+vetor[3])
    segments(xyz1$x, xyz1$y, xyz2$x, xyz2$y,...)
    }
 
R3<-scatterplot3d(0,0,0,xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),zlim=c(-1.5,1.5),type="n",main="")
R3$plane3d(Intercept=c(0,0,0),col="gray",lty=3)
 
vga.vetor(c(1,1,1),length=0.1)
vga.vetor(c(1,0,0),length=0.1,col="blue",lty=3)
vga.vetor(c(0,1,0),length=0.1,col="blue",lty=3)
vga.vetor(c(0,0,1),length=0.1,col="blue",lty=3)
 
R3<-scatterplot3d(x,y,z,xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),zlim=c(-1.5,1.5),type="n",main="")
R3$plane3d(Intercept=c(0,0,0),col="gray",lty=3)
vga.vetor(c(1,0,0),length=0.1,col="blue",lty=3)
vga.vetor(c(0,1,0),c(1,0,0),length=0.1,col="blue",lty=3)
vga.vetor(c(1,1,0),length=0.1,col="red",lty=2)
 
R3<-scatterplot3d(x,y,z,xlim=c(-2,2),ylim=c(-2,2),zlim=c(-2,2),type="n",main="")
R3$plane3d(Intercept=c(0,0,0),col="gray",lty=3)
vga.reta(c(0,0,4),c(1,1,-2),col="gray",lty=3)
vga.vetor(c(1,1,1),length=0.1)
vga.vetor(c(1,0,0),length=0.1,col="blue",lty=3)
vga.vetor(c(0,1,0),length=0.1,col="blue",lty=3)
vga.vetor(c(0,0,1),c(1,1,0),length=0.1,col="blue",lty=3)

Machine Learning – Classificação binária bayesiana

Recentemente eu comecei a ler um livro chamado Machine Learning for Hackers. Ele é bem legal, com o intuito de ensinar por exemplos.
Bem, segundo o wikipedia, Machine learning, é um ramo da inteligencia artificial que estuda sistemas que podem aprender a partir de dados, por exemplo um sistema que pode ser treinado a distinguir mensagens de e-mail, se estas são spam ou não. Após o treinamento, o sistema pode então ser usado para classificar novas mensagens se são spam ou não. Esse treinamente é algo como delimitar as fronteiras de o que é uma classe ou outra:

01

Bem legal né. Então eu decidi reproduzir o primeiro exemplo do livro, caso alguém queira dar uma olhada, todos os exemplos do livros estão num repositório de um dos autores aqui.

O primeiro exemplo explora os dados de mensagem de um projeto chamado SpamAssassin, um filtro de spam opensource, que da pra dar uma olhada nele nesse endereço aqui, mas no repositório do github, esta tudo pronto, os e-mails que serão usados, são diretórios nos quais cada e-mail esta num arquivo txt.

Esses são arquivos de texto, então pra transformar eles em números, a gente utiliza o pacote tm para text mining,ou mineração de texto.

Bem a idéia é criar um modelo com fronteiras não lineares entre o que é spam ou não.
Por exemplo a ocorrência da palavra html é um sinal que a mensagem é um spam, normalmente spam vai ter um link pra você clicar e sei la, instalar um vírus, ver uma propaganda indesejada, abrir alguma janela. Mas não existe uma relação simples, do tipo quanto mais vezes a palavra html aparecer, mais a mensagem é um spam. Existe muitas outras palavras na mensagem que podem diminuir nossa crença de que a mensagem é um spam.

Mas a solução é probabilidade condicional, que a gente ja viu aqui. Ou melhor se a gente vai apostar na soma do resultado de dois dados de 6 lados, e ja saiu um 6 no primeiro dado, ninguém vai querer apostar que a soma vai dar 1, 3 ou 5, ou qualquer valor menor que seis correto?

Para colocar essa idéia em prática, a gente usa o algorítimo conhecido como Naive Bayes classifier. De forma geral ele funciona assim, uma fruta é considerada uma maça se ela é vermelha, arredondada e tem uma tamanho x. A Naive Bayes classifier considera esses três quesitos independentemente para a probabilidade da fruta ser ou não uma maça, independentemente da outra. Ou seja uma maça verde ainda pode ser classificada como maça, dependendo do seu tamanho e forma, e algumas irregularidades na forma são aceitáveis. Mas ele não usa exatamente o teorema de bayes, ja que apesar de multiplicar por um prior, ele não divide por uma probabilidade, por isso a probabilidade vai ficando um número pequeno, já que só temos multiplicações.
Tudo que a gente precisa é um set inicial de dados já classificados para treinar, estabelecer limites para esses parâmetros e depois podemos usar o modelo resultante em novos dados.

No caso de e-mail, algumas palavras são batatas num e-mail de spam. Por exemplo palavras como html e table.

Lembre-se que quando você le nesse post algo assim:

Naive Bayes classifier.

Na verdade o texto esta em html, o computador esta lendo:

< a href="http://en.wikipedia.org/wiki/Naive_Bayes_classifier">Naive Bayes classifier.

Nos dados de teste por exemplo, a ocorrencia dessas palavras é mais o menos assim:

Term | Spam | Ham html | 377 | 9 table | 1,182 | 43

Então podemos usar a ocorrência dessas palavras para delimitar a fronteira do que é spam ou não.

01_init_plot1

Só que mesmo mensagens que não são spam, vira e mexe teremos o envio de links, então podemos usar a ocorrência de outras palavras para corrigir um possível erro de chamar de spam uma mensagem legítima, ou o contrario, que é chamar uma mensagem legítima de spam.

Bem a parte de mineração do texto, eu acho bem complicada, além de o pacote tm ter uma estrutura de dados própria para melhor manipular essas informações, isso precisa de uma boa leitura, então pulamos para depois.

Mas se você pegar todas as mensagens, pode contar a ocorrência de certas palavras, bem com a chance de elas ocorrerem. Fazendo isso para os e-mails classificados como spam temos…

head(spam.df[with(spam.df, order(-occurrence)),]) term frequency density occurrence 2122 html 377 0.005665595 0.338 538 body 324 0.004869105 0.298 4313 table 1182 0.017763217 0.284 1435 email 661 0.009933576 0.262 1736 font 867 0.013029365 0.262 1942 head 254 0.003817138 0.246

Olha ai, em spam, quase todos os termos tem algo a ver com html. Ja para as mensagens que não são spam…

head(easyham.df[with(easyham.df, order(-occurrence)),]) term frequency density occurrence 3553 yahoo 185 0.008712853 0.180 966 dont 141 0.006640607 0.090 2343 people 183 0.008618660 0.086 1871 linux 159 0.007488344 0.084 1876 list 103 0.004850940 0.078 3240 time 91 0.004285782 0.064

A maior ocorrência, a palavra yahoo, não passa de 18% das mensagens e as frequências de ocorrências são muito baixas. Ou seja, ai em cima a gente vê a informação que permite a classificação de mensagens.

Agora que sabemos como são mensagens de spam e mensagens que não são spam, podemos implementar o classificador. Isso é mais simples que parece, se a frequência de ver html em uma mensagem é 0.30 e a chance de ver a palavra table na mensagem é 0.10, então a chance de ver os dois é 0.30 × 0.10 = 0.03. Agora um problema é quando a gente ve uma palavra que a gente nunca viu quando treinamos o classificador, porque se a gente simplesmente multiplicar por 0, 0% de chance, o resultado final do produtorio vai ser zero, ai vai ser ruim, então aqui, quando uma palavra é nova, é atribuído uma probabilidade bem baixinha do tipo 0.00001% de ocorrência dela tanto para spam como não spam, mas existem maneira mais elaboradas de lidar com esse problema.
Outra coisa é que usamos um prior, para esse partimos do principio que existe 50% de chance de cada mensagem ser spam ou não spam, meio a meio, mas isso não precisa ser necessariamente assim.

Agora se a gente, pegar um novo set de mensagens, a gente faz o processo inverso, se antes nos usamos mensagens de spam e não spam para delimitar as fronteiras do que é ou não um spam, agora é so usar essas fronteiras definidas para dar nome a uma mensagem, se ela é ou não um spam, ou mais precisamente, dar uma chance dela ser ou não um spam, a partir daqui temos que dar um limite de o que será um spam, podemos cometer, como na estatística que vemos comumente, tanto erros do tipo 1 como erros do tipo 2, mas aqui estamos comparando se a chance de ser spam é maior que a chance de ser não spam.

Nos três tipos de dados, os spam, o easyham que são mensagens fáceis de identificar como não spam, e o hardham, que são mais difíceis de identificar, por conter mais termos como html e table por exemplo e ainda sim não serem spam.

03_final_classification

Veja que a linha no meio do gráfico representa a linha de 1 pra 1, quando a probabilidade de ser spam ou não spam são iguais, então pra cima temos uma chance maior de uma mensagem ser spam e para baixo da linha a chance de não ser spam é maior que a chance de ser spam.

O classificador funciona relativamente bem, salvo os erros que passam, mas nada é perfeito, claro que num servidor de e-mail sério, não teremos algo assim tão simples, eu espero, e coisas como lista negra para e-mail ja deve dar uma boa limpada também no que é e-mail de spam. Mas visto isso, machine-learning parece ser algo bem legal, e que pode ser muito útil.

Bem, além de ler o livro, outro lugar legal para dar uma olhada é esse curso do coursera aqui sobre Machine Learning, do Andrew Ng.
Eu me inscrevi nele, e espero começar a postar alguma coisa sobre o tópico a partir de agora. Afinal de contas, curiosidade não mata ninguém. E olhando novamente esse post, ele está mais com cara de agenda que sei la o que.

Por enquanto não temos código, mas a tarefa de casa é usar esse classificador com os dados de exemplo que vem no R chamado iris :), tópico para um próximo post com certeza.

Problema – Rosalind – Counting Phylogenetic Ancestors – INOD

Esse problema é bem simples, a pergunta consiste em, dado um número n inteiro, que é o número de espécies que vamos avaliar, quantos ancestrais eles tem, quantos ancestrais precisamos para fazer uma árvore filogenética ligando eles.

A resposta aqui é bem simples, e nem precisa de um programa na verdade. A resposta é simplesmente n-2.

Mas para chegar até essa resposta, a gente tem que olhar algo legal, as árvores binárias sem raiz.

Essa é a estrutura de dados que é gerada quando a gente vai fazer a árvore genealógica de uma família por exemplo. Que é um tipo especial de árvore binária. Poxa, como teorias dos grafos é uma coisa legal, ja esbarramos neles aqui e aqui, mas essas eram estruturas e ideias bem diferentes da daqui, usando grafos bipartidos.

781px-Waldburg_Ahnentafel

Pela imagem , que é da família de Sigmund Christoph von Waldburg-Zeil-Trauchburg de 1776, da para ter a ideia, que mesmo a gente pensando isso hoje, essa estrutura de dados já esta por ai a muito tempo.

02

Em teoria dos grafos, esse grafo, que é uma estrutura de dados de árvore, é caracterizado pelo fato de cada vértice, ou node, ou bolinha se preferir, ter no máximo dois “filhos”, então ele vai ter um grau 3 no máximo, ou seja 3 arestas ou edge que saem da uma bolinha, no máximo, poderiam ser 2 como numa folha porém. Veja que cada bolinha tem no máximo 3 conexões. a gente fez um problema ja comentando sobre isso aqui.

Agora olhando o texto do wikipedia sobre árvores binárias sem raiz, a própria definição diz que para n folhas da árvore, que na biologia são as espécies, haverão n-2 vértices internos, ou ancestrais, e somente assim é possível construir essas árvores. Por isso o exemplo diz que para 4 espécies, precisamos de 2 ancestrais, e assim vai, n-2 vértices internos na árvore.

Agora outra coisa legal, que pelo menos eu sempre tive dificuldade de entender ou pensar, é como as pessoas veem coisas diferentes na mesma figura, comparando por exemplo computação e biologia.

Por exemplo, numa aula de grafos, uma árvore binaria sem raiz vai ser apresentada assim, talvez, num sei, eu imagino, apesar que nunca tive aulas sobres grafos:

01

Numa filogenia, a terminação ali é uma espécie, árvores binárias sem raiz não tem ciclos, ou seja, todas as pontas terminam numa extremidade solta, e além disse se você pensa na união dos vértices ali, todo mundo ta ligado a um ancestral.

O problema é que a primeira vista, isso não se parasse nada com o que vimos nesse post aqui não?

03

A gente quer ver uma filogenia mais desse jeito, uma árvore desse jeito.

Mas a partir daquela figura:

01

A gente tem que pensar como se pegássemos uma parte…

02

Ai a gente puxa essa parte…

03

Dai se apontarmos todas as espécies para o mesmo lado…

04

Fica com uma cara muito melhor de filogenia, ou árvore filogenética não? Da para dar mais uma melhoradinha ainda…

05

E olha ai, agora parece que da para entender, ver que é tudo igual, so depende da perspectiva. Depois de entender que o jeito que a gente visualiza as coisas é relativo, fica mais fácil entender da onde diabos vem o maximum likelihood para desenhar árvores, ja que a gente sempre pensa em mínimos quadrados da perspectiva de uma anova, como aqui onde somamos as distancias da média, mas agora que temos folhas, que são espécies ligadas por risquinhos, se a gente pensar que essa distância entre cada espécie, o caminho de uma folha a outra é o quão igual elas devem ser, da para começar a imaginar como se coloca o conceito de verossimilhança nesse problema. E além disso, se da para usar verissimilhança, da para ter um “pre-conceito”, uma noção de como a árvore deve ser e usar isso de prior para as distâncias das espécies, numa perspectiva bayesiana, mas isso fica para um próximo post.

O algoritimo de Metropolis–Hastings aplicado a distribuição normal

Bem a gente ja viu aqui como é um algorítimo de Metropolis–Hastings usado para a distribuição binomial. Mas para muitas coisas a gente usa a distribuição normal como ele fica para a distribuição normal?

Eu peguei aqui um exemplo da página do Brian Neelon aqui. Tem bastante coisa legal para olhar, mas vamos la.

A primeira parte interessante, é que como já vimos nos modelos usando openbugs que fazemos, é que o desvio-padrão, a gente estima indiretamente, usando o tau que é 1/variância. Aqui tem uma explicação melhor e mais detalhada, mas é simples pensar que a média a gente pode usar uma distribuição normal para representar, agora o desvio-padrão não, já que ele não pode assumir valores negativos, assim como uma distribuição de poisson, mas pode assumir qualquer valor entre 0 e + infinito, dai precisamos de algo como a distribuição gamma que é o conjugate prior do desvio padrão da distribuição normal, nesse post aqui tem algo sobre conjugate prior se alguém quiser olhar.

Mas então a gente ta fazendo isso aqui.

normal

Beleza. A idéia é que o que falei ali em cima é o que esta nessa figura, eu espero, se não estiver entendendo nada errado.

Então precisamos do nosso Prior

###########
# Priors  #
###########
mu.prior<-0
sigma2.prior<-10
tau.prior<-1/sigma2.prior
a<-0.001
mtau<-a+n/2

Além disso começamos a cadeia em qualquer lugar, eu comecei aqui em média igual a 0.1 e tau igual a 0.1 também, só para ver aquela caminhadinha inicial, mas normalmente queremos começar com um valor aleatório aqui. Ja que se de qualquer lugar que começamos, terminamos no mesmo lugar, é muito mais convincente que se delimitamos sempre os lugares onde começamos.

####################
# Inicio da cadeia #
####################
 
mu<-0.1
tau<-0.1
nsim<-5000
 
##########
# Cadeia #
##########
cadeia.mu<-c(mu,rep(0,nsim-1))
cadeia.tau<-c(tau,rep(0,nsim-1))
A<-0

Agora é so rodar o MCMC. Mas antes um pequeno comentário. Os valores de máxima verosimilhança ficam bem pequenos, apesar de que para a distribuição binomial ficou algo com 4 casas depois do zero, aqui o likelihood vai ser um produtório de vários números pequenininhos, que vai resultar num número incrivelmente pequeno.

Então gerando vários valores possíveis para média e calculando o likelihood vemos o seguinte.

01

Agora olhe os valores no eixo y, 3.502837e-75 no topo la, é um número muito muito pequeno. Números grandões gastam assim (grande no sentido de número de casas) gastam muito bytes de memória para serem guardados, inclusive temos um limite que podemos usar, quando tentamos usar mais que o limite temos um estouro da memória ou overflow.

Agora o que a gente vê normalmente é as pessoas trabalhando com loglikelihood, que é meramente calcular o log natural do valor usado ali em cima. E quando a gente faz isso, veja o que acontece com o mesmo gráfico.

02

Os valores são bem menores, já que a gente muda do produtório para o um somatório. Agora séra que isso da problema na estimativa de maximum likelihood? Que nada mais é que o pico desses gráficos ae!

Veja só, a gente fez um vetor chamado valores, e cada elemento desse vetor teve o likelihood calculado, o maior valor que encontramos foi:

> max(likelihood) [1] 3.502837e-75

Agora a gente pode olhar qual elemento do vetor é isso.

> which(likelihood==max(likelihood)) [1] 101

Agora se a gente olhar qual valor representa esse elemento do vetor de valores vemos que:

> valores[101] [1] 5

E olha a média dos dados são, a média é exatamente a estimativa de maxima verosimilhança;

> mean(dados) [1] 5.038861

E se olhamos o proximo, o elemento 102, o elemento 101 é muito mais próximo da média.

> valores[102] [1] 5.1

Agora o valor em log é um valor diferente, ja que fizemos o log do valor acima:

> max(loglikelihood) [1] -171.4403

Mas o pico continua sendo no mesmo lugar:

> which(loglikelihood==max(loglikelihood)) [1] 101

Agora algo legal é como, devido aos números serem pequenos, não da para calcular exatamente, ou seja, o contrario de log, que é o exponencial do loglikelihood não é igual o likelihood.

> exp(loglikelihood)[100]==likelihood[100] [1] FALSE

Mas isso é por causa de aproximação, não tem como representar no computador um número como pi por exemplo, ja que não da para ter infinitas casas, mas para fazer comparações assim, sem sofrer com erros de aproximação, comum quando entramos nesse tipo de conta, da pra usar a função all.equal do R.

> all.equal(exp(loglikelihood)[100],likelihood[100]) [1] TRUE

Mais comentários sobre isso estão no faq do R.

Mas vamos a como fazer o MCMC usando Metropolis-Hastings para uma distribuição normal.

############
#   MCMC   #
############
for (t in 2:nsim){
    mus<-mu+0.3*rt(1,3)
 
    r<- (sum(dnorm(dados,mus,1/sqrt(tau),log=T) - dnorm(dados,mu,1/sqrt(tau),log=T)) +
         dnorm(mus,mu.prior,sqrt(sigma2.prior),log=T) -  dnorm(mu,mu.prior,sqrt(sigma2.prior),log=T))
 
    if (log(runif(1))<r) {
        mu<-mus
        A<-A+1
        }
 
    cadeia.mu[t]<-mu
 
    tau<-rgamma(1,mtau,b+sum((dados-mu)^2)/2)
    cadeia.tau[t]<-tau
}

Bem curtinho não. Mas seguimos o mesmo esquema, a cada iteração, sugerimos uma nova média para a cadeia, aqui usando uma distribuição t, depois disso calculamos uma chance de aceitar, que aqui parece um pouco mais complicado, mas note que estamos usando o log em tudo, log=T, então os produtos viram soma, divisão subtração, depois vemos isso com calma, mas é a mesma coisa feita para a distribuição binomial, o likelihood * prior dividido pelo do valor anterior da cadeia, mas a gente ve a chance de aceitar o novo valor ou não, dai sorteamos se vamos aceitar ele ou não, se sim adicionamos ele a cadeia, se não repetimos o ultimo valor, depois disso usamos a distribuição gamma para ajustar a precisão tau com a nova ou antiga média. Essa parte precisa de um pouco mais de estudo mas vamos indo.

Na verdade ja acabou, temos nossa cadeia final, ou cadeia, uma para média e uma para o desvio padrão. Podemos olhar os nossos resultados.

03

E as cadeia convergem bonitinho, em baixo nos histogramas eu tentei olhar varios pedacinhos da cadeia para a média e o desvio padrão, e vemos como ambas vão para intervalos que cobrem os valores reais usados para gerar os dados.

Podemos também olhar como a cadeia caminha em relação aos dois parametros, média e desvio padrão.

04

Novamente vemos que existem muitas nuances, que espero fazer um post experimentando o que acontece de acordo com a estrategia usada para gerar novos valores de parâmetros propostos, mas ja da para começar a ter uma ideia de como é um MCMC para uma distribuição normal que tem dois parâmetros, diferente do que fizemos para a distruição binomial que so tinha um parâmetro.

Bem é isso, o repositório recologia no github ta atualizado se alguém quiser alguma coisa, e ficamos por aqui. Quem quiser pode olhar o script original desse MCMC na página do Brian Neelon.

#Original
#http://people.duke.edu/~neelo003/r/
set.seed(250)
 
#######
#Dados#
#######
 
n<-100
mu.real<-5
sigma2.real<-2
tau.real<-1/sigma2.real
dados<-rnorm(n,mu.real,sqrt(sigma2.real))
 
a<-0.001
b<-0.001
 
###########
# Priors  #
###########
mu.prior<-0
sigma2.prior<-10
tau.prior<-1/sigma2.prior
mtau<-a+n/2
 
####################
# Inicio da cadeia #
####################
 
mu<-0.1
tau<-0.1
nsim<-5000
 
##########
# Cadeia #
##########
cadeia.mu<-c(mu,rep(0,nsim-1))
cadeia.tau<-c(tau,rep(0,nsim-1))
A<-0
 
#Valores para avaliar a maxima verossimilhança
valores<-seq(-5,15,by=0.1)
likelihood<-vector()
 
for(i in 1:length(valores)) {
    likelihood[i]<-prod(dnorm(dados,valores[i],sqrt(2),log=F))
    }
 
#likelihood
plot(likelihood~valores,type="b",frame=F,ylab="likelihood")
 
loglikelihood<-vector()
 
for(i in 1:length(valores)) {
    loglikelihood[i]<-sum(dnorm(dados,valores[i],sqrt(2),log=T))
    }
 
#Loglikelihood
plot(loglikelihood~valores,type="b",frame=F,ylab="loglikelihood")
 
exp(loglikelihood)[100]==likelihood[100]
all.equal(exp(loglikelihood)[100],likelihood[100])
 
max(likelihood)
 
which(likelihood==max(likelihood))
valores[101]
mean(dados)
valores[102]
 
max(loglikelihood)
which(loglikelihood==max(loglikelihood))
 
 
############
#   MCMC   #
############
for (t in 2:nsim){
    mus<-mu+0.3*rt(1,3)
 
    r<- (sum(dnorm(dados,mus,1/sqrt(tau),log=T) - dnorm(dados,mu,1/sqrt(tau),log=T)) +
         dnorm(mus,mu.prior,sqrt(sigma2.prior),log=T) -  dnorm(mu,mu.prior,sqrt(sigma2.prior),log=T))
 
    if (log(runif(1))<r) {
        mu<-mus
        A<-A+1
        }
 
    cadeia.mu[t]<-mu
 
    tau<-rgamma(1,mtau,b+sum((dados-mu)^2)/2)
    cadeia.tau[t]<-tau
}
 
###########
# Gráfico #
###########
 
# Cadeias
layout(matrix(c(1,1,1,2,2,2,1,1,1,2,2,2,3,4,5,6,7,8),nrow=3,byrow=T))
 
plot(seq(1,nsim) ,cadeia.mu[1:nsim],type="l",xlab="Iteration", ylab="mu",frame=F,
main="Média",col=2,lty=1,ylim=c(0,10))
abline(h=mu.real)
plot(seq(1,nsim) ,cadeia.tau[1:nsim],type="l",xlab="Iteration", ylab="mu",frame=F,
main="Tau",col=2,lty=1,ylim=c(0,1))
abline(h=1/sigma2.real)
 
par(mar=c(2,2,2,2),font=3,cex=0.7)
 
hist(cadeia.mu[1:100],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(3,6),breaks=seq(0,6,by=0.02))
hist(cadeia.mu[2000:2500],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(3,6),breaks=seq(3,6,by=0.02))
hist(cadeia.mu[4500:5000],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(3,6),breaks=seq(3,6,by=0.02))
hist(cadeia.tau[1:500],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(0,1),breaks=seq(0,1,by=0.01))
hist(cadeia.tau[2000:2500],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(0,1),breaks=seq(0,1,by=0.01))
hist(cadeia.tau[4500:5000],xlab="",ylab="Frequência",main="",cex=0.7,xlim=c(0,1),breaks=seq(0,1,by=0.01))
 
#Evolução da cadeia
plot(cadeia.mu[1:500],sqrt(1/cadeia.tau[1:500]),type="b",frame=F,xlab="Média",ylab="Desvio Padrão",cex=0.5,pch=19,
     xlim=c(0,6),ylim=c(0,6))
library(MASS)
contornos<-kde2d(cadeia.mu[1:500],sqrt(1/cadeia.tau[1:500]))
contour(contornos,add=T,lwd=2,nlevels = 6,col=rev(heat.colors(6)))