4. Controlando e acelerando a convergência

No Capítulo 3, o método Monte Carlo foi apresentado e discutido como uma abordagem baseada em simulação para a aproximação de integrais complexas. Embora os princípios já devam ser bem compreendidos, há mais a ser dito sobre a avaliação da convergência; ou seja, quando e por que parar de executar simulações.

Apresentamos neste capítulo as especificidades da estimativa de variância e controle para métodos Monte Carlo, bem como dispositivos de aceleração. Nas Seções 4.2 e 4.5, focamos particularmente na construção de bandas de confiança, enfatizando as limitações das avaliações de base normal na Seção 4.2 e desenvolvendo estimativas de variância para amostradores de importância na Seção 4.3 e ferramentas de avaliação de convergência na Seção 4.4.

Esses são conceitos fundamentais e veremos conexões com desenvolvimentos semelhantes no domínio dos algoritmos MCMC, discutidos nos Capítulos 6–8. A segunda parte do capítulo cobre vários dispositivos de aceleração, como Rao-Blackwellization na Seção 4.6 e correlação negativa na Seção 4.7.



4.1 Introdução

No Capítulo 3 mencionamos que o Teorema do Limite Central se aplica às estimativas Monte Carlo da forma \[ \overline{h}_n = \dfrac{1}{n}\sum_{i=1}^n h(X_i), \qquad X_i \sim f, \] sob condições de integrabilidade e, portanto, que pode ser explorado para avaliar a convergência para a integral de interesse, \[\begin{equation} \tag{4.1} \mathcal{I} = \int h(x) f(x) \mbox{d} x, \end{equation}\] no sentido de que a variável aleatória \(\sqrt{n}(\overline{h}_n - \mathcal{I})\) é assintoticamente normal.

O painel inferior da Figura 3.3 associado ao Exemplo 3.3 fornece uma ilustração direta do uso de um intervalo de confiança normal para essa avaliação. Ele representa, para cada número fixo de iterações, um intervalo de confiança assintoticamente válido no valor de (4.1).

No entanto, essa abordagem tem limitações porque o envelope construído sobre as iterações e representado na Figura 3.3 não tem validade geral como uma faixa de confiança. De fato, se você repetir o experimento Monte Carlo mais uma vez, a sequência \(\{\overline{h}_n\}\) produzida na segunda execução provavelmente não ficará dentro desse envelope e, se você repeti-la várias vezes, a frequência com que ela permanecerá dentro da banda não atendem à probabilidade nominal de 0.95.

A explicação para esta aparente discrepância é que o método de monitoramento ilustrado pelo painel inferior da Figura 3.3 é essencialmente de natureza univariada. Ou seja, os limites de confiança colocados na estimativa \(\overline{h}_k\) na iteração \(k\) dependem apenas dos valores de \(\overline{h}_k\) e da estimativa de variância no tempo \(k\) e, portanto, ignoram qualquer estrutura de correlação nas iterações. Uma banda de confiança válida precisa dar conta da distribuição de toda a sequência, em uma perspectiva multivariada ou mesmo funcional, conforme discutido na Seção 4.5, embora a construção teórica de tal banda global esteja de alguma forma fora de nosso alcance.

Ressaltamos que a faixa de confiança fornecida pelo Teorema do Limite Central é assintoticamente válida. O que estamos buscando é uma avaliação mais global da convergência para a sequência de estimadores conforme varia o número de simulações.

Isso pode ser considerado como uma avaliação de convergência de segunda ordem, necessariamente mais conservador, ou seja, exigindo mais simulações do que o original. Embora não explicitemos a ligação aqui, a largura do lote fixa média e os métodos bootstrap de Jones et al. (2006) descritos na Seção 8.4.4 também se aplicam às configurações iid do capítulo atual.


4.2 Variação de monitoramento

Primeiro, você deve perceber que existe uma solução direta e simples para avaliar a variabilidade de uma sequência de estimativas de Monte Carlo, que é executar várias sequências independentes em paralelo. Isso é mais fácil de derivar e mais amplamente aplicável do que as técnicas baseadas em aproximação assintótica, embora muito mais ávidas em tempo de computação.

Infelizmente, esta última característica é uma característica que encontraremos repetidamente no livro, ou seja, que a validação da avaliação da variação é de ordem superior à convergência do próprio estimador. Ou seja, a avaliação do erro requer muito mais tempo de computação do que a validação da convergência pontual (exceto em casos muito especiais, como a regeneração, abordada em Robert e Casella, 2004). Uma versão aproximada, porém mais barata, dessa estimativa básica de Monte Carlo da variabilidade é inicializar (consulte a Seção 1.5) a amostra atual, conforme já utilizado no Exemplo 3.7 (consulte o painel inferior esquerdo da Figura 3.6).


Exemplo 4.1

Se repetirmos as simulações do Exemplo 3.3, podemos produzir uma matriz de estimadores convergentes como em: 

#x = matrix(h(runif(200*10^4)),ncol=200)
#estint = apply(x,2,cumsum)/(1:10^4)


No Exemplo 4.1, o apelo de usar uma faixa de confiança bootstrap é de alguma forma limitado porque o custo computacional de produzir uma sequência bootstrap é aproximadamente o mesmo que o custo computacional de produzir uma nova sequência. Em configurações mais complexas, no entanto, produzir uma nova sequência pode ser muito mais caro do que reamostrar a partir da sequência original.

Este exemplo simples nos adverte contra o uso cego de uma aproximação normal quando invocada repetidamente em iterações com estimadores dependentes simplesmente porque a aproximação de confiança normal tem apenas uma validação pontual. Usar uma banda de estimadores em paralelo é obviamente mais custoso, mas fornece a avaliação correta sobre a variação desses estimadores.


4.3 Variância assintótica dos estimadores da amostragem de importância

O exemplo a seguir ilustra uma dificuldade básica ao avaliar a convergência para amostragem por importância.


Exemplo 4.2

Para uma observação normal igual a \(x = 2.5\) com uma distribuição a priori de Cauchy em sua média, \[ X\sim N(\theta,1), \qquad \theta \sim C(0,1) \] a esperança a posteriori da média \(\theta\) é dada por \[ \delta^\pi(x) = \int_{-\infty}^\infty \dfrac{\theta}{1+\theta^2}e^{-\frac{1}{2}(x-\theta)^2} \mbox{d}\theta \Big/ \int_{-\infty}^\infty \dfrac{1}{1+\theta^2}e^{-\frac{1}{2}(x-\theta)^2} \mbox{d}\theta\cdot \]

Portanto, uma aproximação de \(\delta^\pi(x)\) pode ser baseada na simulação de variáveis independentes identicamente distribuídas \(\theta_1,\cdots,\theta_n \sim N(x,1)\) como \[ \widehat{\delta}^\pi_n(x)=\sum_{i=1}^n \dfrac{\theta_i}{1+\theta_i^2} \Big/ \sum_{i=1}^n \dfrac{1}{1+\theta_i^2}, \] já que tanto o numerador quanto o denominador são convergentes em \(n\).

Observe que o estimador também pode ser interpretado como uma aproximação da amostragem de importância auto-normalizada com a função de importância sendo a densidade normal e a razão de importância sendo igual a \(1/(1+\theta^2)\).


Uma dificuldade associada a este exemplo e a qualquer outro amostrador de importância auto-normalizado é que o estimador \(\widehat{\delta}^\pi_n(x)\) é de fato uma razão de estimadores e que a variância de uma razão não é a razão das variâncias. Esta é uma ocorrência comum em cálculos Bayesianos, veja os Exercícios 4.1 e 4.2.

Como mencionado anteriormente, as constantes de normalização são supérfluas na inferência bayesiana, exceto no caso em que vários modelos são considerados simultaneamente para serem comparados. Nessas ocorrências, a grandeza fundamental é o fator de Bayes \[ \rho = \dfrac{m_1(x)}{m_2(x)} = \dfrac{\displaystyle \int_{\Theta_1} \pi_1(\theta_1) f_1(x|\theta_1)\mbox{d}\theta_1}{\displaystyle \int_{\Theta_2} \pi_2(\theta_2) f_2(x|\theta_2)\mbox{d}\theta_2}, \] cuja posição contra 1 conduz a comparação ver, por exemplo, Robert (2001).

Neste caso, as distribuições a posteriori para ambos os modelos em comparação estão normalmente disponíveis até uma constante de normalização, \[ \pi_1(\theta|x) = \widetilde{\pi}_1(\theta)/c_1 \qquad \mbox{e} \qquad \pi_2(\theta|x) = \widetilde{\pi}_2(\theta)/c_2, \] onde apenas \(\widetilde{\pi}_1\) e \(\widetilde{\pi}_2\) são conhecidos e onde \(c_1\) e \(c_2\) correspondem às verossimilhanças marginais, \(m_1(x)\) e \(m_2(x)\), a dependência de \(x\) é removida para fins de simplificação. O fator de Bayes é, portanto, idêntico à razão dessas constantes ausentes, \(\varrho = c_1 /c_2\). Técnicas computacionais especiais foram desenvolvidas para a aproximação dos fatores de Bayes, como em Chen et al. (2000), ver também Exercício 4.2.

A transformada \(\psi = \log(\varrho)\) costuma ser vista como mais relevante para a comparação do modelo, devido à escala \(\chi^2\) da razão log-verossimilhança e essa razão log-odds pode ser aproximada por si só.

Considere assim um estimador de razão geral \[ \delta_h^n = \sum_{i=1}^n \omega_i h(x_i) \Big/ \sum_{i=1}^n \omega_i, \] onde os \(x_i\)’s são realizações de variáveis aleatórias \(X_i\sim g(y)\) com \(g\) como uma distribuição candidata para o alvo \(f\). Além disso, os \(\omega_i\)’s são realizações de variáveis aleatórias \(W_i\) tais que \(\mbox{E}\big(W_i |X_i = x\big) = \kappa f(x)/g(x)\), sendo \(\kappa\) uma constante arbitrária que corresponde à falta de constantes de normalização em \(f\) e \(g\).


Exemplo 4.3 (Continuação do Exemplo 4.2)

Se gerarmos uma amostra normal \(N(x, 1)\) para a aproximação da amostragem de importância (4.2), a aproximação de variância acima pode ser usada para avaliar a variabilidade dessas estimativas mas, novamente, a natureza assintótica dessa aproximação deve ser levada em consideração.

Se tomarmos como referência o intervalo de 500 sequências paralelas de estimadores de \(\delta^\pi(x)\), a justaposição da banda produzida em uma única sequência mostra uma subestimação da variância se a estimativa da variância usual for usada, já que os respectivos intervalos são 0.0559 para o primeiro (Monte Carlo) e 0.0539 para o segundo (variância corrigida). A Figura 4.2 compara as bandas de variância aproximadas para \(x = 2.5\) com a variância real das estimativas, avaliadas nas 500 sequências paralelas. 

set.seed(362019)
norma = matrix(rnorm(500*10^4),ncol=500)+2.5
weit = 1/(1+norma^2)
esti = apply(norma*weit,2,cumsum)/apply(weit,2,cumsum)
plot(esti[,1],type="l",col="red",ylim=c(1.5,2.0),xlab = "Iterations",ylab = "Estimators posterior mean")
band = apply(esti,1,quantile,c(.025,.975))
polygon(c(1:10^4,10^4:1),c(band[1,],rev(band[2,])),col = "cyan")
vare = cumsum(weit[,1]*norma[,1]^2)/cumsum(weit[,1])-esti[,1]^2
lines(esti[,1]+2*sqrt(vare/(1:10^4)))
lines(esti[,1]-2*sqrt(vare/(1:10^4)))
varw = cumsum(weit[,1]^2)*(1:10^4)/cumsum(weit[,1])^2
lines(esti[,1]+2*sqrt(varw*vare/(1:10^4)),col="sienna")
lines(esti[,1]-2*sqrt(varw*vare/(1:10^4)),col="sienna")
grid()

Figura. 4.2. Convergência de uma sequência de estimadores Monte Carlo para a média a posteriori no problema Cauchy-Normal quando \(x = 2.5\) e as simulações são de uma função de importância \(N (x, 1)\). A zona sombreada representa a faixa de confiança de 95% em todo o conjunto de 500 sequências paralelas de estimadores Monte Carlo em cada iteração, o limite interno corresponde à banda normal para a estimativa da variância padrão e o limite externo (mais claro) corresponde à estimativa de variância corrigida para a Normal.

A Figura 4.3 reproduz esta avaliação para \(\theta_i\)’s simulados da distribuição \(C(0, 1)\) a priori e para a estimativa de amostragem de importância \[ \widetilde{\delta}_h^n = \sum_{i=1}^n \theta_i e^{-\frac{(x-\theta_i)^2}{2}} \Big/ \sum_{i=1}^n e^{-\frac{(x-\theta_i)^2}{2}}\cdot \] A inversão dos papéis das distribuições \(N (x, 1)\) e \(C(0, 1)\) ilustra mais uma vez a ambigüidade da representação integral e as oportunidades oferecidas pela amostragem por importância.

Neste caso, como as correspondentes funções \(h\) são limitadas para ambas as escolhas, as variabilidades das estimativas são bastante semelhantes, com uma pequena vantagem para a amostragem Normal. O intervalo da avaliação Monte Carlo é de fato 0.0594, enquanto o intervalo da banda assintótica é 0.0800. Deve-se levar em consideração o fato de que este último é baseado em uma única sequência. ,/p>

Outra sequência produziria um intervalo diferente, como pode-se verificar facilmente. Da mesma forma, as curvas de contorno que aparecem nas Figuras 4.2 e 4.3 são produzidas por uma única sequência e, portanto, não devem ser superinterpretadas.



4.4 Tamanho efetivo da amostra e perplexidade

A representação (4.3) proposta para a variância aproximada do estimador de amostragem de importância auto-normalizada conduz a uma ferramenta essencial para avaliar o desempenho de amostradores de importância. Ao expandir a variância empírica, \[ \widehat{\mbox{Var}}(W)=\dfrac{1}{n}\sum_{i=1}^n \omega_i^2 - \dfrac{1}{n^2}\left( \sum_{i=1}^n \omega_i\right)^2, \] o coeficiente \(1+\widehat{\mbox{Var}}_g(W)\) é igual a \[ n^2 \sum_{i=1}^n \omega_i^2 \Big/ \left( \sum_{i=1}^n \omega_i\right)^2\cdot \]

Se agora denotarmos os pesos normalizados por \[ \underline{\omega}_i = \omega_i \Big/ \sum_{j=1}^n \omega_j, \] então definimos o tamanho efetivo da amostra (ESS) por \[ \mbox{ESS}_n = 1\Big/ \sum_{i=1}^n \underline{\omega}_i^2\cdot \]

Além de ser útil para avaliar a perda de variância devido aos pesos de importância, esse fator fornece uma avaliação direta do valor do amostrador de importância, pois é equivalente a um tamanho de amostra. Para uma amostra uniformemente ponderada, \(\mbox{ESS}_n\) é igual a \(n\), enquanto, para uma amostra completamente degenerada onde todos os pesos de importância, exceto um, são zero, \(\mbox{ESS}_n\) é igual a \(1\). O tamanho efetivo da amostra, portanto, avalia o tamanho da amostra independete idênticamente distribuída equivalente à amostra ponderada e permite uma comparação direta de amostradores.

Uma segunda avaliação e quase equivalente é fornecida pela chamada perplexidade (Cappé et al., 2008), \(\exp(\hbar_n)/n\), onde \[ \hbar_n = -\sum_{i=1}^n \underline{\omega}_i \log(\underline{\omega}_i) , \] é a entropia de Shannon dos pesos de amostragem de importância normalizada. Essa ferramenta é usada principalmente na teoria da informação e no reconhecimento de fala, ver Jelinek (1999), mas ainda assim traz um recurso adicional para o tamanho efetivo da amostra que pode ser explorado em todas as configurações.

A perplexidade fornece uma estimativa de \(\exp\big(\Xi(f,g)\big)\), onde \[ \Xi(f,g) = \int \log\left( \dfrac{f(x)}{g(x)}\right) f(x) \mbox{d}x, \] é a divergência de Kullback–Leibler entre o alvo \(f\) e a função de importância \(g\). Portanto, quanto mais próxima a perplexidade estiver de 1, mais apropriada será a função de importância.


Exemplo 4.4. (Continuação do Exemplo 4.3)

div> Se compararmos os tamanhos efetivos de amostra para as simulações normal e de Cauchy, a diferença de eficiência (e perplexidade) é muito mais clara do que na comparação entre as Figuras 4.2 e 4.3. O cálculo dessas quantidades é simples.

Para os tamanhos de amostra efetivos,

#ess = apply(weit,2,cumsum)^2/apply(weit^2,2,cumsum)
#essbo = apply(ess,1,quantile,c(.025,.975))
#ech = apply(wachd,2,cumsum)^2/apply(wachd^2,2,cumsum)
#echbo = apply(ech,1,quantile,c(.025,.975))
#sumweit = apply(weit,2,cumsum)
#plex = (apply(weit*log(weit),2,cumsum)/sumweit)-log(sumweit)
#chumweit = apply(wachd,2,cumsum)
#plech = (apply(wachd*log(wachd),2,cumsum)/chumweit)-log(chumweit)
#plob = apply(exp(plex),1,quantile,c(.025,.975))
#ploch = apply(exp(plech),1,quantile,c(.025,.975))



4.5 Monitoramento simultâneo

Conforme mencionado na introdução deste capítulo, um método válido para anexar variâncias a um gráfico de média e, portanto, para construir um Teorema do Limite Central válido, é derivar os limites usando uma abordagem multivariada. Uma primeira construção, baseada em uma aproximação normal, é encontrada em Robert and Casella (2004).


4.6 Rao–Blackwellização e descondicionamento

Um famoso teorema de estatística matemática, o Teorema de Rao-Blackwell, afirma que a substituição de um estimador por uma expectativa condicional melhora sua variância, \[ \mbox{Var}\Big( \mbox{E}\big( \delta(X) | Y\big)\Big) \leq \mbox{Var}\big( \delta(X)\big) \] se \(Y\) é uma estatística suficiente (Lehmann and Casella, 1998). Esse teorema também tem influência na metodologia computacional, pois fornece uma abordagem genérica para reduzir a variância de um estimador Monte Carlo, que é usar condicionamento.

Essa técnica às vezes é chamada de Rao-Blackwellization (Gelfand and Smith, 1990, Liu et al., 1994, Casella and Robert, 1996), embora o condicionamento nem sempre seja em termos de estatísticas suficientes em configurações Monte Carlo. Ele basicamente afirma que o uso de esperanças condicionais, que podem ser calculadas, expressões Monte Carlo traz uma melhoria na variabilidade desses estimadores Monte Carlo, enquanto não perturba sua imparcialidade inerente.

Rao–Blackwellization significa aproveitar o fato de que, se \(\delta(X)\) é um estimador de \(\Im = \mbox{E}_f\big(h(X)\big)\) e se \(X\) pode ser simulado a partir da distribuição conjunta \(f^*(x, y)\) satisfazendo \[ \int f^*(x, y) \mbox{d}y = f(x), \] então o novo estimador \(\delta^*(Y) =\mbox{E}_f\big(\delta(X)|Y\big)\) domina \(\delta\) em termos de variância, enquanto o viés é o mesmo. Obviamente, este resultado só é útil em configurações onde \(\delta^*(Y)\) pode ser calculado explicitamente.


Exemplo 4.6.

div> Considere calcular a esperança de \(h(x) = \exp(−x^2)\) quando \(X\sim T(\nu,\mu,\sigma^2)\). A distribuição \(t\)-Student padrão foi construída por William Gosset (ver, por exemplo, Stigler, 1986) como a distribuição da razão de uma variável Normal por uma variável normalizada \(\chi\), não \(\chi^2\)!; isto é, se \(X\sim T (\nu,\mu,\sigma^2)\), então \[ X = \mu+ \sigma \dfrac{\epsilon}{\sqrt{\xi/\nu}}, \] com \(\epsilon \sim N(0,1)\) e \(\xi \sim \chi^2_\nu\). 

#y = sqrt(rchisq(Nsim,df=nu)/nu)
#x = rnorm(Nsim,mu,sigma/y)
#d1 = cumsum(exp(-x^2))/(1:Nsim)
#d2 = cumsum(exp(-mu^2/(1+2*(sigma/y)^2))/sqrt(1+2*(sigma/y)^2))/(1:Nsim)


Infelizmente, esse princípio de condicionamento parece ter aplicabilidade limitada, pois envolve um tipo particular de simulação (variáveis conjuntas) e também requer funções integrandos \(h\) suficientemente regulares para que as esperanças condicionais sejam explícitas. Existem, no entanto, situações específicas em que Rao-Blackwellization é sempre possível. Uma ilustração é fornecida em Robert and Casella (2004) na configuração dos métodos Aceitar-Rejeitar, onde as simulações rejeitadas podem ser recicladas integrando sobre as variáveis uniformes que são usadas no estágio de seleção.


4.7 Métodos de aceleração

Semelhante ao dispositivo de descondicionamento acima, existem estratégias de aceleração global que são mais ou menos independentes da configuração da simulação, mas tentam explorar a saída da simulação de maneiras mais eficientes e, portanto, devem ser implementadas sempre que possível. Quando vistos como dispositivos de pós-processamento, esses métodos também podem ser usados para avaliar a convergência, fornecendo soluções alternativas de estimadores da mesma quantidade.


4.7.1 Simulações correlacionadas

Embora os métodos usuais de simulação levem a amostras independentes idênticamente distribuídas, na verdade pode ser preferível gerar amostras de variáveis negativamente ou positivamente correlacionadas ao estimar uma integral \(\Im\), pois elas podem reduzir a variância do estimador correspondente.


4.8 Exercícios

4.1 Assumimos aqui que ambas as a posteriori são absolutamente contínuas um em relação ao outro, ou seja, que os parâmetros para ambos os modelos pertencem ao mesmo espaço paramátrico.

  1. Mostre que o fator de Bayes \(\varphi\) pode ser aproximado por \[ \widehat{\varrho}=\dfrac{1}{n}\sum_{i=1}^n \dfrac{\widetilde{\pi}_1(\theta_i)}{\widetilde{\pi}_2(\theta_i)}, \qquad \theta_1,\cdots,\theta_n \sim \pi_2\cdot \]
  2. Mostre que a identidade \[ \dfrac{\displaystyle \int \widetilde{\pi}_1(\theta)\alpha(\theta) \pi_2(\theta|x)\mbox{d}\theta}{\displaystyle \int \widetilde{\pi}_2(\theta)\alpha(\theta) \pi_1(\theta|x)\mbox{d}\theta} = \dfrac{c_1}{c_2} = \varrho, \] vale para toda função \(\alpha\) tal que ambas as integrais são finitas.

4.2 Na mesma hipótese do Exercício 4.1, quando as prioris π1 e π2 pertencem à mesma família parametrizada (ou seja, quando πi (θ) = π(θ|λi )), as constantes de normalização correspondentes são denotadas por c( λi). O parâmetro λ é um número real.

  1. Quando π(λ) é uma distribuição arbitrária em λ com suporte (λ1 , λ2 ), verifique a identidade

4.3 Mostre que \(\mbox{ESS}_n\) sempre assume valores entre 1 e \(n\).

4.4 Mostre que \(\exp(\hbar_n)\) sempre assume valores entre 1 e \(n\).