aula de hoje:

    1. o método de Monte Carlo
    1. as agulhas de Buffon
    1. simulação de Monte Carlo: distribuição uniforme de pontos em uma esfera
    1. Markov Chain Monte Carlo em uma dimensão
    1. simulação de uma distribuição exponencial por MCMC
    1. simulação de distribuições multivariadas por MCMC
  • exercícios
  • Apêndice:
  • A1 integração por Monte Carlo
  • A2 o viés de Malmquist-Eddington

1. o método de Monte Carlo

  • objetivo: obter amostras de uma distribuição de probabilidades unidimensional \(P(x)\)
  • dado um número aleatório \(\gamma\) uniformemente distribuído entre 0 e 1, obter \(x\) \[ \gamma = \int_{-\infty}^x P(x^\prime) dx^\prime = F(x)\]

\(F(x)\): função cumulativa de \(P(x)\)

\(x(\gamma)\) é obtido resolvendo-se a equação \(F(x)-\gamma = 0\).

\(x(\gamma)=quantil(\gamma)\)

ex.: simulação de uma distribuição exponencial

Vamos considerar 3 maneiras de se amostrar esta distribuição.

Para a distribuição exponencial (que é definida para \(x>0\)): \[ P(x) = e^{-x} \] e a função cumulativa é \[ F(x) = \int_0^x e^{-x^\prime} dx^\prime = 1-e^{-x}\]

  1. Monte Carlo: vamos amostrar 1000 valores de \(P(x)\) resolvendo \(F(x)-\gamma = 0\):
# semente para assegurar a reprodutibilidade dos resultados
set.seed(1234)

# número de simulações
nsim = 1000

# distribuição cumulativa:
cum_prob = function(x){1-exp(-x)}

# números aleatórios:
g = runif(nsim)

# amostragem por Monte Carlo
# inicialização da variável que vai conter os valores simulados
xsim = NULL

# simulações
for(i in 1:nsim){
h = function(x){cum_prob(x)-g[i]}
# se obtém x resolvendo-se a equação h=0
# vamos supor que x está entre 0 e 100
xsim[i] = uniroot(h, interval=c(0, 100))$root
# a saída de uniroot é uma lista e $root pega um dos elementos dessa lista, a raiz da equação
}

# estatísticas de x: os valores esperados para a média e o desvio padrão são 1 e 1
c(mean(xsim),sd(xsim))
## [1] 1.044677 1.061014
# distribuição dos valores simulados:
hist(xsim,col='gold',xlab='x',ylab='frequência',main='')

Para a distribuição exponencial, a equação a ser resolvida fica: \[ \gamma = \int_0^x e^{-x^\prime} dx^\prime = 1-e^{-x}\]
logo \[x = -\ln (1-\gamma)=-\ln \gamma^\prime\] onde \(\gamma^\prime\) é também um número aleatório uniformemente distribuído entre 0 e 1.

  1. Nesse caso, a simulação fica fácil:
# vamos gerar 1000 pontos como acima (por default log é o logarítmo natural):
xsim=-log(runif(1000))

# média e o desvio padrão
c(mean(xsim),sd(xsim))
## [1] 1.0233029 0.9879662
# distribuição dos valores simulados:
hist(xsim,col='gold',xlab='x',ylab='frequência',main='')

  1. Agora com rexp(): esta função gera números aleatórios distribuídos com uma pdf exponencial
xsim=rexp(1000)

# média e desvio padrão
c(mean(xsim),sd(xsim))
## [1] 1.020706 1.001654
# distribuição dos valores simulados:
hist(xsim,col='gold',xlab='x',ylab='frequência',main='')

2. As agulhas de Buffon

Em 1777 o Conde de Buffon calculou que uma agulha de comprimento L jogada ao acaso sobre um plano com linhas paralelas com separação T entre elas, com T > L, tem probabilidade \(2L/(\pi T)\) de cruzar uma linha.

Esta é também uma maneira de se estimar \(\pi\)!

seja x a distância do centro da agulha à linha paralela mais próxima e \(\theta\) o ângulo agudo entre a agulha e as linhas paralelas

a agulha cruza a linha se \(x \le L/2*sin(\theta)\)

https://en.wikipedia.org/wiki/Buffon%27s_needle_problem
https://en.wikipedia.org/wiki/Buffon%27s_needle_problem

vamos verificar isso fazendo uma simulação: para determinar a probabilidade prob de uma agulha jogada ao acaso cruzar uma linha, vamos jogar uma agulha \(n\) vezes e determinar a fração das jogadas que cruza uma linha

# função que calcula a probabilidade de cruzamento
# joga-se uma agulha n vezes e verifica-se a fração que cruza 

buffon <- function(n, L, T) {
# amostra o ângulo entre a agulha e as linhas paralelas
# theta varia entre 0 e pi mas, por simetria podemos amostrar apenas entre 0 e pi/2
  theta = runif(n, 0, pi/2)
# posição do meio da agulha: por simetria, distribuição uniforme entre 0 e T/2 
  x = runif(n, 0, T/2)
# cruza é um vetor true/false: verifica se x <= L/2*sin(theta)
cruza = (x <= L/2 * sin(theta))
# retorna quantas vezes cruza em n simulações:
prob = sum(cruza)/n
return(prob)
}

# vamos fazer uma simulação:
# reprodutibilidade
set.seed(123)
L = 1
T = 2
n = 100000
buffon(n, L, T)
## [1] 0.31934
# vamos agora estimar pi como 2L/(T*prob)
# vamos fazer 1000 simulações, cada uma com n jogadas, e determinar a distribuição de valores calculados de pi
PIcalc = NULL
for(i in 1:1000) {
PIcalc[i] = 2*L/(buffon(n, L, T)*T)
}

#estatísticas
mean(PIcalc)
## [1] 3.141249
quantile(PIcalc,  c(0.025, 0.975))
##     2.5%    97.5% 
## 3.114971 3.169070
#visualização
hist(PIcalc,main='',col='salmon',breaks=50)
abline(v=pi,col='blue',lwd=3)

3. simulação de Monte Carlo: distribuição uniforme de pontos em uma esfera centrada na origem

  • vamos simular pontos dentro de uma esfera de raio 1, com distribuição angular aleatória usando coordenadas esféricas
  • questão: como simular as coordenadas \((r, \theta, \phi)\) de um ponto?

vamos definir explicitamente todas as probabilidades envolvidas:

parâmetros: N pontos dentro de um raio R

densidade media: \(\bar n = 3N/(4 \pi R^3)\)

simulação das coordenadas esféricas \((r,\theta, \phi)\):

  • para o ponto \(i\) geramos 3 números aleatórios uniformemente distribuídos entre 0 e 1: \[ \gamma_{r,i}, \gamma_{\theta,i}, \gamma_{\phi,i}\]

e obtemos valores para \((r_i,\theta_i, \phi_i)\) como: \[r_i = R \gamma_{r,i}^{1/3}\] \[\theta_i = {\rm acos}(1- 2 \gamma_{\theta,i})\] \[\phi_i = 2 \pi \gamma_{\phi,i}\]

# biblioteca gráfica
#install.packages("scatterplot3d")
library(scatterplot3d)

set.seed(1)

# parâmetros do modelo
R=1
N=1000

#densidade media:
dm=3*N/(4*pi*R^3)

# coordenadas esfericas: (r,teta,fi)
teta=acos(1-2*runif(N))
fi=2*pi*runif(N)
r=R*(runif(N))^(1/3)

# coordenadas cartesianas
x=r*sin(teta)*cos(fi)
y=r*sin(teta)*sin(fi)
z=r*cos(teta)

# visualizacao
scatterplot3d(x,y,z,main=" ",asp=1)

# asp=1: (aspect) escala igual nos 3 eixos

Vamos ver se a distribuição de pontos é mesmo uniforme: vamos calcular o número esperado de pontos dentro de um raio r: \[ N(<r) = N \Big(\frac{r}{R}\Big)^3\]

par(mfcol=c(1,1)) 
#ordeno os raios, do menor para o maior
rs=sort(r)
#densidade media dentro do raio r
nn=seq(1:N)
dens.media=3*nn/(4*pi*rs^3)
plot(rs,nn,xlab="raio r",ylab="número cumulativo de pontos",type='l',col='blue',main = 'simulado x esperado',lwd=2)
#numero de pontos esperado dentro do raio r:
rr = seq(0,R,0.01)
nt=N*(rr/R)^3
lines(rr,nt,col="red",lwd=2)

Se voces repetirem algumas vezes a simulação com diferentes sementes iniciais (façam isso com N=100), verão que os gráficos acima variam bastante: isso é um exemplo da chamada “variância estatística”!

4. Markov Chain Monte Carlo

objetivo do MCMC: produzir amostras de uma distribuição de probabilidades \(P(x)\), onde \(x\) pode ser um vetor

exemplo: MCMC em uma dimensão: amostragem de uma gaussiana de média 2 e desvio padrão 2 (equivalente a rnorm(mean=2,sd=2))

Em geral, quando se trata de probabilidades, é melhor trabalhar com \(\log P(x)\) do que com \(P(x)\) diretamente. Vamos considerar primeiro o caso direto e depois o uso do logaritmo.

Vamos gerar amostras por MCMC para a distribuição de probabilidades \(P(x) \sim N(\mu=2,\sigma=2)\):

O primeiro passo é definir \(P(x)\):

p = function(x){dnorm(x, mean = 2, sd = 2, log = FALSE)}

e, em seguida, a função proposta: vamos supor que as propostas sejam também geradas por uma gaussiana, mas de média \(x\) e desvio padrão 1:

funcao_proposta = function(x){rnorm(1,mean = x, sd= 1)}

parâmetros do algoritmo de Metropolis:

  • valor inicial de \(x\): \(xini\)
  • número de iterações desejadas: \(niter\)
# mcmc: função que vai gerar a cadeia via MCMC com o algoritmo de Metropolis-Hastings
mcmc = function(xini, niter){
# defino o vetor que conterá os valores da cadeia  
    cadeia = array(niter+1)
# inicialização da cadeia    
    cadeia[1] = xini
# loop para gerar a cadeia    
    for (i in 1:niter){
# obtenho uma proposta      
        proposta = funcao_proposta(cadeia[i])
# razão entre o valor proposto nesta iteração e o valor anterior da cadeia        
        probab = p(proposta)/p(cadeia[i])
# gero um número aleatório uniformemente distribuído entre 0 e 1 
      runif1=runif(1)
# testo se a proposta é aceita ou não     
        if (probab > runif1) {
# aceita          
            cadeia[i+1] = proposta}
        else{
# não aceita: usa o valor anterior da cadeia         
            cadeia[i+1] = cadeia[i]}
    }
# retorna a cadeia    
    return(cadeia)
}

Esta função devolve uma cadeia com \(niter+1\) valores \(x_i\).

Exemplo:

#inicialização
set.seed(1234)

# vamos começar "longe" do valor médio esperado
xini=-10
niter=1000

#mcmc
y = mcmc(xini,niter)

#visualização da cadeia:
par(mfrow = c(1,2))
x=seq(1,niter+1,1)
plot(x,y,type='l',col='red',ylab='x',xlab='iteração',main='cadeia')

# adiciona uma reta horizontal no valor esperado
abline(h = 2, col="blue",lwd=2)

hist(y,col='green',main='distribuição de x',xlab='x',ylab='frequência')

# adiciona uma reta vertical no valor esperado
abline(v = 2, col="blue",lwd=2)

note que, nas primeiras iterações, o valor de x começa longe do valor médio esperado mas depois de umas \(\sim 50-100\) iterações, fica dentro da região esperada- é isso que justifica o burn-in, a remoção de um certo número de iterações iniciais quando se vai calcular as estatísticas da cadeia.

log de P(x)

Muitas vezes é conveniente se trabalhar com o logaritmo da distribuição de probabilidades. Isso se deve pois em muitos casos os valores de \(P(x)\) podem ficar muito pequenos e causar problemas numéricos.

Vamos então considerar MCMC com o log de \(P(x)\). Notem a diferença na definição da função:

p = function(x){return(dnorm(x, mean = 2, sd = 2, log = TRUE))}

e no algoritmo mcmc:

# vamos chamar mcmc1 a função que vai fazer MCMC com o logaritmo de P
mcmc1 = function(xini, niter){
    cadeia = array(niter+1)
    cadeia[1] = xini
    for (i in 1:niter){
        proposta = funcao_proposta(cadeia[i])
# como estou usando o log(p):        
        probab = exp(p(proposta) - p(cadeia[i]))
      runif1=runif(1)
        if (probab > runif1) {
            cadeia[i+1] = proposta}
        else{
            cadeia[i+1] = cadeia[i]}
    }
    return(cadeia)
}

e, rodando de novo:

#inicialização
set.seed(123)
xini=-10
niter=1000

#mcmc
y = mcmc1(xini,niter)

#visualização da cadeia:
par(mfrow = c(1,2))
x=seq(1,niter+1,1)
plot(x,y,type='l',col='red',ylab='x',xlab='iteração',main='cadeia')
abline(h = 2, col="blue",lwd=2)
hist(y,col='green',main='distribuição de x',xlab='x',ylab='frequência')
abline(v = 2, col="blue",lwd=2)

Burn-in

Um certo número de cadeias iniciais (denominadas burn-in) devem ser removidas da análise por não terem ainda atingido uma “região de convergência”. Por exemplo, a figura abaixo mostra as 200 iterações iniciais de 4 cadeias com condições iniciais diferentes:

#inicialização
ncadeia=4
niter=200
xinic=c(-10,0,10,50)

#inicialização de variáveis
y=rep(0,niter+1)
Y=matrix(nrow=niter+1,ncol=ncadeia)

#mcmc
for (i in 1:ncadeia){
y=mcmc1(xinic[i],niter)

# combina  as colunas do vetor y (cadeias)
Y[,i]=cbind(y)
}

#visualização da cadeia:
par(mfrow = c(1,1))
x=seq(1,niter+1,1)
plot(x,Y[,1],type='l',col='red',ylab='x',xlab='iteração',ylim=c(-11,51))
lines(x,Y[,2],type='l',col='green',ylab='x',xlab='iteração')
lines(x,Y[,3],type='l',col='blue',ylab='x',xlab='iteração')
lines(x,Y[,4],type='l',col='black',ylab='x',xlab='iteração')
abline(h = 2, col="blue" )

Nesse exemplo, vemos que as primeiras ~150 iterações devem ser descartadas como burn-in.

Análise

Vamos supor que fizemos uma ou mais simulações e a cadeia convergiu.

Como podemos exibir os resultados? A primeira coisa a se fazer é definir quanto vamos eliminar da cadeia como burn-in. Com o resto da cadeia podemos:

  • determinar a distribuição de \(x\)
  • saber a taxa de aceitação do algoritmo, isto é, a fração de propostas que foram aceitas
  • calcular estatísticas de interesse: média, mediana, variância, quartis…

Vamos fazer uma simulação:

#inicialização
set.seed(1234)
xini=-10
niter=10000

#mcmc
y = mcmc1(xini,niter)

#visualização da cadeia:
par(mfrow = c(1,2))
x=seq(1,niter+1,1)
plot(x,y,type='l',col='red',ylab='x',xlab='iteração',main='cadeia')
abline(h = 2, col="blue" )
hist(y,col='green',main='distribuição de x',xlab='x',ylab='frequência')
abline(v = 2, col="blue" )

Vamos jogar fora, conservadoramente, as 1000 primeiras iterações como burn-in:

burnin = 1000

#taxa de aceitacao: fração dos valores propostos que são aceitos
aceitacao = 1-mean(duplicated(y[-(1:burnin)]))
aceitacao
## [1] 0.8403511
#media, mediana, desvio padrão
mean(y[-(1:burnin)])
## [1] 2.062565
median(y[-(1:burnin)])
## [1] 2.094431
sd(y[-(1:burnin)])
## [1] 2.002334
#intervalo de confiança de 95%
quantile(y[-(1:burnin)], probs = c(0.025,0.975))
##      2.5%     97.5% 
## -1.983647  5.904849

Convergência

Vamos agora analisar a convergência da cadeia. Para isso vamos jogar fora as iterações de burn-in; em seguida dividimos a cadeia restante em 4 segmentos, calculamos a média e a dispersão dos 4 segmentos e avaliamos.

seg=(niter-burnin)/4
media.seg = c(0,0,0,0)
sd.seg = c(0,0,0,0)
for (i in 1:4){
i1=burnin+1+(i-1)*seg
i2=i1+seg
media.seg[i] = mean(y[i1:i2])
sd.seg[i] = sd(y[i1:i2])
}
#media de cada segmento da cadeia:
media.seg
## [1] 2.109344 2.153425 1.883418 2.105790
#desvio padrão de cada segmento da cadeia
sd.seg
## [1] 1.970430 1.881251 2.054113 2.087099
#media das cadeias
mean(media.seg)
## [1] 2.062994
#desvio padrão das cadeias
sd(media.seg)
## [1] 0.1216623

a semelhança entre as estatísticas de cada segmento sugere que a cadeia convergiu

Um outro procedimento muito comum (principalmente em problemas complexos) é rodar várias cadeias e compará-las, para avaliar a convergência.

Exemplo: vamos comparar resultados de 4 cadeias

#inicialização
ncadeia=4
niter=11000
xinic=c(-1,0,3,5)
#inicialização de variáveis
y=rep(0,niter+1)
Y=matrix(nrow=niter+1,ncol=ncadeia)
#mcmc
for (i in 1:ncadeia){
y=mcmc1(xinic[i],niter)
Y[,i]=cbind(y)}
#visualização da cadeia:
par(mfrow = c(1,1))
x=seq(1,niter+1,1)
plot(x,Y[,1],type='l',col='red',ylab='x',xlab='iteração')
lines(x,Y[,2],type='l',col='green',ylab='x',xlab='iteração')
lines(x,Y[,3],type='l',col='blue',ylab='x',xlab='iteração')
lines(x,Y[,4],type='l',col='black',ylab='x',xlab='iteração')
abline(h = 2, col="yellow",lwd=3)

Podemos, também, calcular algumas estatísticas para essas cadeias (assumindo aqui o mesmo burn-in de antes para todas):

# media de cada cadeia:
media=rep(0,ncadeia)
for (i in 1:ncadeia){media[i]=mean(Y[-(1:burnin),i])}
media
## [1] 1.938203 1.940195 1.981517 2.003466
# desvio padrão de cada cadeia:
dp=rep(0,ncadeia)
for (i in 1:ncadeia){dp[i]=sd(Y[-(1:burnin),i])}
dp
## [1] 2.089776 2.056148 2.138996 1.971478

e, se quisermos, podemos combinar as diferentes cadeias numa “super-cadeia”:

# super cadeia: combinação de todas as cadeias sem o burn-in
Ysuper=c(Y[-(1:burnin),1],Y[-(1:burnin),2],Y[-(1:burnin),3],Y[-(1:burnin),4])
# super media e desvio padrão
mean(Ysuper)
## [1] 1.965845
sd(Ysuper)
## [1] 2.065111

5. simulação por MCMC de uma distribuição exponencial

A distribuição exponencial é definida como \[ P(x) = \left \{ \begin{matrix} e^{-x}, & \mbox{se }x \ge 0 \\ 0, & \mbox{se } x<0 \end{matrix} \right. \]

Para amostrar uma função por MCMC é preciso definir a distribuição de probabilidades da função a ser amostrada e a função de propostas.

# reprodutibilidade
set.seed(1234)
# função a ser amostrada
p = function(x) dexp(x, log = FALSE)

A função proposta deve propor valores dentro do espaço amostral da distribuição de probabilidades considerada.

A distribuição exponencial é definida para \(x> 0\), enquanto que uma função proposta gaussiana pode ter valores negativos.

Assim, precisamos modificar a função proposta gaussiana para ela propor apenas valores positivos:

funcao_proposta = function(x){
a = -1
while(a < 0){a = rnorm(1,mean = x, sd= 1)}
return(a)}

Vamos começar “longe” do valor esperado:

xini=10
niter=2000

#vamos usar a função mcmc sem o logaritmo de P:
y = mcmc(xini,niter)

#visualização da cadeia:
par(mfrow = c(1,2))
x=seq(1,niter+1,1)
plot(x,y,type='l',col='red',ylab='x',xlab='iteração',main='cadeia')
# adiciona uma reta horizontal com o valor esperado da média (1)
abline(h = 1, col="blue" )
hist(y,col='green',main='distribuição de x',xlab='x',ylab='frequência',breaks=20)
# adiciona uma reta vertical com o valor esperado da média (1)
abline(v = 1, col="blue",lwd=3)

# estatísticas da amostra simulada, sem o burn-in
# vamos supor um burn-in de 200 iterações
burnin=200

# média e desvio padrão da amostra
mean(y[-(1:burnin)])
## [1] 1.191461
sd(y[-(1:burnin)])
## [1] 0.9768215
# compare com os valores esperados (1 e 1)

6. Simulações de distribuições multivariadas

MCMC é muito útil também para simular distribuições multivariadas.

Vamos ilustrar isso com a simulação de pares (x,y) de acordo com uma distribuição gaussiana bivariada com coeficiente de correlação \(\rho\): \[ P(x,y|\mu_x,\mu_y,\sigma_x,\sigma_y,\rho) = {1 \over 2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \times \] \[ \exp\Bigg[{-1 \over 2(1-\rho^2)} \Bigg({(x-\mu_x)^2 \over \sigma_x^2} + {(y-\mu_y)^2 \over \sigma_y^2} -{2\rho (x-\mu_x) (y-\mu_y) \over \sigma_x \sigma_y} \Bigg) \Bigg] \]

  • se \(\rho \approx 0\) a correlação é nula ou fraca; se \(\rho \approx \pm 1\) a correlação/anti-correlação é forte.
  • aqui definimos os parâmetros da distribuição que vamos simular:
rho=0.5
media.x=3
sd.x=2
media.y=7
sd.y=1

# constantes
c1=1/(2*pi*sd.x*sd.y*sqrt(1-rho^2))
c2=1/(2*(1-rho^2))
c3=2*rho/(sd.x*sd.y)
  • Vamos gerar um conjunto de pontos amostrados de uma gaussiana bivariada por MCMC
  • Ao invés de x e y vamos considerar um vetor com as variáveis da distribuição: param[1]=x e param[2]=y
npar=2

#função a ser amostrada
p = function(param){
z=(param[1]-media.x)^2/sd.x^2-c3*(param[1]-media.x)*(param[2]-media.y)+(param[2]-media.y)^2/sd.y^2
pxy=c1*exp(-c2*z)

# vamos trabalhar com o logaritmo da distribuição de probabilidades
return(log(pxy))
}

# função proposta gaussiana para cada parâmetro
funcao_proposta = function(param){rnorm(2,mean = param, sd= c(0.1,0.3))}

# mcmc adaptado para mais de uma variável
mcmc2 = function(xini, iterations){
    cadeia = array(dim = c(iterations+1,npar))
    cadeia[1,] = xini
    for (i in 1:iterations){
        proposta = funcao_proposta(cadeia[i,])
# como a probabilidade está em log         
        probab = exp(p(proposta) - p(cadeia[i,]))
        if (runif(1) < probab){
            cadeia[i+1,] = proposta}
        else{
            cadeia[i+1,] = cadeia[i,]}
    }
    return(cadeia)
}
  • Agora inicializamos e rodamos a simulação:
# para medir o tempo de processamento deste módulo
t0 = Sys.time()

xini = c(0,10)
set.seed(234234)
cadeiamc = mcmc2(xini, 100000)

# burn-in
burnin = 1000
aceitacao = 1-mean(duplicated(cadeiamc[-(1:burnin),]))
aceitacao
## [1] 0.8886072
# visualização sem o burn-in
# notem como se pode remover o burn-in da cadeia
par(mfrow = c(2,2))
hist(cadeiamc[-(1:burnin),1],nclass=30, , main="distribuição de x", xlab="valor do input = linha vermelha",col='gold')
abline(v = mean(cadeiamc[-(1:burnin),1]),lwd=2)
abline(v = media.x, col="red",lwd=2)
hist(cadeiamc[-(1:burnin),2],nclass=30, main="distribuição de y", xlab="valor do input = linha vermelha",col='gold')
abline(v = mean(cadeiamc[-(1:burnin),2]),lwd=2)
abline(v = media.y, col="red",lwd=2)
plot(cadeiamc[-(1:burnin),1], type = "l", xlab="valor do input = linha vermelha" , main = "Cadeia de  valores de x", col='gray27')
abline(h = media.x, col="red",lwd=2)
plot(cadeiamc[-(1:burnin),2], type = "l", xlab="valor do input = linha vermelha" , main = "Cadeia de  valores de y", col='gray27')
abline(h = media.y, col="red",lwd=2)

# essa figura sugere que a cadeia de $x$ pode não ter convergido e que deve-se aumentar o número de simulações

# outra visualização
par(mfrow = c(1,1))
plot(cadeiamc[-(1:burnin),1],cadeiamc[-(1:burnin),2],xlab="x", ylab="y" , col="red", type = "l")

# tempo de processamento deste módulo
Sys.time()-t0
## Time difference of 4.53287 secs
  • Outra visualização das amostras simuladas:
library(KernSmooth)
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
ca=cadeiamc[-(1:burnin),1]
cb=cadeiamc[-(1:burnin),2]
smoothScatter(ca,cb,nrpoints=0,add=FALSE,xlab="x",ylab="y")
# contornos
Ndim = 101
z = cbind(ca,cb)
xmax = max(ca)
xmin = min(ca)
ymax = max(cb)
ymin = min(cb)
deltax = (xmax-xmin)
deltay = (ymax-ymin)
bw.x = 6*deltax/Ndim
bw.y = 6*deltay/Ndim
dens = bkde2D(z,bandwidth=c(bw.x,bw.y),gridsize=c(Ndim,Ndim)) 
contour(dens$x1,dens$x2,dens$fhat,add=T,drawlabels=F,nlevels=5)
abline(v = media.x, col="red",lwd=3)
abline(h = media.y, col="red",lwd=3)

  • Os valores médios dos x e y simulados podem ser comparados com os parâmetros iniciais:
# valor médio de x
mediax=mean(cadeiamc[-(1:burnin),1])
mediax
## [1] 2.734159
# valor inicial da média de x
media.x
## [1] 3
# valor médio de y
mediay=mean(cadeiamc[-(1:burnin),2])
mediay
## [1] 6.925175
# valor inicial da média de y
media.y
## [1] 7

Exercícios

Para cada exercício inicialize o gerador de números aleatórios com seu número USP.

  1. Estime a integral \(I = \int_0^1 x^3 dx\) via MC usando n=10,100,1000 e 10000 pontos. Para estudar o erro da média de estimativas deste tipo, faça \(N=10000\) simulações para cada \(n\) e determine a variância nos valores simulados de \(I\). Faça um gráfico log-log dessa variância versus \(n\) e verifique se é a variância é proporcional a \(1/n\).
  1. Simule por MCMC a distribuição de probabilidades \(P(x) \propto x e^{-x}\) (\(x \ge 0\)). Faça simulações com 10000 iterações. Use uma função de aceitação gaussiana com \(\sigma=\) 0.1, 1 e 10. Faça figuras relevantes e comente como a escolha de \(\sigma\) afeta o algoritmo. Comente sobre o ‘burn-in’ nesses casos.

**Apêndice.

integração por Monte Carlo

MC oferece uma forma simples para se integrar uma função positiva \(f(x)\): \[ I = \int_a^b f(x) dx,~~~~f(x)>0\]

Algoritmo do método da rejeição: vamos supor que a função \(f(x)\) possa ser ``coberta’’ por uma função mais simples \(g(x)\) (por exemplo \(g(x)=1\) no intervalo [0,1]), de área \(A\) no intervalo de integração

  • sorteamos N pontos uniformemente dentro de A
  • contamos quantos pontos cairam “dentro” de \(f(x)\), \(N_{ac}\)
  • integral: \[ I \simeq {N_{ac} \over N} A\]
Sobol, A Primer for the Monte Carlo Method (1968…)
Sobol, A Primer for the Monte Carlo Method (1968…)

integração em uma dimensão

qual é a integral de uma função gaussiana, de média 2 e desvio padrão 10, entre 3 e 6? \[ \int_3^6 N(\mu = 2, \sigma = 10)dx \]

usando R, essa integral pode ser calculada diretamente a partir da função cumulativa da gaussiana:

pnorm(6,mean=2,sd=10)-pnorm(3,mean=2,sd=10)
## [1] 0.1155939

podemos também estimar essa integral gerando um número muito grande de amostras de \(N(\mu = 2, \sigma = 10)\) e determinando a fração de pontos que caem entre 3 e 6:

# reprodutibilidade
set.seed(123)

# vamos considerar simulações com 100, 10000 e 1000000 de pontos:
nsim=c(10^2,10^4,10^6)

# simulações
for(i in 1:3){

# número de pontos na simulação
n = nsim[i]

# amostragem da gaussiana
x=rnorm(n,mean=2,sd=10)

# estimativa da integral: fração de pontos no intervalo desejado
integral = length(x[x >= 3 & x <= 6])/n
print(integral)
}
## [1] 0.13
## [1] 0.1157
## [1] 0.115726
# e também podemos fazer integração numérica para comparação
f = function(x) dnorm(x,mean=2,sd=10)
integral = integrate(f,lower=3,upper=6)
integral
## 0.1155939 with absolute error < 1.3e-15

cálculo de \(\pi\) por integração de Monte Carlo

Vamos calcular \(\pi\) pela área de um quarto de círculo unitário (\(A=\pi/4\)):

# inicialização
set.seed(123)

# número de pontos na simulação
npt=1000

# simulação de x e y em um quadrado de lado 1
x=runif(npt)
y=runif(npt)

# distância de (x,y) ao centro (0,0)
d=sqrt(x^2+y^2)

# número de pontos 'aceitos', que caem dentro do círculo
nac=length(which(d<=1)) 

# valor de pi na simulação: 4 vezes a fração de pontos dentro da área do quarto de círculo
meu_pi=4*nac/npt 
meu_pi
## [1] 3.22
#visualização
plot(x[which(d<=1)],y[which(d<=1)],xlab="x",ylab="y",main=sprintf(paste("pi por Monte Carlo ", npt," pontos")),col='red',asp=1,pch=20)
points(x[which(d>1)],y[which(d>1)],col='blue',pch=20)
xl=seq(0,1,0.01)
yl=sqrt(1-xl^2)
lines(xl,yl,col='green',lwd=3)
abline(v=c(0,1))
abline(h=c(0,1))

O bias de Malmquist-Eddington**

1 contagens de galáxias

uma importante informação em estudos extragalácticos é a contagem de galáxias em função da magnitude

  • contagens cumulativas: número de galáxias por grau quadrado com magnitudes mais brilhantes que uma certa magnitude \(m\);
  • contagens diferenciais: número de galáxias por grau quadrado com magnitudes entre \(m\) e \(m+dm\) as contagens diferenciais podem ser obtidas como a derivada das contagens cumulativas
  • em 1936 Hubble propôs usar contagens para testar modelos cosmológicos
http://star-www.dur.ac.uk/~nm/pubhtml/counts/counts.html
http://star-www.dur.ac.uk/~nm/pubhtml/counts/counts.html

exemplo: universo euclidiano e uniforme

  • suponha que as galáxias têm todas a mesma luminosidade, \(L\), e estão uniformemente distribuidas com densidade volumétrica \(n\)
  • queremos \(N(<m)\): número de galáxias mais brilhantes que a magnitude \(m\)
  • como supomos que todas as galáxias têm a mesma luminosidade, este número será igual ao número de galáxias a uma distância menor que um certo \(d\)
  • número de galáxias a uma distância \(d\): \[ N(<d) = {4 \over 3} \pi d^3 n \]
  • fluxo de uma galáxia a uma distância \(d\): \[ f = {L \over 4 \pi d^2} \]
  • magnitude aparente dessa galáxia: \[ m = m_0 -2.5 \log f = {\rm cte}+5 \log d \] e, portanto, \[ d \propto 10^{0.2 m} \]
  • como todas as galáxias \(N(<d)\) terão magnitudes menores que \(m\) e \(N(<d) \propto d^3\), \[ N(<m) \propto 10^{0.6 m} \] este é o comportamento esperado das contagens em um universo euclidiano uniforme
m = seq(14,28,0.1)
contagens = 2*10^(0.6*(m-14))
plot(m,contagens,type='l',log='y',main='modelo euclidiano uniforme',col='salmon',lwd=3)

# neste modelo, o log das contagens varia linearmente com as magnitudes
  • vocês podem verificar que se este universo é infinito, o fluxo que chega no observador também é infinito: a noite não poderia ser escura (o paradoxo de Olbers, 1826)
  • note que nesse caso tanto as contagens cumulativas quanto as diferenciais têm a mesma forma!

2 o viés de Malmquist-Eddington

aplicação do conceito de probabilidades condicionais

viés de um estimador: a diferença entre o valor esperado e o valor medido

viés de M-E:

  • como a distribuição das contagens é muito não uniforme (as contagens crescem rapidamente conforme \(m\) aumenta) e como as magnitudes têm erros, uma amostra limitada em magnitudes tem um excesso de galáxias com magnitudes intrínsecas mais fracas que a magnitude limite
  • em outras palavras: é mais provável uma galáxia mais fraca que a magnitude limite ficar com magnitude menor (mais brilhante) que a magnitude limite que o contrário
  • vamos supor que \(P(m)\) seja a probabilidade de que uma galáxia tenha magnitude aparente intrínseca \(m\); vamos adotar a probabilidade esperada para uma distribuição uniforme, euclidiana: \[ P(m) \propto 10^{0.6 m} \]
  • \(P(m_{obs}|m)\): probabilidade de que uma galáxia de magnitude intrínseca \(m\) seja observada com magnitude \(m_{obs}\) devido aos erros fotométricos
  • usando a regra do produto, \[ P(m_{obs},m) = P(m_{obs}|m)P(m) \] e usando a regra da soma \[ P(m_{obs}) = \int P(m_{obs}|m)P(m)dm \]
  • para erros gaussianos,
    \[P(m_{obs}|m) = f(m-m_{obs}) \propto {\rm exp} \Big[{(m-m_{obs})^2 \over 2 \sigma^2} \Big],\] e, portanto, \[ P(m_{obs}) = \int f(m-m_{obs})P(m)dm \] que é equivalente a uma convolução da distribuição verdadeira de magnitudes com um kernel gaussiano.

para ver como isso leva ao viés de M-E, vamos simular esse problema, examinando a distribuição que seria observada, produzida por uma população com magnitude verdadeira no intervalo \(0 < m < m_{lim} = 10\) com erros nas magnitudes que vamos supor crescentes com a magnitude: \[ \sigma = 0.1+0.05*m \]

# reprodutibilidade: 
# devemos estabelecer a semente do gerador de números aleatórios para termos resultados reprodutíveis
set.seed(123)

# magnitude limite
mlim = 10

# distribuição de probabilidades não normalizada 
prob = function(x){10^(0.6*x)}

# distribuição de probabilidades normalizada 
probn = function(x){prob(x)/integrate(prob,lower= 0, upper= mlim)$value}

# distribuição cumulativa
cum_prob = function(x){integrate(probn,lower= 0, upper= x)$value}
 
# número de simulações 
nsim = 100000

# inicialização de variáveis: magnitude verdadeira e observada
xv = rep(0,nsim)
xo = rep(0,nsim)

# vetor com números aleatórios entre 0 e 1
g = runif(nsim,min=0,max=1)

# vamos ver quanto tempo demora para rodar esta simulação
t0 = Sys.time()

for(i in 1:nsim){
# amostragem das magnitudes por Monte Carlo
# a magnitude é aquela cujo valor corresponde à probabilidade g na distribuição cumulativa
h = function(x){cum_prob(x)-g[i]}
# que se obtém resolvendo-se a equação h=0
xv[i] = uniroot(h, interval=c(0, 10))$root

# vamos adotar erros gaussianos, com um desvio padrão que aumenta com a magnitude:
xo[i] = xv[i] + rnorm(1,mean = 0, sd = 0.1+0.05*xv[i])
}

# tempo de processamento:
Sys.time()-t0
## Time difference of 48.5653 secs
# sumário das magnitudes observadas simuladas
summary(xo)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.809   8.781   9.380   9.279   9.896  12.382
# visualização
par(mfrow=c(1,1))
hist(xv,xlim=c(5,11),ylim=c(0,25000),breaks=50,col=rgb(0,0,1,0.3),main="",xlab='x',cex.lab=1.5)
par(new=TRUE)
hist(xo,xlim=c(5,11),ylim=c(0,25000),breaks=50,col=rgb(1,0,1,0.3),main="",xlab="",ylab="")
legend(5,25000,legend='100,000 simulations',bty="n",cex=1.5)
legend(5.5,20000,legend=c('true','observed'),col=c(rgb(0,0,1,0.3),rgb(1,0,1,0.3)),  text.col = c(rgb(0,0,1,0.3),rgb(1,0,1,0.3)),bty="n",cex=1.5)

note que para \(m < 9\) as contagens observadas estão acima das verdadeiras:

naturalmente, devido ao desenho da simulação, as contagens não são confiáveis próximas a mlim (devido à não inclusão de galáxias simuladas mais fracas)