aula de hoje:

    1. o método de Monte Carlo
    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 distribuições multivariadas por MCMC
  • exercícios

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 \(u\) uniformemente distribuído entre 0 e 1, obter \(x\) \[ u = \int_{-\infty}^x P(x^\prime) dx^\prime = C(x)\]

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

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

\(x(u)\): função quantil\((u)\)

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

distribuição exponencial: \[ P(x) = e^{-x}, ~~~ x \ge 0 \] função cumulativa: \[ C(x) = \int_0^x e^{-x^\prime} dx^\prime = 1-e^{-x}\]

Note que esta distribuição tem média e variância iguais a 1.

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

Método 1: Monte Carlo: vamos amostrar 1000 valores de \(P(x)\) resolvendo numericamente \(C(x)-u = 0\):

# semente para assegurar a reprodutibilidade dos resultados
set.seed(1234)

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

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

# números aleatórios:
u = 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){
# defino uma função com a equação a ser resolvida
h = function(x){cum_prob(x)-u[i]}
# se obtém x resolvendo-se a equação h=0
# vamos supor que x está entre 0 e 100 e usar a função *uniroot()* do R
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='')

Método 2: Monte Carlo: vamos amostrar 1000 valores de \(P(x)\) resolvendo analiticamente \(C(x)-u = 0\):

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

Simulação:

# 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='')

Método 3: 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. 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

modelo: esfera com distribuição uniforme de pontos

parâmetros: N e R: N pontos dentro de uma esfera de raio R

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: \[ u_{r,i}, ~ u_{\theta,i}, ~ u_{\phi,i}\]

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

e \[x_i=r_i\sin(\theta_i) \cos(\phi_i)\] \[y_i=r_i\sin(\theta_i)\sin(\phi_i)\] \[z_i=r_i\cos(\theta_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”!

3. 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(N,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, 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
# intervalo de credibilidade: intervalo contendo uma certa percentagem dos valores do posterior
# exemplo: o intervalo de credibilidade de 95% é a região central do posterior que contém 95% dos valores
ic95 = quantile(Ysuper,probs=c(0.025,0.975))
ic95
##      2.5%     97.5% 
## -2.118080  5.987682

exemplo: 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)

# intervalo de credibilidade de 95%
ic95 = quantile(y[-(1:burnin)],probs=c(0.025,0.975))
ic95
##       2.5%      97.5% 
## 0.04713444 3.60291802

4. 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 3.987616 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
# intervalo de credibilidade de 95%
ic95x = quantile(cadeiamc[-(1:burnin),1],probs=c(0.025,0.975))
ic95x
##       2.5%      97.5% 
## -0.8235299  6.6370032
ic95y = quantile(cadeiamc[-(1:burnin),2],probs=c(0.025,0.975))
ic95y
##     2.5%    97.5% 
## 4.969951 8.876661

Exercícios

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

  1. Use o método de Monte Carlo para obter 1000 amostras da distribuição de probabilidades \(P(x) \propto x e^{-x}\) (\(x \ge 0\)), calculando a distribuição cumulativa numericamente. Ilustre o resultado com um histograma.
  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.