3 Integração Monte Carlo

Enquanto o Capítulo 2 enfocou as técnicas de simulação úteis para produzir variáveis aleatórias por computador, este capítulo apresenta os principais conceitos dos métodos de Monte Carlo; isto é, aproveitando a disponibilidade de variáveis aleatórias geradas por computador para aproximar integrais univariados e multidimensionais.

Na Seção 3.2, introduzimos a noção básica de aproximações de Monte Carlo como um subproduto da Lei dos Grandes Números, enquanto a Seção 3.3 destaca a universalidade da abordagem enfatizando a versatilidade da representação de uma integral como uma esperança. O Capítulo 5 também tratará da resolução de problemas de otimização por técnicas de simulação.



3.1 Introdução

Duas classes principais de problemas numéricos que surgem na inferência estatística são problemas de otimização e problemas de integração. Com efeito, numerosos exemplos (ver Rubinstein, 1981, Gentle, 2002, ou Robert, 2001) mostram que nem sempre é possível calcular analiticamente os estimadores associados a um dado paradigma (máxima verossimilhança, Bayes, método dos momentos, etc.).

Assim, seja qual for o tipo de inferência estatística, muitas vezes somos levados a considerar soluções numéricas. O capítulo anterior introduziu uma série de métodos para a geração computacional de variáveis aleatórias com qualquer distribuição dada e, portanto, fornece uma base para a construção de soluções para nossos problemas estatísticos. Uma solução geral é, de fato, usar a simulação, das distribuições verdadeiras ou de algumas substitutas, para calcular as quantidades de interesse. Na configuração da teoria da decisão, seja ela clássica ou Bayesiana, esta solução é natural uma vez que os riscos e os estimadores de Bayes envolvem integrais em relação às distribuições de probabilidade.

Observe que a possibilidade de produzir um número quase infinito de variáveis aleatórias distribuídas de acordo com uma determinada distribuição nos dá acesso ao uso de resultados freqüentistas e assintóticos com muito mais facilidade do que nas configurações inferenciais usuais, onde o tamanho da amostra é mais frequentemente fixo. Pode-se, portanto, aplicar resultados probabilísticos como a Lei dos Grandes Números ou o Teorema do Limite Central, pois permitem avaliar a convergência dos métodos de simulação o que equivale aos limites determinísticos utilizados pelas abordagens numéricas.

Antes de iniciar a descrição das técnicas de Monte Carlo, observe que uma alternativa aparentemente óbvia ao uso de métodos de simulação para aproximar integrais da forma \[ \int_\mathcal{X} h(x) f(x)\mbox{d}x, \] onde \(f\) é uma densidade de probabilidade, seria confiar em métodos numéricos como o de Simpson e as regras do trapézio. Por exemplo, R oferece duas funções relacionadas que executam integração unidimensional, area, na biblioteca MASS e integrate. No entanto, area não pode lidar com limites infinitos na integral e, portanto, requer algum conhecimento prévio da região de integração. A outra função, integrate, aceita limites infinitos, mas infelizmente é muito frágil e pode produzir uma saída não confiável.


Exemplo 3.1.

Como teste, comparamos o uso de integrate na integral \[ \int_0^\infty x^{\lambda-1}\exp(-x)\mbox{d}x, \] com o cálculo de \(\Gamma(\lambda)\) via função gamma. Implementando esta comparação como: 

ch = function(la){ integrate(function(x){x^(la-1)*exp(-x)},0,Inf)$val}
plot(lgamma(seq(.01,10,le=100)),
     log(apply(as.matrix(seq(.01,10,le=100)),1,ch)),xlab="log(integrate(f))",
     ylab=expression(log(Gamma(lambda))),pch=19,cex=.6)
grid()

Figura. 3.1. Avaliação de integrate da integral \(\Gamma(\lambda)\), com seu valor verdadeiro.

Obtemos a sequência representada na Figura 3.1, que não apresenta discrepância mesmo para valores muito pequenos de \(\lambda\).


Uma dificuldade principal com os métodos de integração numérica, como integrate, é que eles muitas vezes falham em identificar a região de importância para a função ser integrada. Em contraste, os métodos de simulação visam naturalmente essa região explorando as informações fornecidas pela densidade de probabilidade associada às integrais.


Exemplo 3.2.

Considere uma amostra de dez variáveis aleatórias \(X_i\sim Cauchy(\theta)\), \(1\leq i\leq 10\) com parâmetro de localização \(\theta = 350\). A função pseudo-marginal amostral, sob um plano anterior é então \[ m(x)=\int_{-\infty}^\infty \prod_{i=1}^{10} \dfrac{1}{\pi}\dfrac{1}{1+(x_i-\theta)^2}\mbox{d}\theta\cdot \]

No entanto, integrate retorna um valor numérico errado.

cac = rcauchy(10)+350
lik = function(the){ 
  u = dcauchy(cac[1]-the)
  for (i in 2:10)
    u = u*dcauchy(cac[i]-the)
  return(u)}
integrate(lik,-Inf,Inf)
## 5.164922e-44 with absolute error < 1e-43
integrate(lik,200,400)
## 6.54421e-08 with absolute error < 1.2e-07

e deixa de sinalizar a dificuldade já que a avaliação do erro é absurdamente pequena. Além disso, o resultado não é comparável à área:

library(MASS)
cac=rcauchy(10)
nin=function(a){integrate(lik,-a,a)$val}
nan=function(a){area(lik,-a,a)}
x=seq(1,10^3,le=10^4)
y=log(apply(as.matrix(x),1,nin))
z=log(apply(as.matrix(x),1,nan))
plot(x,y,type="l",ylim=range(cbind(y,z)),lwd=2)
lines(x,z,lty=2,col="sienna",lwd=2)
grid()
Figura. 3.2. Comparação dos resultados com as funções integrate e area na integral de uma verossimilhança da função de densidade Cauchy em escala logarítmica, o resultado da área corresponde à curva tracejada acima.

O uso da função area nesse caso produz uma avaliação mais confiável, conforme mostrado na Figura 3.2, uma vez que area(lik,-a,a) se achata à medida que aumenta, mas isso obviamente requer algum conhecimento prévio sobre a localização do modo do integrando.


Por último, as ferramentas de integração numérica não podem enfrentar facilmente as integrais altamente ou moderadamente multidimensionais que são a regra em problemas estatísticos. Desenvolver ferramentas de integração específicas para esses problemas seria muito caro, especialmente porque podemos tirar proveito da natureza probabilística dessas integrais.


3.2 Integração clássica de Monte Carlo

Antes de aplicar nossas técnicas de simulação a problemas práticos, vamos relembrar as propriedades que justificam seu uso, consultando Robert and Casella (2004) para mais detalhes. O problema genérico é sobre avaliar a integral \[\begin{equation} \tag{3.1} \mbox{E}_f\big(f(X)\big)=\int_\mathcal{X} h(x)f(x)\mbox{d}x \end{equation}\] onde \(\mathcal{X}\) denota o conjunto onde a variável aleatória \(X\) assume seus valores, que normalmente é igual ao suporte da densidade \(f\).

O princípio do método de Monte Carlo para aproximar (3.1) é gerar uma amostra \((X_1,\cdots, X_n)\) a partir da densidade \(f\) e propor como aproximação a média empírica \[ \overline{h}_n=\dfrac{1}{n}\sum_{i=1}^n h(x_i), \] calculado por mean(h(x)) em R, uma vez que \(\overline{h}_n\) converge quase certamente, ou seja, converge para quase todas as sequências geradas para \(\mbox{E}_f\big(f(X)\big)\), pela Lei Forte dos Grandes Números.

Além disso, quando \(h^2(X)\) tem esperança finita sob \(f\), a velocidade de convergência de \(\overline{h}_n\) pode ser avaliada, pois a convergência ocorre a uma velocidade \(O(\sqrt{n})\) e a variância assintótica da aproximação é \[ \mbox{Var}\big(\overline{h}_n\big)=\dfrac{1}{n}\int_\mathcal{X} \left( h(x)-\mbox{E}_f\big(f(X)\big)\right) f(x)\mbox{d}x, \] que também pode ser estimada a partir da amostra \((X_1,\cdots, X_n)\) através de \[ v_n =\dfrac{1}{n^2}\sum_{i=1}^n \big( h(x_i)-\overline{h}_n\big)^2 \cdot \]

Mais especificamente, devido ao Teorema do Limite Central, para \(n\) grande, \[ \dfrac{\overline{h}_n-\mbox{E}_f\big(f(X)\big)}{\sqrt{v_n}}, \] é aproximadamente distribuído como uma variável \(N(0,1)\), e isso leva à construção de um teste de convergência e limites de confiança aproximados de \(\mbox{E}_f\big(f(X)\big)\).


Exemplo 3.3

Para a função, que tem integral analítica, \[\begin{equation} \tag{3.2} h(x) = \big(\cos(50x) + \sin(20x)\big)^2, \end{equation}\] representada no painel superior da Figura 3.3, considere calcular sua integral sobre [0,1].

Pode ser vista como uma esperança uniforme e, portanto, em R geramos variáveis aleatórias \(U_1, U_2,\cdots,U_n\) independentes identicamente distribuídas \(U(0,1)\) e aproximadas 

set.seed(174)
h = function(x){(cos(50*x)+sin(20*x))^2}
par(mar=c(4,4,2,1),mfrow=c(2,1))
curve(h,xlab="Function",ylab="",lwd=2)
grid()
integrate(h,0,1)
## 0.9652009 with absolute error < 1.9e-10
x = h(runif(10^4))
estint = cumsum(x)/(1:10^4)
esterr = sqrt(cumsum((x-estint)^2))/(1:10^4)
plot(estint, xlab="Mean and error range",type="l",lwd=2,
     ylim=mean(x)+20*c(-esterr[10^4],esterr[10^4]),ylab="")
lines(estint+2*esterr,col="gold",lwd=2)
lines(estint-2*esterr,col="gold",lwd=2)
abline(h=0.9652009,col="red",lwd=2,lty=1)
grid()

Figura. 3.3. Aproximação da integral da função (3.2): (painel superior) função (3.2) e (painel inferior) média \(\pm\) dois erros padrão contra iterações para uma única sequência de simulações.

Observe que a faixa de confiança produzida nesta figura não é uma faixa de confiança de 95% no sentido clássico, ou seja, não corresponde a uma faixa de confiança no gráfico de estimativas, mas sim à avaliação de confiança que você pode produzir para cada número de iterações, se você parar neste número de iterações.


Embora o bônus trazido pela avaliação simultânea do erro da estimativa Monte Carlo não possa ser contestado, você deve estar ciente de que ele só é confiável na medida em que \(v_n\) é uma estimativa adequada da variância de \(\overline{h}_n\). Em situações críticas, em que \(v_n\) não converge ou nem mesmo converge rápido o suficiente para que o Teorema do Limite Central seja aplicado, essa estimativa e a região de confiança associada a ela não são confiáveis.

Ao monitorar a convergência Monte Carlo, uma questão que será totalmente abordada no próximo capítulo, o comando R cumsum é bastante útil, pois calcula todas as somas parciais de uma sequência de uma só vez e, assim, permite a representação imediata da sequência de estimadores.

A metodologia Monte Carlo ilustrada pelo exemplo acima pode ser implementada com sucesso em uma ampla gama de casos onde as distribuições envolvidas no modelo podem ser simuladas.

Por exemplo, poderíamos usar as somas Monte Carlo para calcular uma função de distribuição normal, embora esta distribuição agora possa ser encontrada em todos os softwares e em muitas calculadoras de bolso.


Exemplo 3.4

Dada uma amostra normal \(N(0,1)\) de tamanho \(n\), \(x_1,\cdots,x_n\), a aproximação de \[ \Phi(t)=\int_{-\infty}^t \dfrac{1}{\sqrt{2\pi}}e^{-z^2/2}\mbox{d}z, \] pelo método Monte Carlo é \[ \widehat{\Phi}(t)=\dfrac{1}{n}\sum_{i=1}^n \pmb{1}_{\{x_i\leq t\}}, \] com variância exata \[ \mbox{Var}\big(\widehat{\Phi}(t)\big)(t) = \Phi(t)\big([1-\Phi(t)\big)/n, \] uma vez que as variáveis \(\pmb{1}_{\{x_i\leq t\}}\) são Bernoulli independentes com probabilidade de sucesso \(\Phi(t)\). 

x=rnorm(10^8)                     # whole sample
bound=qnorm(c(.5,.75,.8,.9,.95,.99,.999,.9999))
res=matrix(0,ncol=8,nrow=7)
for (i in 2:8)                    # lengthy loop!!
 for (j in 1:8)
 res[i-1,j]=mean(x[1:10^i]<bound[j])
matrix(as.numeric(format(res,digi=4)),ncol=8)
##        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]
## [1,] 0.5900 0.8000 0.8400 0.8800 0.9300 0.9800 1.0000 1.0000
## [2,] 0.5150 0.7630 0.8110 0.9130 0.9540 0.9940 0.9990 1.0000
## [3,] 0.5002 0.7504 0.7995 0.9017 0.9517 0.9900 0.9989 0.9997
## [4,] 0.5000 0.7475 0.7962 0.8974 0.9490 0.9898 0.9990 0.9999
## [5,] 0.4996 0.7496 0.7996 0.8997 0.9503 0.9901 0.9990 0.9999
## [6,] 0.4997 0.7498 0.7999 0.8999 0.9500 0.9900 0.9990 0.9999
## [7,] 0.5000 0.7500 0.8000 0.9000 0.9500 0.9900 0.9990 0.9999

Tabela 3.1. Avaliação de algumas probabilidades normais Pr(X t) por um experimento regular de Monte Carlo baseado em n replicações de uma geração normal. A última linha alcança os valores exatos.

Para valores de \(t\) em torno de \(t = 0\), a variância é, portanto, aproximadamente \(1/4n\) e, para obter uma precisão de quatro casas decimais, precisamos de \(2\times \sqrt{1/4n}\leq 10^{-4}\) simulações, ou seja, cerca de \(n = (10^4)^2 = 10^8\) simulações.

A Tabela 3.1 apresenta a evolução dessa aproximação para vários valores de \(t\) e mostra uma avaliação precisa para 100 milhões de iterações. Observe que maior precisão absoluta é alcançada nas caudas e que métodos de simulação muito mais eficientes poderiam ser usados.


Como você deve ter notado, as saídas em R são representadas com todos os dígitos disponíveis, como em

rnorm(1)
## [1] -2.430835
Embora isso seja lógico do ponto de vista informático, não é recomendável produzir todos esses dígitos em ambientes estatísticos e de simulação porque a maioria deles não é significativa e também porque prejudica a legibilidade da saída. A função format é bastante útil para reduzir o número de dígitos representados, conforme mostrado na última linha do programa R acima.

A aproximação Monte Carlo de uma função de distribuição ilustrada pelo Exemplo 3.4 tem aplicações não triviais, pois pode ser usada para avaliar a distribuição de uma estatística de teste, como um teste de razão de verossimilhança sob uma hipótese nula, conforme ilustrado em Robert and Casella (2004), bem como seu poder sob alternativas.

Neste estágio, pode parecer que a metodologia Monte Carlo introduzida nesta seção é suficiente para aproximar integrais como (3.1) de maneira controlada. No entanto, embora o método direto de Monte Carlo realmente forneça boas aproximações de (3.1) na maioria dos casos regulares, existem alternativas mais eficientes que não apenas evitam uma simulação direta de \(f\), mas também podem ser usadas repetidamente para várias integrais da forma (3.1).

O uso repetido pode ser tanto para uma família de funções \(h\) quanto para uma família de densidades \(f\). Além disso, problemas de simulação de cauda como no Exemplo 3.4 podem ser processados de forma muito mais eficiente do que a simulação de \(f\), pois simular eventos com uma probabilidade muito pequena requer um número muito grande de simulações sob \(f\) para atingir uma dada precisão relativa.


3.3 Amostragem de importância


O método que agora estudamos é chamado de amostragem por importância porque se baseia nas chamadas funções de importância, que são distribuições instrumentais, em vez das distribuições originais. De fato, uma avaliação de (3.1) baseada em simulações de \(f\) quase nunca é ótima no sentido de que o uso de distribuições alternativas pode melhorar a variância do estimador resultante de (3.1).


3.3.1 Uma mudança arbitrária de medida de referência

O método de amostragem por importância é baseado em uma representação alternativa de (3.1). Dada uma função de densidade arbitrária \(g\), que é estritamente positiva quando \(h\times f\) é diferente de zero, podemos realmente reescrever (3.1) como \[\begin{equation} \tag{3.3} \mbox{E}_f\big(h(X)\big)=\int_\mathcal{X} h(x)\dfrac{f(x)}{g(x)}g(x)\mbox{d}x = \mbox{E}_g\left(\dfrac{h(X)f(X)}{g(X)}\right), \end{equation}\] isto é, como uma esperança sob a densidade \(g\). Observe que \(\mathcal{X}\) é novamente o conjunto onde \(X\) assume seu valor e que pode, portanto, ser menor que o suporte da densidade \(g\).

Essa identidade de amostragem é de fundamental importância e justifica o uso do estimador \[\begin{equation} \tag{3.4} \dfrac{1}{n}\sum_{i=1}^n \dfrac{f(X_i)}{g(X_i)}h(X_i) \longrightarrow \mbox{E}_f\big(h(X)\big), \end{equation}\] com base em uma amostra \(X_1,\cdots, X_n\) gerada de \(g\), não de \(f\)!.

De fato, como (3.1) pode ser escrito como uma esperança sob \(g\), (3.4) converge para (3.1) pela mesma razão que o estimador regular Monte Carlo \(\overline{h}_n\) converge, qualquer que seja a escolha da distribuição \(g\), contanto que \(\mbox{supp} (g) \supset \mbox{supp} (h\times f)\).

Esta propriedade onipresente está relacionada ao fato de que (3.1) pode ser representada de infinitas maneiras por pares \((h,f)\) e, portanto, uma dada integral não está intrinsecamente associada a uma dada distribuição. Ao contrário, há uma liberdade quase absoluta em sua representação como esperança.

A restrição de \(\mbox{supp} (g) \supset \mbox{supp} (h\times f)\) é absoluta porque usar um suporte menor trunca a integral (3.3) e assim, produz um resultado enviesado. Isso significa, em particular, que ao considerar soluções não paramétricas para \(g\), o suporte do kernel deve ser irrestrito.


Exemplo 3.5

Conforme mencionado no final do Exemplo 3.4, a aproximação da probabilidades na cauda usando somas Monte Carlo padrão falham quando se avança o suficiente nas caudas.

Por exemplo, se \(Z\sim N (0,1)\) e estamos interessados na probabilidade \(P(Z > 4.5)\), que é muito pequena,

pnorm(-4.5,log=T)
## [1] -12.59242
simular \(Z^{(i)}\sim N(0,1)\) só produz um acerto uma vez em cerca de 3 milhões de iterações!

Obviamente, o problema é que agora estamos interessados na probabilidade de um evento muito raro e, portanto, a simulação ingênua de \(f\) exigirá um grande número de simulações para obter uma resposta estável. No entanto, graças à amostragem de importância, podemos melhorar muito nossa precisão e, assim, reduzir o número de simulações em várias ordens de grandeza.

Por exemplo, se considerarmos uma distribuição com suporte restrito a \((4.5,\infty)\), a variação adicional e desnecessária do estimador Monte Carlo devido à simulação de zeros, ou seja, quando \(x < 4.5\) desaparece.

Uma escolha natural é tomar \(g\) como a densidade da distribuição exponencial Exp(1) truncada em 4.5, \[ g(y)=\exp(-y) \Big/ \int_{4.5}^\infty \exp(-x)\mbox{d}x=\exp\big(-(y-4.5)\big) \] e o estimador amostral de importância correspondente da probabilidade na cauda é \[ \dfrac{1}{n} \sum_{i=1}^n \dfrac{f(Y^{(i)}}{g(Y^{(i)})} = \dfrac{1}{n} \sum_{i=1}^n \dfrac{\exp\big( −Y_i^2/2+Y_i −4.5\big)}{\sqrt{2\pi}} \] onde os \(Y_i\)’s são variáveis aleatórias independentes identicamente distribuídas geradas de \(g\).

O código correspondente é

set.seed(78666)
Nsim = 10^3
y = rexp(Nsim)+4.5
weit = dnorm(y)/dexp(y-4.5)
plot(cumsum(weit)/1:Nsim, type = "l")
abline(a = pnorm(-4.5), b = 0, col="red")
grid()

Figura. 3.4. Convergência da aproximação à probabilidade da cauda normal \(P(Z\geq 4.5)\) por amostragem de importância, com base em uma sequência simulada a partir de uma distribuição exponencial truncada. A linha reta corresponde ao verdadeiro valor da integral.

O valor final é então \[ 0.003130292 = 3.130292\times 10^{-3}, \] a ser comparado com o valor verdadeiro de \[ 3.397673\times 10^{-6}\cdot \] Conforme mostrado na Figura 3.4, a precisão da aproximação é notável, especialmente quando comparada com os requisitos de tamanho originais impostos por uma simulação normal.


A amostragem de importância é, portanto, de considerável interesse, uma vez que coloca poucas restrições na escolha da distribuição instrumental \(g\), que pode ser escolhida entre distribuições fáceis de simular ou eficientes na aproximação da integral. Além disso, a mesma amostra gerada a partir de \(g\) pode ser usada repetidamente, não apenas para diferentes funções \(h\), mas também para diferentes densidades \(f\).


Exemplo 3.6

Este exemplo decorre de uma configuração bayesiana: ao considerar uma observação \(X\) da distribuição \(Beta(\alpha,\beta)\), \[ X\sim \dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-1}(1-x)^{\beta-1} \pmb{1}_{(0,1)}(x), \] existe uma família de apriores conjugados em \((\alpha,\beta)\) da forma \[ \pi(\alpha,\beta) \approx \left(\dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^\lambda x_0^\alpha y_0^\beta, \] onde \(\lambda\), \(x_0\), \(y_0\) são hiperparâmetros, já que a distribuição a posteriori é então igual a \[ \pi(\alpha,\beta \, | \, x) \approx \left(\dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^{\lambda+1} \big(x x_0\big)^\alpha \big((1-x) y_0\big)^\beta\cdot \]

Essa família de distribuições é intratável apenas por causa da dificuldade de lidar com funções gama. Simular diretamente de \(\pi(\alpha,\beta | x)\) é, portanto, impossível.

Precisamos, portanto, usar uma distribuição substituta \(g(\alpha,\beta)\) e podemos ter uma ideia preliminar observando uma representação de imagem de \(\pi(\alpha,\beta | x)\). Se tomarmos \(\lambda = 1\), \(x_0 = y_0 = 0.5\) e \(x = 0.6\), o código R é: 

par(mfrow=c(1,2))
f = function(a,b){
  exp(2*(lgamma(a+b)-lgamma(a)-lgamma(b))+a*log(.3)+b*log(.2))}
aa = 1:150
#alpha grid for image
bb = 1:100
#beta grid for image
post = outer(aa,bb,f)
image(aa,bb,post,xlab=expression(alpha),ylab=" ")
contour(aa,bb,post,add=TRUE)
grid()
x = matrix(rt(2*10^4,3),ncol=2)
#T sample
E = matrix(c(220,190,190,180),ncol=2) #Scale matrix
image(aa,bb,post,xlab=expression(alpha),ylab=" ")
y = t(t(chol(E))%*%t(x)+c(50,45))
points(y,cex=.6,pch=19)
grid()

Figura. 3.5. (esquerda) Representação da distribuição a posteriori \(\pi(\alpha,\beta | x)\) nos parâmetros de uma distribuição \(B(\alpha,\beta)\) para \(x\) = 0.6. (à direita) Superposição de uma amostra de \(10^3\) pontos de uma distribuição \(t\)-Student \(T(3,\mu,\Sigma)\) usada como uma função de importância.

O comando outer é uma abreviação útil para calcular uma matriz

A = outer(aa, bb, FUN="f")

de dimensão c(dim(a),dim(b)) cujo elemento A[i,j] é igual a f(a[i],b[j]). Embora seja muito mais rápido que o laço (loop) de alocação dupla básico:

system.time(outer(aa,bb,f))
##   usuário   sistema decorrido 
##      0.02      0.00      0.02
system.time(for (j in 1:100){for (i in 1:150) post[i,j]=f(a=aa[i],b=bb[j])})
##   usuário   sistema decorrido 
##      0.04      0.00      0.04

compara rapidamente com um único loop de alocação

system.time(outer(aa,bb,f))
##   usuário   sistema decorrido 
##         0         0         0
system.time(for (j in 1:100){post[,j]=f(a=aa,b=bb[j])})
##   usuário   sistema decorrido 
##      0.02      0.00      0.02
system.time(for (i in 1:150){post[i,]=f(a=aa[i],b=bb)})
##   usuário   sistema decorrido 
##      0.02      0.00      0.01
e, portanto, não oferece uma maneira supereficiente de alocar valores a uma matriz.

O exame da Figura 3.5 (esquerda) mostra que uma distribuição normal ou \(t\)-Student no par \((\alpha,\beta)\) pode ser apropriada. Escolhendo uma distribuição \(T(3,\mu,\Sigma)\) (\(t\)-Student) com \(\mu = (50, 45)\) e \[ \Sigma = \begin{pmatrix} 220 & 190 \\ 190 & 180 \end{pmatrix} \] produz um ajuste razoável, conforme mostrado na Figura 3.5 (direita) usando a superposição da simulação desta distribuição \(T(3,\mu,\Sigma)\) com a superfície da distribuição a posteriori. A matriz de covariância acima foi obtida por tentativa e erro, modificando as entradas até que a amostra na Figura 3.5 (à direita) se ajustasse bem.

Observe o uso de t(chol(E)) para garantir que a matriz de covariância seja E até um fator de 3, devido ao uso da distribuição \(t_3\)-Student.

Se a quantidade de interesse é a verossimilhança marginal, como na comparação do modelo bayesiano (Robert, 2001), \[ m(x) \, = \, \displaystyle \int_{\mathbb{R}^2_+} f(x|\alpha,\beta) \pi(\alpha,\beta) \mbox{d}\alpha \, \mbox{d}\beta \, = \, \displaystyle \dfrac{\displaystyle \int_{\mathbb{R}^2_+} \left(\dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^{\lambda+1} \big(x x_0\big)^\alpha \big((1-x) y_0\big)^\beta \mbox{d}\alpha \, \mbox{d}\beta}{x(1-x) \displaystyle \int_{\mathbb{R}^2_+} \left(\dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right)^{\lambda+1} x_0^\alpha y_0^\beta \mbox{d}\alpha \, \mbox{d}\beta} \] precisamos aproximar ambas as integrais e a mesma amostra \(t\)-Student pode ser usada para ambas, pois o ajuste é igualmente razoável na superfície a priori.

Esta aproximação \[\begin{equation} \tag{3.5} \begin{array}{rcl} \widehat{m}(x) & = & \displaystyle \sum_{i=1}^n \left( \dfrac{\Gamma(\alpha_i+\beta_i)}{\Gamma(\alpha_i)\Gamma(\beta_i)}\right)^{\lambda+1} \big(x x_0\big)^{\alpha_i} \big( (1-x)y_0\big)^{\beta_i}/ g(\alpha_i,\beta_i) \\ & & \qquad \qquad \Big/ \displaystyle \sum_{i=1}^n \left( \dfrac{\Gamma(\alpha_i+\beta_i)}{\Gamma(\alpha_i)\Gamma(\beta_i)}\right)^{\lambda+1} x_0^{\alpha_i} y_0^{\beta_i} / g(\alpha_i,\beta_i) \end{array} \end{equation}\] onde \((\alpha_i,\beta_i)_{1\leq i\leq n}\) são realizações iid de \(g\).

Mostramos a seguir a implementação no R. 

ine = apply(y,1,min)
y = y[ine>0,]
x = x[ine>0,]
normx = sqrt(x[,1]^2+x[,2]^2)
f = function(a) exp(2*(lgamma(a[,1]+a[,2])-lgamma(a[,1]) - lgamma(a[,2]))+a[,1]*log(.3)+a[,2]*log(.2))
h = function(a) exp(1*(lgamma(a[,1]+a[,2])-lgamma(a[,1]) - lgamma(a[,2]))+a[,1]*log(.5)+a[,2]*log(.5))
den = dt(normx,3)
mean(f(y)/den)/mean(h(y)/den) 
## [1] 0.2584168

Nossa aproximação da verossimilhança marginal, com base nessas simulações é, portanto, 0.2584168. Da mesma forma, as esperanças a posteriori dos parâmetros \(\alpha\) e \(\beta\) são obtidas por

mean(y[,1]*f(y)/den)/mean(h(y)/den)
## [1] 24.0176
mean(y[,2]*f(y)/den)/mean(h(y)/den)
## [1] 20.54095
ou seja, são aproximadamente iguais a 24.0176 e 20.54095, respectivamente.



3.3.2 Reamostragem da importância da amostragem

A técnica de amostragem de importância faz mais do que aproximar integrais, pois fornece uma maneira alternativa de simular distribuições complexas. Lembre-se de que o método produz uma amostra \(X_1,\cdots,X_n\) simulada de \(g\), junto com seus pesos de importância \(f(X_i)/g(X_i)\). Esta amostra pode então ser reciclada por reamostragem multinomial em uma amostra que é, quase, de \(f\).

De fato, se pudéssemos amostrar com reposição da população ponderada \(\{X_1,\cdots,X_n\}\), escolhendo \(X_i\) com probabilidade \(f(X_i)/ng(X_i)\), teríamos um resultado \(X^∗\) com distribuíção \[ P\big(X^* \in A \big) \, = \, \displaystyle \sum_{i=1}^n P\Big(X^* \in A \bigcap X^* = X_i\Big) \, = \, \int_A \dfrac{f(x)}{g(x)}g(x)\mbox{d}x \, = \, \int_A f(x) \mbox{d} x, \] e o método então produziria uma simulação exata de \(f\)! Infelizmente, as probabilidades \(f(X_i)/ng(X_i)\) não somam 1, pior, algumas podem até ser maiores que 1 e precisam ser renormalizadas em \((i = 1,\cdots, n)\) \[\begin{equation} \tag{3.6} \omega_i = \dfrac{1}{n} \left( f(X_i)/g(X_i)\right) \Big/ \dfrac{1}{n}\sum_{i=1}^n \left( f(X_i)/g(X_i)\right), \end{equation}\] que também pode ser usado em situações em que \(f\) ou \(g\) não possuem uma constante de normalização.

O denominador de (3.6) também está estimando a(s) constante(s) ausente(s). Este é, por exemplo, o caso do Exemplo 3.6: A constante de normalização ausente no a priori é estimada por mean(h(y)/den) no código acima.

Enquanto o denominador está convergindo quase certamente para um, a renormalização induz um viés na distribuição dos valores reamostrados. No entanto, para tamanhos de amostra grandes, esse viés é desprezível e, portanto, usaremos a reamostragem multinomial ou uma versão melhorada; consulte os Exercícios 3.6 e 3.12, para aproximar as amostras geradas a partir de \(f\).

Os pesos de importância fornecem apenas uma avaliação relativa da adequação da amostra simulada à densidade alvo, pois indicam o quanto é mais provável que \(X_i\) seja simulado a partir de \(f\) em comparação com \(X_j\), mas não devem ser superinterpretados. Por exemplo, se \(X_i\) tem um peso auto-normalizado próximo de 1, isso não significa que esse valor seja muito provável de ser gerado a partir de \(f\), mas simplesmente que é muito mais provável do que os outros valores simulados! Mesmo quando o ajuste entre \(f\) e \(g\) é muito ruim, essa ocorrência está fadada a acontecer. Portanto, mais confiável indicadores devem ser usados para julgar a adequação de \(g\) contra \(f\).


Exemplo 3.7 (Continuação do Exemplo 3.6)

A validade da aproximação A estimativa (3.5) da verossimilhança marginal, (ou seja, a convergência da solução de amostragem de importância) pode ser avaliada por meios gráficos como segue: 

par(mfrow=c(2,2),mar=c(4,4,2,1))
weit = (f(y)/den)/mean(h(y)/den)
image(aa,bb,post,xlab=expression(alpha), ylab=expression(beta))
grid()
points(y[sample(1:length(weit),10^3,rep=TRUE,pro=weit),],cex=.6,pch=19)
boxplot(weit,ylab="importance weight")
grid()
plot(cumsum(weit)/(1:length(weit)),type="l", xlab="simulations", ylab="marginal likelihood")
grid()
boot = matrix(0,ncol=length(weit),nrow=100)
for (t in 1:100) boot[t,]=cumsum(sample(weit))/(1:length(weit))
uppa = apply(boot,2,quantile,.95)
lowa = apply(boot,2,quantile,.05)
polygon(c(1:length(weit),length(weit):1),c(uppa,rev(lowa)), col="gold")
lines(cumsum(weit)/(1:length(weit)),lwd=2)
plot(cumsum(weit)^2/cumsum(weit^2),type="l", xlab="simulations", ylab="Effective sample size",lwd=2)
grid()

Figura. 3.6. (canto superior esquerdo) Superposição de \(10^3\) pontos reamostrados sobre a distribuição a posteriori \(\pi(\alpha,\beta|x)\) nos parâmetros de uma distribuição \(Beta(\alpha,\beta)\) para \(x = 0.6\). (canto superior direito) Boxplot dos pesos de importância. (canto superior direito) Convergência da aproximação \(\widehat{m}̂(x\)) e renderização bootstrap de sua variabilidade. (canto superior direito) Evolução do tamanho efetivo da amostra.

Não discutiremos detalhadamente todos esses indicadores, pois alguns serão explicados no próximo capítulo. O gráfico superior esquerdo na Figura 3.6 mostra que a amostra ponderada usando o peso de importância \(\pi(\alpha_i,\beta_i | x)/g(\alpha_i,\beta_i)\) produz uma renderização justa de uma amostra de \(\pi(\alpha,\beta|x)\).

Os pontos reamostrados não degeneram em poucos pontos, mas cobrem, com alta densidade, o intervalo correto para a distribuição alvo (compare com o lado direito da Figura 3.5). O gráfico superior direito dá uma representação da distribuição dos pesos de importância. Embora existam simulações com pesos muito maiores do que as outras, a dispersão do peso não é tão extrema a ponto de significar uma degenerescência do método.

Por exemplo, o ponto reponderado mais alto representa apenas 1% de toda a amostra. O gráfico inferior esquerdo representa a convergência do estimador \(\widehat{m}(x)\) a medida que \(n\) aumenta. A faixa colorida ao redor da sequência é uma renderização bootstrap (Seção 1.5) da variabilidade desse estimador que imita a faixa de confiança representada na Figura 3.3 a baixo custo. A curva inferior direita representa a perda de eficiência no uso da amostragem de importância pelo tamanho efetivo da amostragem; consulte a Seção 4.4, \[ \left( \sum_{i=1}^n \pi(\alpha_1,\beta_1|x)/g(\alpha_i,\beta_i)\right)^2 \Big/ \sum_{i=1}^n \left(\pi(\alpha_1,\beta_1|x)/g(\alpha_i,\beta_i)\right)^2, \] que deve ser igual a \(n\), se os \((\alpha_i,\beta_i)\)’s forem gerados a partir da a posteriori. O gráfico atual mostra que a amostra produzida tem uma eficiência de cerca de 6%. Consideraremos mais este indicador na Seção 4.2.



3.3.3 Seleção da função de importância

A versatilidade da técnica de amostragem de importância é alta, mas a desvantagem dessa versatilidade é que uma má escolha da função de importância \(g\) pode produzir resultados muito ruins. Embora a escolha ótima de \(g\) seja mais um exercício teórico (ver Rubinstein, 1981 ou Robert and Casella, 2004) do que algo útil, uma questão de relevância primária é considerar a variância do estimador resultante (3.3) quando julgar a adequação da função de importância correspondente \(g\).

De fato, enquanto (3.4) converge quase certamente para (3.1), dado que a esperança (3.1) existe, a variância desse estimador é finita somente quando a esperança \[ \mbox{E}_g\left( h^2(X) \dfrac{f^2(X)}{g^2(X)}\right) \, = \, \mbox{E}_f\left( h^2(X) \dfrac{f(X)}{g(X)}\right) \, = \, \int_\mathcal{X} h^2(x) \dfrac{f^2(x)}{g(x)} \mbox{d} x < \infty \] é finita. Embora não proíba exatamente funções de importância com caudas mais leves do que aquelas de \(f\), que levam a razões ilimitadas \(f/g\), essa condição enfatiza que essas funções têm muito mais probabilidade de levar a estimadores de variância infinita.

Antes de discutir esta questão com mais detalhes, vamos considerar um exemplo simples para ilustrar o impacto desastroso de um estimador de variância infinita.


Exemplo 3.8

Uma configuração simples em que ocorre variância infinita ao função de importância normal \(N(0,1)\) destinada a um alvo Cauchy \(C(0, 1)\). A razão \[ f(x)/g(x) \approx \exp(x^2/2)/(1 + x^2) \] é explosiva, pois mesmo valores moderadamente altos de \(x\) obtêm pesos de importância muito grandes.

Se você executar o código algumas vezes, você deve ver gráficos como o da Figura 3.7 ocorrendo, ou seja, padrões com grandes saltos na média acumulada, mesmo com um grande número de termos na média. 

par(mfrow=c(1,2),mar=c(4,4,2,1))
set.seed(76436)
x = rnorm(10^6)
wein = dcauchy(x)/dnorm(x)
boxplot(wein/sum(wein))
grid()
plot(cumsum(wein*(x>2)*(x<6))/cumsum(wein),type="l")
abline(a=pcauchy(6)-pcauchy(2),b=0,col="sienna")
grid()

Figura. 3.7. Evolução do estimador da amostragem de importância da probabilidade \(P(2\leq Z\leq 6)\) contra índices de iteração, quando \(Z\) tem distribuído Cauchy e a função de importância é Normal. A linha reta é o valor exato, 0.095.

Os saltos acontecem em valores das simulações para as quais \(\exp(x^2/2)/(1 + x^2)\) é grande, o que significa quando \(x\) é grande. A razão para esse fenômeno é que, como esses valores são raros na distribuição de importância normal, ou seja, mais raros do que no alvo Cauchy, eles precisam compensar sua raridade assumindo pesos altos.

Por exemplo, na Figura 3.7, o maior salto é devido a um valor de \(x = 5.49\) associado a um peso normalizado de \(ω_i = 0.094\). Isso significa que esse único ponto tem um peso de cerca de 10% em uma amostra de um milhão de pontos! Obviamente, é impossível confiar no resultado desta simulação, uma vez que o tamanho da amostra é irrelevante, ou seja, a maioria dos valores simulados tem um peso insignificante!.


Quando a razão \(f/g\) é ilimitada, os pesos de importância \(f(x_j)/g(x_j)\) frequentemente variam amplamente, dando muita importância a alguns valores \(x_j\) e assim degradando a eficiência do estimador (3.4). Como no exemplo acima, pode acontecer que a estimativa mude bruscamente de uma iteração para a próxima, mesmo após muitas iterações, devido a uma única simulação. Por outro lado, distribuições de importância \(g\) com caudas mais grossas que \(f\) garantem que o comportamento da razão \(f/g\) não seja a causa da divergência de \[ \mbox{E}_f\left( h^2(X) \dfrac{f(X)}{g(X)}\right)\cdot \]

O uso de distribuições de proposta de amostragem de importância de cauda mais grossa é quase uma obrigação ao considerar a aproximação de funções \(h\) tal que (3.1) existe, mas \(\mbox{E}\big(h^2(X)\big)\) não. Nesses casos, o uso de Monte Carlo regular não é possível, pois a média empírica dos \(h(X_i)\)’s não possui variância.

O estimador auto-normalizado (3.7) requer a mesma condição do caso não-normalizado para que a variância seja finita. Mas, conforme detalhado no Capítulo 4, a expressão da variância não está disponível de forma fechada e precisa ser aproximada pelos métodos de Monte Carlo.

Como recomendação genérica, sugerimos nesta fase procurar distribuições \(g\) para as quais \(|h|f/g\) seja quase constante ou pelo menos tenha um comportamento de cauda controlada, uma vez que é mais provável que produzam estimadores com variância finita.

Um requisito básico para funções \(h\) com suportes restritos como no Exemplo 3.5 é que \(g\) adote o mesmo suporte que \(h\), a menos que isso seja impedido pela complexidade de \(h\). Obviamente, isso requer o ajuste de uma nova função de importância para cada integrando \(h\) a ser considerado, mas esse é o preço a se pagar para se obter mais eficiência, como mostra o Exemplo 3.5.

Dado que a amostragem de importância se aplica principalmente em configurações onde \(f\) não é fácil de estudar, essa restrição nas caudas de \(f\) geralmente não é fácil de implementar, especialmente quando a dimensionalidade é alta. No entanto, existe uma solução genérica baseada na incorporação artificial de um componente de cauda gorda na função de importância \(g\). Essa solução é chamada de amostragem defensiva por Hesterberg (1995) e pode ser obtida substituindo a densidade \(g\) por uma densidade de mistura, \[\begin{equation} \tag{3.8} \rho g(x) + (1-\rho)\ell(x), \qquad 0<\rho<1, \end{equation}\] onde \(\rho\) é próximo de 1 e a densidade \(\ell\) é escolhida por suas caudas pesadas, por exemplo, uma distribuição Cauchy ou de Pareto, não necessariamente em conjunto com o problema em questão.

Assumindo que \(g\) é fornecido pela configuração, escolher a função de cauda pesada \(\ell\) é potencialmente delicado. No caso especial da inferência bayesiana, quando a distribuição alvo \(f\) é a distribuição a posteriori é, entretanto, natural escolher a priori se for o caso. De fato, essa função tem caudas mais pesadas do que \(f\) por construção e geralmente é uma distribuição padrão fácil de simular.

Usar a priori como a principal função de importância \(g\) não faria sentido por causa do desperdício induzido, supondo que os dados sejam informativos. Mas usá-lo como um fator estabilizador faz sentido.

Do ponto de vista operacional, gerar de (3.8) significa que as observações são geradas com probabilidade \(\rho\) de \(g\) e com probabilidade \(1-\rho\) de \(\ell\), usando um código como

mix = function(n=1,p=0.5){ 
  m=rbinom(1,size,pro=p) 
  c(simg(m),siml(n-m))}
se simg e siml denotam geradores de \(g\) e \(\ell\), respectivamente. Ressaltamos que o fato de alguns pontos serem gerados a partir de \(g\) e outros de \(\ell\) não impacta o peso de importância que é igual a \[ \dfrac{f(x)}{\rho g(x) + (1-\rho)\ell(x)} \] para todos os valores gerados.

Por construção, o estimador de amostragem de importância integra a variável uniforme usada para decidir entre \(g\) e \(\ell\). O condicionamento nessa variável uniforme induziria mais variabilidade e destruiria o propósito de usar uma mistura dividindo novamente por \(g(x)\) no peso de importância. Discutimos em detalhes essa perspectiva de marginalização no próximo capítulo, na Seção 4.6, onde as variáveis uniformes envolvidas na simulação são integradas no estimador.

Observe que a seleção de um número aleatório de simulações de \(g\) e \(\ell\) é desnecessária, no entanto, uma vez que gerar uma quantidae exatamente \(\rho n\) de variáveis \(x_i\)’s de \(g\) e uma quantidade \((1-\rho)n\), \(y_i\)’s de \(\ell\) produz um estimador não viciado (imparcial), sob a suposição de que \(\rho n\) seja um número inteiro, no sentido de que o estimador de amostragem de importância \[ \dfrac{1}{n} \sum_{i=1}^{\rho n} h(x_i)\dfrac{f(x_i)}{\rho g(x_i) + (1-\rho)\ell(x_i)}+\dfrac{1}{n} \sum_{i=1}^{(1-\rho) n} h(y_i)\dfrac{f(y_i)}{\rho g(y_i) + (1-\rho)\ell(y_i)} \] tem uma esperança global, se não temporal, igual a \(\mbox{E}_f\big(h(X)\big)\), ver Owen and Zhou (2000), para mais detalhes. Assim, simular um número fixo de pontos de cada distribuição é válido e interessante na medida em que elimina completamente a variabilidade devido à amostragem binomial acima.


Exemplo 3.9

Conforme indicado no Exercício 3.7, o cálculo da integral \[\begin{equation} \tag{3.9} \int_1^\infty \sqrt{\dfrac{x}{1-x}} t_2(x)\mbox{d}x \, = \, \int_1^\infty \sqrt{\dfrac{x}{1-x}} \dfrac{\Gamma(3/2)/\sqrt{2\pi}}{(1+x^2/2)^{3/2}} \mbox{d}x, \end{equation}\] é delicado porque a função \[ h(x) = \sqrt{1/(x-1)} \] não é integrável ao quadrado e, portanto, usar simulações da distribuição \(T_2\) produzirá uma variância infinita para o estimador Monte Carlo da integral.

Esta característica faz com que seja necessária uma mistura da densidade \(T_2\) com um \(\ell\) bem comportado. Alcançar a integrabilidade do quociente de funções \(h^2(x)f(x)/\ell(x)\) exige que \(\ell\) seja divergente em \(x = 1\) e que \(\ell\) diminua mais rápido que \(x^5\) quando \(x\) vai para o infinito.

Essas condições de contorno sugerem que: \[ \ell(x) \approx \dfrac{1}{\sqrt{x-1}}\dfrac{1}{x^{3/2}}\pmb{1}_{x>1}, \] que, definido até uma constante normalizadora, é uma densidade aceitável.

Para caracterizar essa densidade, você pode verificar que \[ \begin{array}{rcl} \displaystyle \int_1^y \dfrac{1}{\sqrt{x-1}}\dfrac{1}{x^{3/2}}\mbox{d}x & = & \displaystyle \int_0^{y-1} \dfrac{1}{\sqrt{\omega}}\dfrac{1}{(\omega+1)^{3/2}}\mbox{d}\omega + \int_0^\sqrt{y-1} \dfrac{1}{(\omega^2+1)^{3/2}}2\mbox{d}\omega \\ & = & \displaystyle \int_0^\sqrt{2(y-1)} \dfrac{1}{(1+t^2/2)^{3/2}}2\mbox{d}t\cdot \end{array} \]

Isso implica que \(\ell(x)\) corresponde à densidade de \((1+T^2/2)\) quando \(T\sim T_3\), ou seja \[ \ell(x) = \dfrac{\sqrt{2}\Gamma(3/2)/\sqrt{2}}{\sqrt{x-1}x^{3/2}} \pmb{1}_{(0,\infty)}(x)\cdot \] Pode-se verificar se esta é a constante de normalização correta executando integrate. A comparação da amostragem defensiva com o amostrador de importância original consiste, portanto, em adicionar uma pequena amostra de \(\ell\) à amostra original de \(g = f\): 

set.seed(28947)
sam1 = rt(.95*10^4,df=2)
sam2 = 1+.5*rt(.05*10^4,df=2)^2
sam = sample(c(sam1,sam2),.95*10^4)
h=function(x){(x>1)/sqrt(abs(x-1))}
weit = dt(sam,df=2)/(0.95*dt(sam,df=2)+
                       0.05*(sam>0)*dt(sqrt(2*abs(sam-1)),df=2)*sqrt(2)/sqrt(abs(sam-1)))
plot(cumsum(h(sam1))/1:length(sam1),type="l",ylim=c(0.2,0.6),ylab="",xlab="Iterations")
lines(cumsum(weit*h(sam))/(1:length(sam1)),col="blue")
grid()

Figura. 3.8. Convergência de dois estimadores da integral (3.9) do Exemplo 3.9 com base em uma amostra de \(T_2\) (linha escura) e uma versão defensiva (linha azul).

Observe que as simulações menores que 1 obtêm um peso igual a 1/0.95 na versão defensiva, pois \(\ell(x) = 0\) para \(x\leq 1\). Como no Exemplo 3.8 e na Figura 3.7, a amostra original pode apresentar saltos importantes na média acumulada que são avisos de problemas de variância infinita.

A solução de amostragem defensiva produz uma avaliação muito mais estável da integral. Em simulações alternativas, ambos os gráficos de convergência também podem ser bastante semelhante se nenhuma simulação for próxima o suficiente de 1 para induzir um grande valor de \(1/\sqrt{x−1}\). Na Figura 3.8, selecionamos uma ocorrência de discrepância entre ambas as amostras onde a amostragem defensiva traz um claro elemento de estabilização.


O exemplo acima ilustra claramente o impacto da amostragem defensiva quando o componente de cauda pesada da mistura está de alguma forma relacionado ao problema em questão. Escolhas genéricas de \(\ell\) geralmente levam a soluções menos eficientes, mesmo quando garantem uma variância finita para o estimador Monte Carlo.


Exemplo 3.10

Este exemplo considera um modelo probit de um ponto bayesiano de vista. Lembramos que o modelo probit é um caso particular de modelo linear generalizado onde os observáveis \(y\) são variáveis binárias, assumindo valores 0 e 1, e as covariáveis são vetores \(x\in\mathbb{R}^p\) tais que \[ P(y=1 \, | \, x) = 1-P(y=0 \, | \, x) = \Phi(x^\top \beta), \qquad \beta\in\mathbb{R}^p\cdot \]

Os dados para este modelo podem ser facilmente simulados, mas em vez disso usamos um conjunto de dados R chamado Pima.tr que também está disponível na biblioteca MASS. Este conjunto de dados pesquisa 200 mulheres indianas Pima em termos de presença ou ausência de diabetes Pima.tr$type, esta é a variável binária \(y\) para explicar, e várias covariáveis fisiológicas. Para fins de ilustração, consideramos apenas a covariável do índice de massa corporal Pima.tr$bmi, com um intercepto.

Uma estimativa GLM padrão do modelo é fornecida por:

Pima.tr=read.csv("http://leg.ufpr.br/~lucambio/SSim/Pima.csv", sep=";")
head(Pima.tr)
##   npreg glu bp skin  bmi   ped age type
## 1     5  86 68   28 30.2 0.364  24   No
## 2     7 195 70   33 25.1 0.163  55  Yes
## 3     5  77 82   41 35.8 0.156  35   No
## 4     0 165 76   43 47.9 0.259  26   No
## 5     0 107 60   25 26.4 0.133  23   No
## 6     5  97 76   27 35.6 0.378  52  Yes
Pima.tr$type = factor(Pima.tr$type)
ajuste=glm(type~bmi, data=Pima.tr, family=binomial(link="probit"))
summary(ajuste)
## 
## Call:
## glm(formula = type ~ bmi, family = binomial(link = "probit"), 
##     data = Pima.tr)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5788  -0.9242  -0.6442   1.2506   1.9661  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.54303    0.54211  -4.691 2.72e-06 ***
## bmi          0.06479    0.01615   4.011 6.05e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 256.41  on 199  degrees of freedom
## Residual deviance: 239.55  on 198  degrees of freedom
## AIC: 243.55
## 
## Number of Fisher Scoring iterations: 4
o que indica que a covariável do índice de massa corporal tem um impacto significativo na possível presença de diabetes.

De uma perspectiva Bayesiana, introduzimos um a priori vago em \(\beta=(\beta_1,\beta_2)\) que é uma distribuição \(N(0,100)\) normal. A distribuição a posteriori é então o produto deste essencialmente a priori pela verossimilhança, definida na função like abaixo.

Experimentar com a função de imagem e essa verossimilhança indica que a parte central da verossimilhança está localizada perto do estimador de verossimilhança máxima (MLE) com um intervalo de \(-0.6/-0.3\) para o intercepto \(\beta_1\) e um intervalo de \(0.04/0.09\) para \(\beta_2\). Usar uma proposta normal centrada no MLE com uma matriz de covariância diagonal correspondente à estimativa fornecida por glm é uma escolha natural para \(g\), embora isso não garanta uma variância finita para todos os propósitos.

Essa implementação mostra que os pesos de importância são bastante desiguais, se não degenerados, você pode verificar usando boxplot(weit), por exemplo. Uma representação de \(10^4\) pontos reamostrados com base nesses pesos na Figura 3.9 confirma esse padrão.

Para avaliar o baixo impacto de uma implementação de amostragem defensiva, também criamos uma amostra de importância que inclui simulações da a priori com probabilidade 0.05. A diferença na eficiência não é visível, no entanto. Ao usar o critério de tamanho amostral efetivo, definido na Seção 4.4, a diferença em \(10^3\) simulações é um tamanho amostral efetivo de 302 para a amostra normal versus um tamanho amostral efetivo de 283 para a defensiva.

As estimativas produzidas por ambos os métodos são \((-0.452,0.0653)\) e \((-0.452,0.0652)\), respectivamente. Observe a proximidade com o MLE se incorporarmos a média de Pima.tr$bmi. A razão para essa forte semelhança é que o termo adicional no denominador devido à inclusão da densidade a priori é praticamente zero. 

N = 10^3
like = function(beda) {
  mia = mean(Pima.tr$bm)
  prod(exp(beda[1] + (Pima.tr$bm[Pima.tr$t == "Yes"] - mia) * beda[2]))/
    prod(1 + exp(beda[1] - (Pima.tr$bm - mia) * beda[2]))
}
library(MASS)
bmi = Pima.tr$bm - mean(Pima.tr$bm)
glm(Pima.tr$t ~ bmi, family = binomial)
## 
## Call:  glm(formula = Pima.tr$t ~ bmi, family = binomial)
## 
## Coefficients:
## (Intercept)          bmi  
##     -0.7249       0.1048  
## 
## Degrees of Freedom: 199 Total (i.e. Null);  198 Residual
## Null Deviance:       256.4 
## Residual Deviance: 240   AIC: 244
be1 = seq(-1.2, -0.3, le = 100)
be2 = seq(0.02, 0.18, le = 130)
li = matrix(0, ncol = 130, nro = 100)
for (i in 1:100) for (j in 1:130) li[i,j] = like(c(be1[i], be2[j]))
image(be1, be2, li, xlab=expression(beta[1]),ylab=expression(beta[2]))
contour(be1, be2, li,add=TRUE)

sim = cbind(rnorm(N, m = -0.72, sd = 0.55), rnorm(N, m = 0.1, sd = 0.212))
weit = apply(sim, 1, like)/(dnorm(sim[, 1], m = -0.72, sd = 0.55) * 
                              dnorm(sim[, 2], m = 0.1, sd = 0.2121))
weit = weit/sum(weit)
vari1 = (1/(1 + exp(-sim[, 1] - sim[, 2] * bmi))) - sum((Pima.tr$t == "Yes"))/length(Pima.tr$bmi)
vari2 = (bmi/(1 + exp(-sim[, 1] - sim[, 2] * bmi))) - sum(bmi[Pima.tr$t == 
        "Yes"])/length(Pima.tr$bmi)
resim = sample(1:N, N, rep = TRUE, pro = weit)
reg = as.vector(lm(sim[resim, 1] ~ t(rbind(vari1[resim], vari2[resim])) - 1)$coef)
par(mfrow = c(2, 1), mar = c(4, 2, 2, 1))
plot(cumsum(weit * sim[, 1])/cumsum(weit), type = "l", ylim = c(-0.8, -0.6), 
     lwd = 2, xlab = "Iterations")
grid()
lines(cumsum(weit * (sim[, 1] - reg[1] * vari1 - reg[2] * 
        vari2))/cumsum(weit), col = "sienna", lwd = 2, lty = 2)
reg = as.vector(lm(sim[resim, 2] ~ t(rbind(vari1[resim], vari2[resim])) - 1)$coef)
plot(cumsum(weit * sim[, 2])/cumsum(weit), type = "l", ylim = c(0.08, 0.12), 
     lwd = 2, xlab = "Iterations")
grid()
lines(cumsum(weit * (sim[, 2] - reg[1] * vari1 - reg[2] * vari2))/cumsum(weit), 
      col = "sienna", lwd = 2, lty = 2)

Figura. 3.9. Distribuição a posteriori do parâmetro \((\beta_1,\beta_2)\) para a regressão do diabetes no índice de massa corporal no conjunto de dados Pima.tr com valores reamostrados de uma proposta normal sobreposta.



3.4 Exercícios