Markov Chain Monte Carlo

Aproveitando meu maior tempo livre nesse meses de agosto e setembro, irei iniciar uma série de notas sobre Markov Chain Monte Carlo, com base em materiais da internet. Assim, creio que essa série de posts será desinteressante para a esmagadora maioria dos meus 3 ou 4 leitores.

A série contará com códigos em R e, se eu tiver muito inspirado, códigos em C. A depender da minha paciência, posso colocar os links para vídeos e outros materiais que achar na internet. Mas a idéia é que quem se interessar possa acompanhar o conteúdo por aqui. Tentarei ser o mais didático possível.

É possível que, em algum momento, eu coloque todo o conteúdo em pdf, como notas de aula. Vamos ver.

******************

Monte Carlo

Para falar de Markov Chain Monte Carlo (MCMC) é preciso começar pela segunda parte do MC, que é Monte Carlo.

Aparentemente, físicos como Von Neuman e Ulan, envolvidos no projeto Manhatan (pra criar a bomba atômica) tinham que calcular integrais cabulosas. Como eram difíceis de serem calculadas, inventaram um método de aproximar os valores que envolve a simulação de número aleatórios. Por usar números aleatórios, batizaram esse método de Método de Monte Carlo, em homenagem aos jogos de azar do principado de Mônaco.

Assim, por exemplo, imagine que você quer calcular a seguinte integral:

\int^{a}_{b} h(x) dx

Se h(x) puder ser decomposta em h(x) = f(x) * p(x), em que p(x) é uma função densidade de probabilidade, então:

\int^{a}_{b} h(x)p(x) dx \approx E _{p(x)} [f(x)]

Ou seja, a integral pode ser aproximada pela esperança de f(x) com valores de x ocorrendo com probabilidade p(x). No caso mais simples em que

f(x)=x, temos que a esperança é simplesmente 1/n * \sum^{n}_{i=1}x_{i}, dado uma amostra independente de n observações de x_{1}, x_{2}, ..., x_{n}.

Em outras palavras, se for possível simular valores de x, então é possível computar a integral de interesse por meio de um simples cálculo de média.

Eu vou ilustrar agora em R, como calcular aproximadamente o valor de \pi por meio de simulação de Monte Carlo.

Como nós sabemos, a área de círculo, A, é dada por \pi r^{2}. Se considerarmos um círculo de raio r inscrito num quadrado de lado 2r, a área do quadrado é 4r^{2}. A razão entre a área do círculo m e a área do quadrado q, dado aqui por p, é \pi / 4. Assim, \pi = 4 * p.

Ora, sendo assim, é fácil ver que se eu puder gerar números aleatórios dentro do quadrado de raio 1, a proporção de números que cair dentro do círculo será a área do círculo. Ou ainda, de forma mais fácil, posso olhar apenas no primeiro quadrante, tal como na figura acima (um quarto de círculo inscrito num quadrado de raio 1). E como eu faço para gerar números aleatórios dentro do quadrado e checo se os números estão dentro do (um quarto) do círculo?

Simples, gere pares de números aleatórios (distribuição uniforme), x e y, no intervalo entre 0 e 1. E cheque se estão dentro do semi (semi círculo) testando a seguinte condição: x^{2} + y^{2} \leq 1

Abaixo segue o código em R que faz a simulação.

n = 10000
x = runif(n)
y = runif(n)
# figura

plot(y~x, pch=".")
curve(sqrt(1-x^2), 0,1, add=T, col="red")

m=sum(ifelse(x^2 + y^2 <= 1, 1, 0))
p = m/n
res.pi = p*4
res.pi

Created by Pretty R at inside-R.org

update: Faltou escrever a integral que permitiria calcular analiticamente a área do círculo e, portanto, \pi.

\pi = 4 \int \int I((x^{2} + y^{2}) < 1))* P(x,y) dx dy, em que I é uma função indicadora se a condição do parêntesis for satisfeita, e P(x,y) é a densidade uniforme conjunta de x e y no intervalo 0-1.

Sobre Manoel Galdino

Corinthiano, Bayesiano e Doutor em ciência Política pela USP.
Esse post foi publicado em estatística, Política e Economia e marcado , , , , , , , . Guardar link permanente.

3 respostas para Markov Chain Monte Carlo

  1. asesa disse:

    E onde entra a parte de markov chain?

  2. Ficou faltando a parte de MArkopv Chain. Talvez eu retome essas notas no segundo semestre.

  3. Fernando disse:

    Por favor termine de explicar a cadeia de Markov! Sua didatica é ótima!

Deixe uma resposta

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s