Sampling Importance Resampling

Mais uns exemplos

1 Sampling importance resampling (SIR)

Em geral, o problema computacional da inferência Bayesiana é conseguir amostrar da posteriori \(p(\theta|dados)\), que é obtida pela combinação de uma verossimilhanca \(L(\theta|dados)\) com uma priori \(p(\theta)\), i.e. \[ p(\theta|dados) = \frac1kp(\theta)L(\theta|dados). \] Como sabemos avaliar \(p(\theta)\) e \(L(\theta|dados)\) pontualmente, sabemos também avaliar \(p(\theta|dados)\) pontualmente, a menos de uma constante (\(1/k\) aqui!).

O que o algoritmo SIR faz é amostrar de uma distribuição auxiliar (ou proposta ou candidata) \(q(\theta)\) que i) seja fácil de amsotrar; ii) possa ser avaliada pontualmente e iii) se “aproxime” da posterior \(p(\theta|dados)\).

Daí o SIR funciona assim:

1. Uma amostra \(\{\theta^{(1)},\ldots,\theta^{(M)}\}\) é obtida da densidade proposta \(q(\theta)\). O tamanho da amostra Monte Carlo pode ser \(M=1.000\) ou \(M=10.000\) ou ainda \(M=100.000\), ou qualquer outra quantidade dependendo da complexidade contextual.

2. Cada componente \(\theta^{(i)}\) da amostra recebe um peso \[ w^{(i)} = \frac{p(\theta^{(i)}|dados)}{q(\theta^{(i)})}, \] ou seja, \[ w^{(i)} = \frac1k \frac{p(\theta^{(i)}L(\theta^{(i)}|dados)}{q(\theta^{(i)})}. \] A constante \(k\) é irrelevante pois os pesos podem ser recalculados dividindo-os pela sua soma; daí o \(k\) é cancelado no numerador e no denominador.

3. Reamostra-se \(N \leq M\) elementos do conjunto \(\{\theta^{(1)},\ldots,\theta^{(M)}\}\) com pesos \(\{w^{(1)},\ldots,w^{(M)}\}\) e com reposição. Assim, obtem-se uma amostra da posterior \(p(\theta|dados)\):

\[ \left\{\theta_*^{(1)},\ldots,\theta_*^{(N)}\right\}. \]

1.1 Caso particular: Utilização da priori como proposta

Nesse caso \(q(\theta) = p(\theta)\) e os pesos ficam \[ w^{(i)} = \frac1k \frac{p(\theta^{(i)})L(\theta^{(i)}|dados)}{p(\theta^{(i)})} = \frac1k L(\theta^{(i)}|dados). \] Portanto, os pesos do SIR só serão a verossimilhança se a densidade proposta for a própria priori.

Observação: Apesar da convenienca de se amostrar da priori e calcular pesos só com a verossimilhanca, o grande problema reside no fato que em geral priori e verossimilhança podem estar em regiões distintas e portanto amostrar da priori pode não ser uma boa estratégia. Lembre-se que eu disse acima que desejamos uma proposta que se “aproxime” da posterior \(p(\theta|dados)\).

2 Exemplos

2.1 Proposta envelopa a distribuição de interesse

Distribuição de interesse: \(h(\theta) = \frac1k exp(-0.5\theta^2)\)

Aqui \(k\) é a constante normalizadora.

Distribuição proposta: \(N(0,5^2)\)

M     = 10000
N     = 10000
theta = rnorm(M,0,5)
w     = exp(-0.5*theta^2)/dnorm(theta,0,5)
w     = w/sum(w) 
thetastar = sample(theta,size=N,replace=TRUE,prob=w)

par(mfrow=c(1,2))
thetas = seq(-5,5,length=1000)
hist(theta,prob=TRUE,main="Amostras de q(theta)",ylim=c(0,0.4))
lines(thetas,dnorm(thetas),lwd=2,col=2)
hist(thetastar,xlab="theta",prob=TRUE,main="Amostras de h(theta)")
lines(thetas,dnorm(thetas),lwd=2,col=2)

2.2 Proposta NÃO envelopa a distribuição de interesse

Distribuição de interesse: \(h(\theta) = \frac1k exp(-0.5\theta^2)\)

Novamente, \(k\) é a constante normalizadora.

Distribuição proposta: \(N(-15,5^2)\)

M     = 10000
N     = 10000
theta = rnorm(M,-15,5)
w = exp(-0.5*theta^2)/dnorm(theta,-15,5)
w = w/sum(w) 
thetastar = sample(theta,size=N,replace=TRUE,prob=w)

par(mfrow=c(1,2))
thetas = seq(-5,5,length=1000)
hist(theta,prob=TRUE,main="Amostras de q(theta)",ylim=c(0,0.4))
lines(thetas,dnorm(thetas),lwd=2,col=2)
hist(thetastar,xlab="theta",prob=TRUE,main="Amostras de h(theta)")
lines(thetas,dnorm(thetas),lwd=2,col=2)

2.3 SIR na inferência Bayesiana

Dados = \(\{1,0,1,1,0,1,0,0,1,1\}\)

Modelo: \(x_i|\theta \ iid \ Bernoulli(\theta)\) for \(i=1,\ldots,n\) and \(\theta \in (0,1)\).

Priori: \(\theta \sim Beta(6,2)\), such that \(E(\theta)=0.75\) and \(V(\theta)=0.144^2\).

Posteriori (distribuição de interesse) \[ p(\theta|dados) = \frac1k \theta^{6-1}(1-\theta)^{2-1} \theta^6(1-\theta)^4 \] Novamente, \(k\) é a constante de normalização.

Distribuição proposta: \(Beta(2,2)\)

Você pode trocar essa proposta por qualquer distribuição no intervalo \((0,1)\).

M     = 10000
N     = 10000
theta = rbeta(M,2,2)
w = dbeta(theta,6,2)*theta^6*(1-theta)^4/dbeta(theta,2,2)
w = w/sum(w) 
thetastar = sample(theta,size=N,replace=TRUE,prob=w)

par(mfrow=c(1,2))
thetas = seq(0,1,length=1000)
hist(theta,prob=TRUE,main="Amostras da proposta",ylim=c(0,4),xlim=c(0,1))
lines(thetas,dbeta(thetas,12,6),lwd=2,col=2)
hist(thetastar,xlab="theta",prob=TRUE,main="Amostras da posteriori",ylim=c(0,4),xlim=c(0,1))
lines(thetas,dbeta(thetas,12,6),lwd=2,col=2)

c(mean(thetastar), sqrt(var(thetastar)))

## [1] 0.6659571 0.1076641

quantile(thetastar,c(0.025,0.5,0.975))

##      2.5%       50%     97.5% 
## 0.4426397 0.6704314 0.8562359

2.4 Vamos repetir 2.3) mas usando a priori como proposta

Distribuição proposta = Priori: \(\theta \sim Beta(6,2)\).

M     = 10000
N     = 10000
theta = rbeta(M,6,2)
w = theta^6*(1-theta)^4
w = w/sum(w) 
thetastar = sample(theta,size=N,replace=TRUE,prob=w)

par(mfrow=c(1,2))
thetas = seq(0,1,length=1000)
hist(theta,prob=TRUE,main="Amostras da proposta",ylim=c(0,4),xlim=c(0,1))
lines(thetas,dbeta(thetas,12,6),lwd=2,col=2)
hist(thetastar,xlab="theta",prob=TRUE,main="Amostras da posteriori",,ylim=c(0,4),xlim=c(0,1))
lines(thetas,dbeta(thetas,12,6),lwd=2,col=2)

c(mean(thetastar), sqrt(var(thetastar)))

## [1] 0.667744 0.107167

quantile(thetastar,c(0.025,0.5,0.975))

##      2.5%       50%     97.5% 
## 0.4460907 0.6722582 0.8584376

2.5 Vamos repetir 2.4) mas agora com MUITO mais dados e priori MUITO longe da verossimilhança

Dados = \(\{1,0,1,1,0,1,0,0,1,1,1,0,1,1,0,1,0,0,1,1,1,0,1,1,0,1,0,0,1,1\}\)

Priori: \(\theta \sim Beta(3,12)\)

Verossimilhança: muito mais concentrada! \[ L(\theta|dados) = \theta^{18}(1-\theta)^{12} \]

M     = 10000
N     = 10000
theta = rbeta(M,3,12)
w = theta^(18)*(1-theta)^(12)
w = w/sum(w) 
thetastar = sample(theta,size=N,replace=TRUE,prob=w)

par(mfrow=c(1,2))
thetas = seq(0,1,length=1000)
hist(theta,prob=TRUE,main="Amostras da proposta",ylim=c(0,6),xlim=c(0,1))
lines(thetas,dbeta(thetas,21,24),lwd=2,col=2)
hist(thetastar,xlab="theta",prob=TRUE,main="Amostras da posteriori",ylim=c(0,6),xlim=c(0,1))
lines(thetas,dbeta(thetas,21,24),lwd=2,col=2)

c(mean(thetastar), sqrt(var(thetastar)))

## [1] 0.45943438 0.07241472

quantile(thetastar,c(0.025,0.5,0.975))

##      2.5%       50%     97.5% 
## 0.3217701 0.4584958 0.6163900