Hidden Markov Chain

Rapaz, se uma cadeia de Markov já era uma ideia difícil de digerir, imagina ela escondida?
Eu me perguntei por um tempo, o que diabos era isso?

Mas uma Cadeia de Markov Oculta, aka Hidden Markov Chain está mais próxima da biologia que imaginamos.

Bem, vamos continuar o exemplo anterior, quanto tratamos de Cadeias de Markov num post anterior.

Então temos a matrix de transição, o clima vai de “estados”, entre sol e chuva. Mas e se eu não puder saber o clima? Por exemplo, esse clima é o da cidade da minha amiga, e eu não tenho como saber se esta um dia de sol ou um dia de chuva la na cidade dela.

O processo contínua acontecendo, independente se eu posso ver ele ou não, se eu posso coletar amostras deles ou não. Mas nesse caso eu posso coletar amostras de algo relacionado a ele. Amostras do que minha amiga esta fazendo na cidade.

Eu sei que minha amiga faz quatro coisas. Ela deixa no status do msn dela o que ela esta fazendo.
As atividades dela são:

[1] "Dormir" "Limpar" "Shopping" "Caminhar"

E eu sei que o que ela faz depende se o dia é de sol ou de chuva, da seguinte forma:

Dormir Limpar Shopping Caminhar Sol 0.05 0.15 0.40 0.40 Chuva 0.80 0.10 0.05 0.05

Dia de sol, ela preferencialmente (maior probabilidade) ela caminha ou vai ao shooping, mas ela também limpa a casa e dorme se estiver com sono.
Agora dia de chuva, hummmm chuvinha é bom para dormir, mas as vezes não tem jeito e ela tem que limpar a casa, ir ao shooping ou caminha na chuva com seu guardachuvinha.

Então esta rolando um processo Markoviano aqui, que é o clima, e ele esta refletindo no comportamento da minha amiga, mas eu so vejo o que minha amiga faz, não o clima da cidade dela.

Vamos simular alguns dados e olhar na forma de uma figura:

01

Aqui temos 50 dias, na parte superior temos como no post anterior uma simulação do clima, na parte inferior temos o comportamento da minha amiga.
Olha so, normalmente coincide de o dia de chuva ela dormir, mas olhem que la pelo dia 45 ela dormiu num dia de sol, ela devia estar muito cansada coitada. Mas durante o sol vemos que ela oscila mais comumente entre compras e shooping.
Mas de modo geral ha algum reflexo do que clima no que ela esta fazendo, que é a emissão vista naquelas probabilidades que vimos la em cima.

Mas vamos olhar uma quantidade de dias um pouco maior:

02

Então olhando 500 dias, parece que continuamos vendo um reflexo do comportamento da minha amiga em relação ao clima. Mas vimos que uma propriedade da cadeia de Markov, assegurando algumas premissas, que podemos prever a porcentagem de cada estado ocorrer, proporção entre chuva e sol. Poxa se o resultado é um processo markoviano que conhecemos o comportamento, sera que se olharmos so para o comportamento da minha amiga podemos prever o tempo?

Então tudo que sabemos é o comportamento dela ao longo dos dias(primeiros seis casos abaixo):

head(amostras) [1] "Caminhar" "Caminhar" "Caminhar" "Caminhar" "Dormir" "Dormir"

E as como é o processo do clima na cidade e o comportamento da minha amiga:

> prob_trans Sol Chuva Sol 0.9 0.1 Chuva 0.5 0.5 > prob_emiss Dormir Limpar Shopping Caminhar Sol 0.05 0.15 0.40 0.40 Chuva 0.80 0.10 0.05 0.05

Um cara chamado Andrew Viterbi falou que isso é informação suficiente.
A idéia é que pelo comportamento podemos chutar o estado do dia, sabendo o estado do dia podemos chutar o dia posterior e anterior, e comparar com o comportamento dos dias anterior e posterior, e assim chutar e chutar, e buscar a máxima verosimilhança da configuração de dias mais provável de fazer esses comportamentos acontecerem.

E seguindo a proposta do Viterbi temos o seguinte para os primeiros dias visto acima:

head(ifelse(apply(v, 1, function(x) which.max(x))==1,"Sol","Chuva")) [1] "Sol" "Sol" "Sol" "Sol" "Chuva" "Chuva"

Hummm, mas melhor que isso, como viemos de uma simulação, sabemos a “realidade”, então podemos comparar os resultados e ver se deu certo essas previsões.

prev Chuva Sol Chuva 70 18 Sol 13 399

E temos que em 399 casos, era sol nos dados originais e essa foi a previsão feita, muito bom. Prevemos também 70 dias de chuva corretamente. Mas houve 13 dias que eram sol e falamos que eram dias de chuva e 18 casos que foi o contrario, era chuva e falamos que era sol.

Mas se somarmos 70 e 399, e dividir por 500, que foi o total de observações, vemos que:

sum(diag(prev))/sum(prev) [1] 0.938

Poxa, é um índice de acerto de quase 94%.

Bem eu sei tudo isso parece num ter nenhuma relação com biologia, com ecologia, mas estamos muito longe estar certos ao pensar assim, isso esta mais perto do que nunca, e não pelo fato de usarmos MCMC para estatística Bayesiana apenas. Por exemplo, eu não tenho culhões ainda para programar sozinho esse algorítimo do Viterbi, mas olha o titulo do livro da onde tirei ele :).
Mas outra coisa de explodir a cabeça é pensar, como muitas vezes desprezamos os dados que temos. Aqui é um exemplo onde observações podem ser usadas para inferir o estado do clima, sem nem mesmo observar o clima. Isso tudo é bem legal não?

#Hidden Markov Chain
set.seed(123)
#Wim P. Krijnen – Applied Statistics for Bioinformatics using R",pag.209
#cran.r-project.org/doc/contrib/Krijnen-IntroBioInfStatistics.pdf
 
#função que aplica o algoritimo de Viterbi
makeViterbimat <- function(sequence,transitionmatrix,emissionmatrix) {
sequence<-sequence
numstates<-dim(transitionmatrix)[1]
v<-matrix(NA, nrow = length(sequence), ncol = dim(transitionmatrix)[1])
v[1, ] <- 0
v[1,1] <- 1
for (i in 2:length(sequence)) {
for (l in 1:numstates) {
statelprobnucleotidei <- emissionmatrix[l,sequence[i]]
v[i,l] <- statelprobnucleotidei * max(v[(i-1),] * transitionmatrix[,l])
}
}
return(v)
}
 
#estados possiveis
estados<-c("Sol","Chuva")
#observaçoes possiveis
observacao<-c("Dormir","Limpar","Shopping","Caminhar")
#probabilidades iniciais para a simulação
prob.inicial<-c(Sol=0.6,Chuva=0.4)
 
#Probabilidades de transmissão
prob_trans<-matrix(c(0.9,0.1,0.5,0.5),ncol=2,nrow=2,byrow=T,
dimnames=list(estados,estados))
#Probabilidades de Emissão
prob_emiss<-matrix(c(0.05,0.15,0.4,0.4,0.8,0.1,0.05,0.05),ncol=4,nrow=2,byrow=T,
dimnames=list(estados,observacao))
 
#Criando objeto para receber simulação
simul<-matrix(NA,ncol=2,nrow=500)
#O primeiro dia temos que sortear
simul[1,1]<-sample(estados,1,prob=prob.inicial)
 
#simulando o clima
for(i in 2:500) {
simul[i,1]<-sample(estados,1,prob=prob_trans[simul[i-1,1],])
}
#simulando o comportamento
for(i in 1:500) {
simul[i,2]<-sample(observacao,1,prob=prob_emiss[simul[i,1],])
}
 
#Mudando os dados de observação para numeros
#para poder fazer os graficos
recod <- 1:4
obs.rec <- recod[match(simul[,2],observacao)]
 
#Figura 1
par(mfrow=c(2,1),mar=c(2,5,2,2))
plot(x=1:50,y=ifelse(simul[,1]=="Sol",1,0)[1:50],type="l",yaxt="n",
frame.plot=F,ylab="",main="Real")
axis(2,at=c(0,1),label=rev(estados),las=2)
par(mar=c(4,5,2,2))
plot(x=1:50,y=obs.rec[1:50],type="b",frame.plot=F,yaxt="n",ylab="",
xlab="Dias",main="Observado")
axis(2,at=c(1,2,3,4),label=observacao,las=2)
 
#Figura 2
par(mfrow=c(2,1),mar=c(2,5,2,2))
plot(x=1:500,y=ifelse(simul[,1]=="Sol",1,0),type="l",yaxt="n",
frame.plot=F,ylab="",main="Real")
axis(2,at=c(0,1),label=rev(estados),las=2)
par(mar=c(4,5,2,2))
plot(x=1:500,y=obs.rec,type="l",frame.plot=F,yaxt="n",ylab="",xlab="Dias"
,main="Observado")
axis(2,at=c(1,2,3,4),label=observacao,las=2)
 
#olhando os dados
amostras<-simul[,2]
head(amostras)
 
#usando a função makeViterbimat
v<-makeViterbimat(sequence=amostras,transitionmatrix=prob_trans,
emissionmatrix=prob_emiss)
#Vendo como ficaram as previsões
head(ifelse(apply(v, 1, function(x) which.max(x))==1,"Sol","Chuva"))
 
#fazendo tabela de previsões
prev<-table(ifelse(apply(v, 1, function(x) which.max(x))==1,"Sol","Chuva"),simul[,1])
prev
 
# % de Acerto
sum(diag(prev))/sum(prev)

Crescimento populacional de aves e algumas simulações.

Bem já tivemos alguns posts sobre crescimento populacional por aqui. Muita teoria, mas sempre eu ficava com dúvida quando e como esses modelos seriam confrontados com dados reais? O que adianta saber esse monte de formula se são inúteis na vida real?
Bem vamos ver um pequeno exemplo que peguei do livro “Primer of Ecology with R” do M. Henry H. Stevens. Ele tem cara de ser um cara legal hehehe.

Primeiro temos que instalar o pacote primer, que ele fez junto do livro para servir de apoio para o livro, que vem com várias funções legais, dados e outras coisas para explorarmos. Mas por agora vamos olhar o dataset chamado sparrows que vem no pacote.

01

O dados consistem da contagem de passarinhos ao longo de alguns anos. Bem como vimos o crescimento da população pode ser medido como o tamanho da população de um ano dividido pela população do ano anterior, ja que o lambda vai ser a divisão.

Vamos lembrar que podemos descrever o crescimento discreto como:

 N_{t+1}= N_t\cdot(1+r)

E esse 1+r é exatamente o

 N_{t+1}= N_t\cdot\lambda

Assim se ele é maior que 1, a população ta crescendo, se ele é menor que 1, a população esta diminuindo e caso ele seja exatamente 1, a população não mudaria de um ano para o outro.

No parte inferior do gráfico anterior estamos vendo exatamente isso, o \lambda da mudança entre os anos, mas ele parece variar bastante de um ano para o outro.

Mas se esses são os valores de \lambda, e eles se repetem, poderiamos pensar que o próximo ano, seria a população no último ano vezes um valor parecido com algum desses ae. Então podemos tentar ter uma idéia de como será a população no próximo ano se sorteamos um desses valores de crescimento e multiplicar última contagem da população.

Mas não precisamos parar por ai, podemos pegar o resultado dessa simulação para o ano seguinte e fazer denovo, sortear outro valor de crescimento e ver o que acontece, e denovo e denovo, umas 50 vezes por exemplo.

Então teríamos isso:

02

Uma previsão de como estará a população daqui a 50 anos, mas essa é apenas uma simulação dependendo da nossa sorte de ter pego somente número altos ou azar de pego muitos números baixos, veríamos a população subir ou descer, então deveríamos fazer isso denovo para ver o que vai acontecer, repetir o mesmo processo outra vez.
Alias, como as contas residem nas costas do computador, vamos fazer isso 10 vezes, e olhar um gráfico como o acima para ver o que esperar.

03

Agora vejam como aconteceu um caso, onde a população cresceu muito, muito mesmo, a linha la em cima, mas notem que ela esta la isolada, a maioria esta dentro de um mesmo patamar, mas temos também outra simulação que parece que decaiu muito, culminando na extinção da população. Então temos vários cenários possíveis, a população ficando gigante, a população extinguindo, mas em media não é isso que vemos, vemos algo entre 10 e mil indivíduos como expectativa.
Mas olhômetro tem seu limites, então o que podemos fazer é realmente muitas, mas muitas simulações mesmo e ver em média, qual o tamanho da população devera ter daqui 50 anos.

04

Então utilizamos uma função para fazer esse serviço e geramos 1000 simulações, agora é so pegar o último valor da simulação, o tamanho da população que esperamos para daqui 50 anos e olhar. Vamos fazer um histograma para entender melhor o resultado.
Bem primeiro vemos que olhar diretamente o histograma é ruim, ja que existem alguns valores gigantes, lembram do primeiro gráfico que tinha um valor de lambda de mais ou menos 3? Imagina que esse cara sai várias vezes na simulação, a população ficaria enorme, e esse são os extremos que estamos vendo. Isso faz o histograma ficar ruim de ver daquele jeito.

Mas podemos usar a escala logaritima. So uma coisa é que aqui usamos log do tamanho da população + 1, isso porque logo de 0 não existe, mas log de 1 é 0, então não temos problemas com os zeros dos dados e conseguimos ver como esta a população numa escala logaritima.

Agora vemos algo como uma distribuição normal, mas alguns valores gigantes, so que esses são poucos. Aqui podemos calcular o intervalo de confiança de 95% dos dados e colocar uma linhas ali no histograma identificando ele, o intervalo de confiança de 95%. Então esperamos que daqui 50 anos a população esteja por ai. Vemos muitos casos onde a população foi para zero, isso pode ser perigoso, isso quer dizer que a qualquer momento que a espécie baixar muito a população deveríamos preocupar, poderíamos nessa simulação ainda ver por exemplo quantos anos para frente começam a aparecer muitos casos de extinção da espécie, para ter uma idéia o intervalo mínimo que deveríamos ir olhar como a espécie esta.

Agora fizemos uma simulação, considerando o crescimento como um efeito estocástico, mas qual a confiança nessas simulações?

Bem podemos tentar olhar como é a distribuição dos crescimentos que usamos na simulação, e melhor que isso, podemos usar uma distribuição t, que leva em consideração o número de valores que temos, para ver se uma distribuição t pode descrever bem esses dados.
Então calculamos o quantiles e determinamos, usando a distribuição t pelos motivos que colocamos nesse post, o intervalo de confiança do que deve ser o lambda.

05

Mas quando fazemos um Q_Qplot vemos que temos um valor muito grande e dois valores muitos pequenos, que parecem não se enquadrar na distribuição, talvez outliers, e esses podem ter grande impacto na simulação, tanto gerando os valores extremos como aumentando em demasia os casos em que a população foi para zero, lembrando que uma vez que a população chega a zero, zero vezes qualquer número é sempre zero, então ela fica zero para o resto da simulação. Então temos que olhar com cuidado esses resultados.

Mas ainda sim temos aqui uma previsão nua do que deve acontecer nos próximos anos dessa espécie. Por um lado podemos encarar essa simulação como irreal, já que não consideramos um limite para o crescimento, como no caso do crescimento logístico, ou competição ou qualquer outra coisa. Mas essa simulação feita aqui pode ser legal exatamente por isso, estamos pensando o vai acontecer sem fazer muitas premissas, assumindo que o mundo é extremamente simples, e temos um primeiro resultado, uma primeira perspectiva para o futuro. E essa medida aqui adquirida pode ainda ser confrontada com o que se conhece da espécie, o que achamos que pode acontecer quanto a interações bióticas e como a espécie deveria estar, Algo como quando usamos o teorema de Hardy-Weinberg, onde o interessante não é o resultado em sim, mas porque o que observamos é diferente desse resultado e como é essa diferença, e isso nos traz muita informação.

Referencia: Stevens, M. H. H. A Primer of Ecology with R. UseR! Series, Springer (2009)

#####################################################################
 
set.seed(3)
library(primer)
data(sparrows)
names(sparrows)
 
#Figura 1
jpeg("01.jpg")
layout(matrix(1:2, nc=1))
plot(Count ~ Year, type="b",data=sparrows,ylab="Contagem",xlab="Ano"
,ylim=c(0,125),frame.plot=F)
obs.R <- sparrows$Count[-1]/sparrows$Count[-length(sparrows$Count)]
plot(obs.R ~ sparrows$Year[-length(sparrows$Count)],ylab="Crescimento",
xlab="Ano",ylim=c(0,3),frame.plot=F)
abline(h=1, lty=3)
dev.off()
 
#Primeira Simulação
years<-50
sim.Rs<-sample(x=obs.R,size=years,replace=TRUE)
output<-numeric(years+1)
output[1]<-sparrows$Count[sparrows$Year==max(sparrows$Year)]
 
for(t in 1:years ) {
output[t+1]<-output[t] * sim.Rs[t]
}
 
plot(0:years, output, type="l",xlab="Anos",ylab="Contagem",main="Simulação",
frame.plot=F)
 
#10 Simulações
sims<-10
sim.RM<-matrix(sample(obs.R, sims * years, replace=TRUE),
nrow=years, ncol=sims)
output[1]<-sparrows$Count[sparrows$Year==max(sparrows$Year)]
outmat<-sapply(1:sims,function(i) {
for( t in 1:years ) {
output[t+1] <- output[t] * sim.RM[t,i]
}
output
}
)
 
options(scipen=10)
matplot(0:years,outmat, type="l",xlab="Anos",ylab="Contagem",
frame.plot=F,log="y")
 
#Função de fazer Simulações
PopSim <- function(Rs, N0, years=50, sims=10) {
sim.RM<-matrix( sample(Rs,size=sims*years,replace=TRUE),nrow=years,
ncol=sims)
output<-numeric(years+1)
output[1]<-N0
outmat<-sapply(1:sims, function(i) {
for( t in 1:years ) {
output[t+1]<-round(output[t]*sim.RM[t,i],0)
}
output
})
return(outmat)
}
 
#Agora Gerando 1000 simulações
output<-PopSim(Rs=obs.R, N0=43, sims=1000)
dim(output)
 
#Expectativa da população daqui 50 anos
N.2053<-output[51,]
summary(N.2053, digits=6)
quantile(N.2053, prob=c(0.0275, .975) )
 
#Figura 4, histogramas
par(mfrow=c(2,1))
hist(N.2053, main="N",ylab="Frequencia")
hist(log10(N.2053+1), main="log(N+1)",ylab="Frequencia" )
abline(v=log10(quantile(N.2053, prob=c(0.0275, .975) )+1), lty=3)
 
#Calculando o intervalo de confiança com distribuição t
logOR<-log(obs.R)
n<-length(logOR)
t.quantiles<-qt(c(0.025, 0.975),df=n-1)
 
se<-sqrt(var(logOR)/n)
CLs95<-mean(logOR)+t.quantiles*se
 
R.limits<-exp(CLs95)
R.limits
 
N.Final.95<-sparrows$Count[sparrows$Year==max(sparrows$Year)]*R.limits^50
round(N.Final.95)
 
#Figura 5 QQplot
qqplot(qt(ppoints(n), df=n-1), scale(logOR))
qqline(scale(logOR))
 
#########################################################################

Rosalind – Independent Alleles – LIA

E la vamos nos denovo, num problema que a explicação mais complica que explica, pelo menos para mim hehe.

Então temos o Tom de genotipo AaBb que tem filhos com Maria também é AaBb, eles tem 2 filhos. Qual a chance de um filho ser do mesmo genotipo do pai AaBb?
Aqui é só lembrar de Mendel e fazendo os cruzamentos e vemos que é 0.25, vale uma olhadinha nos punnet square se ficar com dúvida desse número, mas é 50% de ser Aa vezes 50% de ser Bb que da 0.25.

Mas e nas próximas gerações? Essa é a pergunta!

Então primeiro abordamos o aspecto do crescimento, serão 2 filhos que terão 2 filhos cada gerando 4 indivíduos na próxima população, que esses quatro indivíduos tem 2 filhos cada e assim vai. Logo temos que o tamanho da população vai ser  2^k , onde 2 é o número de filhos que cada indivíduo tem e k é a geração que queremos saber o tamanho da população.

Olhando isso vemos o bom e velho crescimento exponencial.

01

Até aqui tudo bem.
Agora vamos pensar na outra parte que é, quantos caras AaBb vão estar ai.
Bem a gente tem a chance de 0.25 na primeira geração de termos um AaBb, na próxima teremos que esse 0.25 de AaBb teve filhos com outro AaBb e no final das contas a chance continua 0.25, sai fazendo cruzamento que nem um louco que deve chegar nesse número.

Então a população vai ir mudando de acordo com a distribuição binomial.

Ou seja, a chance de um evento x acontecer numa população de N indivíduos com probabilidade de 0.25. Essa frase ficou ruim mas vamos colocar num gráfico:

02

Aqui a gente ve que numa população de 5 pessoas (número que eu tirei da minha cabeça, eu não pensei em seguir o crescimento como ali em cima), temos uma chance alta de somente uma pessoa ser AaBb, agora 2 pessoas serem AaBb ja é mais baixa e assim vai.

No segundo quadro, o que a gente ve é que dificilmente somente uma pessoa vai ser AaBb, porque tem 25 pessoas na população oras, num tem como pelo menos um cara não ser AaBb, mas temos uma chance razoável de uns 5 caras serem AaBb, mas assim como 1 cara ser AaBb é quase nula ja que a população é grande, todo mundo ser AaBb também é pouco provável.

Certo, isso é a PDF (probability density function), que tem uma formula para chegar nesses números ali do lado do quadrinho que explica distribuição binomial no Wikipédia, melhor que isso é que o R tem uma função pronta para isso, mas isso nos da um insight de como será a população provavelmente, mas não é a pergunta ainda, queríamos saber a chance de ter pelo menos n pessoas AaBb nessa população com sei la quantos indivíduos, depende da geração, então queremos somar essas chances até o n que o exercício mandar.

Então vamos olhar aqui a probabilidade acumulada de vermos esses AaBb:

03

Aqui sim, observando a população com 25 pessoas por exemplo, se a gente quer saber a chance de pelo menos 1 pessoa ser AaBb, é praticamente 100%, como a gente tinha visto no gráfico anterior, era pouco provável que somente 1 pessoa fosse AaBb, logo ter pelo menos 1 pessoa AaBb tem que ser muito provável ora bolas :).
O mesmo esquema se repete para todas as distribuições, notem como elas parecem a mesma figurinha (mesma função) mesmo mudando o tamanho da população.

Agora sabemos os passos necessários, que são saber o tamanho da população, ver qual a chance de alguém ser AaBb e quantos caras desse tipo queremos, então entramos com os dados do exercício e usamos a função pbinom para saber qual a chance de termos n AaBb na população na geração k

k <- 2
n <- 1
N <- 2^k ## Tamanho da população na geração K
round(1 - pbinom( n-1, N, 0.25 ),digits=3)
[1] 0.684

E agora la vamos nos para o python, bem eu fui olhar as soluções, já que tenho que estudar mais python ainda, que até hoje não progredi nele 🙁
Mas uma solução legal é essa, onde definimos a pdf e a cdf, notem essas formulas do lado ali no wikipedia , é so usar, não precisa reinventar a roda, outra coisa seria usar de algum pacote pronto, mas é legal ver como elas são definidas.
Dai é so calcular o tamanho da população e retornar a probabilidade. Detalhe para a função map() no python, que é mais ou menos um apply do R para todos os itens

def binomialCoefficient(n,k):
    return factorial(n) / (factorial(k) * factorial(n-k))
 
def binompdf(n, p, r):
    return binomialCoefficient(n,r) * pow(p,r) * pow(1-p,n-r)
 
def binomcdf( n, p, r):
    return sum(map( lambda i: binompdf(n,p,i), range(r+1)))
 
def lia(k,r):
    n = 2*pow(2,k-1)
    return 1-binomcdf( n, .25, r-1 )

E por enquanto é isso :), eu achei esse problema foi um uso legal da distribuição binomial ^^

Codigo das figuras:

k <- 2
n <- 1
N <- 2^k ## Tamanho da população na geração K
 
#Figura 1
plot(x=0:10,y=c(0,2^c(1:10)),pch=19,xlab="Geração",
ylab="Tamanho da população",frame.plot=F,
main="Tamanho da população ao longo das gerações")
 
#Figura 2
par(mfrow=c(2,2))
plot(x=1:5,y=dbinom(c(1:5),5,prob=0.25),frame.plot=F,ylim=c(0,0.4),pch=19,
ylab="",xlab="",main="População com 5 pessoas",xlim=c(0,5))
plot(x=1:25,y=dbinom(c(1:25),25,prob=0.25),frame.plot=F,ylim=c(0,0.4),pch=19,
ylab="",xlab="", main="População com 25 pessoas",xlim=c(0,25))
plot(x=1:50,y=dbinom(c(1:50),50,prob=0.25),frame.plot=F,ylim=c(0,0.4),pch=19,
ylab="Probabilidade",xlab="Número de pessoas AaBb na População",
main="População com 50 pessoas",xlim=c(0,50))
plot(x=1:100,y=dbinom(c(1:100),100,prob=0.25),frame.plot=F,ylim=c(0,0.4),pch=19,
ylab="Probabilidade",xlab="Número de pessoas AaBb na População",
main="População com 100 pessoas",xlim=c(0,100))
 
#Figura 3
par(mfrow=c(2,2))
plot(x=1:5,y=c(1-pbinom(c(1:5),5,prob=0.25)),frame.plot=F,ylim=c(0,1),pch=19,
ylab="",xlab="",main="População com 5 pessoas",xlim=c(0,5))
plot(x=1:25,y=c(1-pbinom(c(1:25),25,prob=0.25)),frame.plot=F,ylim=c(0,1),
pch=19,ylab="",xlab="",main="População com 25 pessoas",xlim=c(0,25))
plot(x=1:50,y=c(1-pbinom(c(1:50),50,prob=0.25)),frame.plot=F,ylim=c(0,1),
pch=19,ylab="",xlab="",main="População com 25 pessoas",xlim=c(1,50))
plot(x=1:100,y=c(1-pbinom(c(1:100),100,prob=0.25)),frame.plot=F,ylim=c(0,1),
pch=19,ylab="",xlab="",main="População com 25 pessoas",xlim=c(1,100))
 
round(1pbinom( n-1, N, 0.25 ),digits=3)

Você consegue resolver isso?

Você consegue resolver isso?

 y=Y+2

Vi no 4chan e achei legal.
Legal pensar como a gente precisa de um mínimo de abstracção para resolver alguns problemas.

1362952035878

Jogando no WolframAlpha da para ter uma idéia da coisa 😉

Mais legal ainda é essa questão aqui:

pergunta

Tem uma discussão gigantesca sobre como responder essa pergunta no blog FlowingData, mas uma boa explicação que achei foi a do blog Understanding Uncertainty.

São boas perguntas para demonstrar como a gente é condicionado a tentar responder por associação, sempre tentando tomar atalhos, mesmo não notando, para as respostas…

E um brinde final, a gente sempre fala de Bugs por aqui, que é Bayesian Using Gibbs Sampler, mas outro algoritmo de mcmc é o Hamiltonian Monte Carlo (HMC) sampler, que é usado pelo Stan

Essa deveria ser a propaganda deles:

Um exemplo de Markov Chain

Usando estatística bayesiana, o termo MCMC é uma constante. MCMC é a sigla de Markov Chain Monte Carlo.
A gente fica usando isso direto, mas talvez valha a pena de vez em quando regredir um pouco e pensar sobre o que diabos são essas coisas.

Então a pergunta é:

“Que diabos é uma Markov Chain?”

Bem o nome é assim por causa do cara que inventou isso, que se chamava Andrey Markov. Se tudo ja é complicado, imagina matemática na Reversal Russia, mas resumindo ele viu que alguns processos estocásticos tinham uma propriedade interessante quando eram um processo sem memoria dos acontecimentos passados, e isso gerava alguns desdobramentos interessantes.

Antes de mais nada devemos lembrar que existem processos determinísticos, que são aquelas que sempre tem o mesmo resultado, o resultado é determinado. Por exemplo misturar cores, se misturar azul e amarelo da verde, sempre, toda vez, sempre, sempre, sempre, é determinado que isso vai acontecer, sempre, nunca muda, nunca, sempre verde. Isso é um processo determinístico.

Um processo estocástico é aqueles que a gente sabe quais podem ser os resultados, atribui chances para cada resultado, mas nunca podemos prever exatamente o resultado antes dele acontecer, por exemplo jogar um dado, a gente tem certeza que o resultado vai ser um número de 1 a 6, a gente tem certeza (se o dado não é roubado) que cada número tem uma chance de 1/6 de sair como resultado e ainda não temos ideia do resultado antes de jogar o dado. Logo processo estocástico é o processo que não é determinístico.

Certo, mas como diabos é um processo que não tem memoria? Vamos la pegar o exemplo do wikipedia de Markov Chain. Um modelo muito simples de previsão do tempo.

A idéia é a seguinte, todo dia, ou o dia esta chuvoso ou o dia está ensolarado, então temos 2 estados, sol ou chuva.

Não sabemos a chance do dia estar ensolarado ou não, mas sabemos que se esse é um dia de sol, temos 90% de chance de o próximo dia continuar ensolarado e 10% de chance do próximo dia chover.
Agora se o dia estivar chuvoso, o negocio funciona assim, temos 50% de chance do próximo dia ser chuvoso e 50% do próximo dia ser de sol.

Podemos resumir isso em uma figura que demonstra os estados e as transições possíveis.

mc

Bem simples não? Mas como essa figura é um grafo (eu acho que é), podemos representar ele na forma de uma matriz assim:

Sol Chuva Sol 0.9 0.1 Chuva 0.5 0.5

Aqui o nome da linha é o estado atual, e a coluna é a chance de ele ir para o estado que é o nome da coluna, é preciso estar vendo a figura acima nessa matriz.

Agora nos vamos fazer algumas simulações, a gente vai sortear um estado inicial qualquer, ou seja, se um dia esta sol ou chuva, e vamos simular dia após dia o que vai acontecer, de acordo com essa matriz. Ou seja se o dia começa com Sol, no dia 2 eu vou sortear 90% de chance de sol e 10% de chance de chuva, se choveu, no terceiro dia choveu, dai eu vou sortear o próximo dia como 50% de chance de chuva e 50% de chance de sol, e assim vai por 50 dias. Eu só não vou usar a matriz no primeiro dia porque não tem como, temos que começar de algum lugar.

Então para fazer 5 simulações eu faço uma matriz assim e sorteio o primeiro dia

> head(simul) [,1] [,2] [,3] [,4] [,5] [1,] "Chuva" "Chuva" "Chuva" "Chuva" "Sol" [2,] NA NA NA NA NA [3,] NA NA NA NA NA [4,] NA NA NA NA NA [5,] NA NA NA NA NA [6,] NA NA NA NA NA

Eu usei o comando head() aqui para a gente ver que ta funcionando as coisas, então eu tenho essa matriz ai (as 6 primeiras linhas demonstradas aqui), e para cada coluna eu vou seguir a próxima linha simulando as chances do estado de acordo com aquela matriz de mudanças.
Assim após a simulação, eu fico com isso:

> head(simul) [,1] [,2] [,3] [,4] [,5] [1,] "Chuva" "Chuva" "Chuva" "Chuva" "Sol" [2,] "Chuva" "Sol" "Sol" "Sol" "Sol" [3,] "Chuva" "Sol" "Sol" "Sol" "Sol" [4,] "Sol" "Sol" "Chuva" "Sol" "Sol" [5,] "Sol" "Sol" "Sol" "Sol" "Sol" [6,] "Chuva" "Sol" "Sol" "Sol" "Sol"

Então cada simulação saiu de um jeito, legal, é um processo estocástico, não deveria ser diferente, agora vamos olhar o resultado da primeira simulação com um gráfico.

01

Então olhando esse gráfico a gente vê que os primeiros dias foram de chuva, ai fez sol la pelo 4 e 5 dias, ai chuva denovo, então períodos de chuva e períodos de sol mais longos. E vemos mais dias de sol que dias de chuva no gráfico de baixo, que é o total de dia de chuva dividido por 50 (número de dias que observamos) e dias de sol também dividido por 50.

Esses números são exatamente:

Chuva Sol 0.24 0.76

24% de dias de chuva e 76% de dias de dias de sol. Mas agora vamos olhar as outras 4 simulações.

02

Apesar de cada processo ser de um jeito, dias de sol e dias de chuva, no final sempre há mais dias de sol que dias de chuva. Olhando a proporção entre eles vamos que:

Simulaçao 2 0.16 0.84 Simulaçao 3 Chuva Sol 0.2 0.8 Simulaçao 4 Chuva Sol 0.26 0.74 Simulação 5 Chuva Sol 0.2 0.8

Nossa, ficam parecidos não? E olha que so simulamos 50 dias, e não tinhamos nada estabelecido a priori que a chance de Sol deveria ser de 80%, ou algo por ai. Mas a gente pode fazer outra coisa, como quem faz as contas é o computador, vamos mandar ele simular dez mil dias (10000) e olhar o que acontece, vamos fazer exatamente o mesmo esquema, mas agora 10 mil vezes ao invés de 50.
Ah e vamos largar mão de sortear um início, vamos começar tudo com Chuva.

03

Agora ficou mais difícil de ver algo no gráfico, mas a gente ve que cada simulação continua diferente, apesar que la em baixo, olhando os gráficos de probabilidade (gráficos de barrinhas), todos parecem iguais, então vamos fazer a mesma conta acima, dividir os dias de sol e de chuva por 10 mil e ver o que acontece.

> table(simul[,1])/N Chuva Sol 0.1625 0.8375 > table(simul[,2])/N Chuva Sol 0.164 0.836 > table(simul[,3])/N Chuva Sol 0.1648 0.8352 > table(simul[,4])/N Chuva Sol 0.174 0.826

Eita, agora os números são mais iguais ainda. Algo em torno de 83% de chance de sol e 17% de chance de chuva. Mas o estranho de pensar é que esses números não servem para prever o dia 10001, eles são pior que aquelas chances que usamos para a simulação. Mas por ser um processo estocástico, podemos garantir que se observamos dias suficientes, cairemos sempre nessa probabilidades, nessa estabilidade. Essa é a propriedade que o Markov viu, e que é o princípio do MCMC. Quem manja de álgebra de matriz, tem como chegar a esses números com calculo, mas como eu malema consigo achar o limite de uma função, fico feliz em ver o negócio funcionando numa simulação por enquanto.

Varias coisas são processos com essa propriedade, jogos de tabuleiro que usam dados e, mais importante, valores aleatórios de distribuições estatísticas, que é onde entra o Monte Carlo do MCMC.

Por enquanto ficamos por aqui, mas acho que olhando melhor o que é uma Markov Chain da para ter uma idéia de como diabos é possível um algoritmo de Metropolis-Hasting sortear números e ainda sim estimar quanto deve ser um parâmetro de um modelo nosso.

#Independente da semente, o resultado é sempre o mesmo, pode testar :) 
set.seed(12345)
#Matriz de mudança de estado
matriz<-matrix(c(0.9,0.1,0.5,0.5),nrow=2,ncol=2,byrow=T,dimnames=
list(c("Sol","Chuva"),c("Sol","Chuva")))
matriz
 
#Estados possiveis
estado<-c("Sol","Chuva")
estado
 
#Matriz para receber os resultados das simulações
simul<-matrix(NA,ncol=5,nrow=50)
 
#Escolhendo um inicio ao acaso
simul[1,]<-sample(estado,5,replace = T)
head(simul)
 
#Loop para as simulações
for(j in 1:ncol(simul)) {
for(i in 1:(nrow(simul)-1)) {
simul[i+1,j]<-sample(estado,1,prob=matriz[simul[i,j],])
}
}
 
head(simul)
 
#Transformando em 0 e 1 apenas para fazer o grafico
simul2<-ifelse(simul=="Sol",1,0)
simul2
 
#Grafico 1
par(mfrow=c(2,1))
plot(1:nrow(simul2),simul2[,1],type="l",frame.plot=F,yaxt="n",xlab="Dias",
ylab="Tempo do Dia",main="Simulação 1")
axis(2,at=c(0,1),labels=c("Chuva","Sol"),las=2)
barplot(table(simul[,1])/50,ylab="Probabilidade",ylim=c(0,1))
 
table(simul[,1])/50
 
#Grafico 2
layout(matrix(1:8,ncol=4,nrow=2))
for(i in 2:ncol(simul)) {
plot(1:nrow(simul2),simul2[,i],type="l",frame.plot=F,yaxt="n",xlab="Dias",
ylab="Tempo do Dia",main=paste("Simulação ",i))
axis(2,at=c(0,1),labels=c("Chuva","Sol"),las=2)
barplot(table(simul[,i])/50,ylab="Probabilidade",ylim=c(0,1))
}
 
#Porcentagens
table(simul[,2])/50
table(simul[,3])/50
table(simul[,4])/50
table(simul[,5])/50
 
#Simulação com 10 mil dias
N<-10000
simul<-matrix(NA,ncol=5,nrow=N)
 
#Iniciando tudo com chuva, tente sua combinação
simul[1,]<-estado[2]
head(simul)
 
#Simulação
for(j in 1:ncol(simul)) {
for(i in 1:(nrow(simul)-1)) {
simul[i+1,j]<-sample(estado,1,prob=matriz[simul[i,j],])
}
}
 
simul2<-ifelse(simul=="Sol",1,0)
 
#Grafico 3
layout(matrix(c(rep(1:5,each=5),6:10),,ncol=5,nrow=6,byrow=T))
par(mar=c(1,4,1,1))
for(i in 1:ncol(simul)) {
plot(1:nrow(simul2),simul2[,i],type="l",frame.plot=F,yaxt="n",xlab="Dias",
ylab="Tempo do Dia",main=paste("Simulação ",i),lwd=0.05,xaxt="n")
axis(2,at=c(0,1),labels=c("Chuva","Sol"),las=2)
}
axis(1,cex.axis=0.5)
par(mar=c(5, 4, 2, 0) + 0.1)
barplot(table(simul[,1])/N,ylab="Probabilidade",ylim=c(0,1),yaxt="n")
axis(2,at=c(0,1))
for(i in 2:ncol(simul)) {
barplot(table(simul[,i])/N,yaxt="n",ylim=c(0,1))
}
 
#Porcentagens
table(simul[,1])/N
table(simul[,2])/N
table(simul[,3])/N
table(simul[,4])/N
table(simul[,5])/N