aula5

aula de hoje:

o método de Monte Carlo

amostragem por aceitação-rejeição

simulação de Monte Carlo: distribuição uniforme de pontos em uma esfera

Markov Chain Monte Carlo em uma dimensão

simulação de distribuições multivariadas por MCMC

exercícios da semana

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]}

# x é obtido resolvendo-se a equação h=0

# vamos supor que x está entre 0 e 100 e usar a função *uniroot()* do R 
# para resolver a equação
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. amostragem por aceitação-rejeição

vamos exemplificar o procedimento amostrando a função beta com alfa = 6 e beta = 3.

# função que queremos amostrar:
funcao_f <- function(x) {
  return(dbeta(x, shape1 = 6, shape2 = 3))
}

x <- seq(0,1,0.01)
y <- funcao_f(x)
plot(x,y,type='l',col='salmon',lwd=3,ylab='f(x)')

vamos adotar como função proposta a distribuição uniforme no intervalo [0, 1]:

# vamos adotar a distribuição uniforme no intervalo [0, 1]
funcao_proposta_g <- function(x) {
  return(dunif(x, min = 0, max = 1))
}

# gerador de números aleatórios para a distribuição proposta
g_aleatorio <- function(n) {
  return(runif(n, min = 0, max = 1))
}

vamos encontrar uma constante M para cobrir a função que se quer amostrar, f(x): \(M g(x) \ge f(x)\).

para isso podemos achar o máximo de f(x) com a função optimize() no intervalo [0, 1]:

max_f <- optimize(f = funcao_f, interval = c(0, 1), maximum = TRUE)$objective
max_f

[1] 2.549958

# vamos adotar um valor um pouco maior que o máximo:
M <- 2.56

vamos agora criar uma função que implementa o algoritmo aceitação-rejeição:

amostragem_aceitacao_rejeicao <- function(n_amostras) {
  amostras <- numeric(n_amostras)
  n_aceitas <- 0
  n_tentativas <- 0
  
  while (n_aceitas < n_amostras) {
    # gera um candidato a partir de g(x)
    z <- g_aleatorio(1)
    # gera um número aleatório uniformemente distribuído u
    u <- runif(1)
    
    # Calcula a probabilidade de aceitação: f(z) / (M * g(z))
    p_aceitacao <- funcao_f(z) / (M * funcao_proposta_g(z))
    
    # Aceita ou Rejeita:
    if (u <= p_aceitacao) {
      amostras[n_aceitas + 1] <- z
      n_aceitas <- n_aceitas + 1
    }
    n_tentativas <- n_tentativas + 1
  }
  
# eficiência do algoritmo:
print(paste("número total de amostras", 
n_amostras, " número de tentativas:", 
n_tentativas))
print(paste("taxa de aceitação:", round(n_amostras / n_tentativas * 100, 2), "%"))
  
  return(amostras)
}

vamos gerar 1000 amostras de f(x):

# reprodutibilidade
set.seed(123)

amostras <- amostragem_aceitacao_rejeicao(1000)

[1] "número total de amostras 1000  número de tentativas: 2467"
[1] "taxa de aceitação: 40.54 %"

# visualização com histograma
hist(amostras, breaks = 30, main = "Amostras por Aceitação-Rejeição", 
xlab = "x",ylab='frequência')

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

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)) 
# distribuição cumulativa observada
# ordeno os raios, do menor para o maior
rs=sort(r)
nn=seq(1:N)
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 distribuição beta com parâmetros de forma alfa = 6 e beta = 3

x <- seq(0,1,0.01)
plot(x,dbeta(x, shape1 = 6, shape2 = 3),type='l',col='salmon',lwd=3,ylab='f(x)')

Em geral, quando se trata de probabilidades, é melhor trabalhar com \(\log P(x)\) do que com \(P(x)\) diretamente, 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)\).

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

P = function(x){return(dbeta(x, shape1 = 6, shape2 = 3, log = TRUE))}

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

Mas note que os valores propostos devem ser consistentes com o intervalo permitido para a distribuição beta: [0,1]

Assim, não consideraremos valores fora desse intervalo:

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

parâmetros do algoritmo de Metropolis:

valor inicial de \(x\)

número de iterações desejadas

algoritmo mcmc:

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)
}

Esta função devolve uma cadeia com niter+1 valores de x.

Exemplo:

#inicialização
set.seed(123)

# vou começar longe da média
xini=0.1
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')
hist(y,col='green',main='distribuição de x',xlab='x',ylab='frequência')

# vamos adicionar uma reta horizontal no valor esperado para a média de P(x)
# mu = alfa/(alfa+beta)
alfa = 6
beta = 3
mu = alfa/(alfa+beta)
abline(v = mu, 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 100 iterações iniciais de 4 cadeias com condições iniciais diferentes:

#inicialização
set.seed(321)

ncadeia=4
niter=100
xinic=c(0.1,0.5,0.7,0.99)

#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(0,1))
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 = mu, col="blue" )

Nesse exemplo, vemos porque as primeiras 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=0.1
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 = mu, col="blue", lwd = 3)
hist(y,col='green',main='distribuição de x',xlab='x',ylab='frequência')
abline(v = mu, col="blue", lwd = 3)

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

burnin = 1000

# taxa de aceitação: fração dos valores propostos que são aceitos
aceitacao = 1-mean(duplicated(y[-(1:burnin)]))
aceitacao

[1] 0.8214643

# média, mediana, desvio padrão
mean(y[-(1:burnin)])

[1] 0.6739694

median(y[-(1:burnin)])

[1] 0.6859825

sd(y[-(1:burnin)])

[1] 0.1391998

# média teórica
mu

[1] 0.6666667

# intervalo de confiança de 95%
quantile(y[-(1:burnin)], probs = c(0.025,0.975))

     2.5%     97.5% 
0.3689704 0.9020958

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] 0.6644368 0.6867003 0.6828202 0.6620090

#desvio padrão de cada segmento da cadeia
sd.seg

[1] 0.1367991 0.1377622 0.1354616 0.1449225

#media das cadeias
mean(media.seg)

[1] 0.6739916

#desvio padrão das cadeias
sd(media.seg)

[1] 0.01257421

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(0.1,0.2,0.3,0.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 = mu, col="yellow",lwd=3)

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

# média de cada cadeia:
media=rep(0,ncadeia)
for (i in 1:ncadeia){media[i]=mean(Y[-(1:burnin),i])}
media

[1] 0.6610813 0.6582730 0.6587578 0.6591039

# desvio padrão de cada cadeia:
dp=rep(0,ncadeia)
for (i in 1:ncadeia){dp[i]=sd(Y[-(1:burnin),i])}
dp

[1] 0.1475402 0.1447327 0.1470697 0.1507975

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] 0.659304

sd(Ysuper)

[1] 0.1475492

# 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% 
0.3442906 0.9070596

5. 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)
z = z+(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 2.782054 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

6. exercício da semana

inicialize o gerador de números aleatórios com seu número USP.

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.