2 Muestreo por Importancia

Supongamos ahora que $f:\Rex^d \rightarrow \Rex$ y que se requiere evaluar la integral

$$I = \int_{\Theta} f(\bmath{\theta}) \, d \bmath{\theta}.$$

Claramente $I$ también puede escribirse como

$$I = \int_{\Theta} \left\{ \frac{f(\bmath{\theta})}{s(\bmath{\theta})} \right\} \, s(\bmath{\theta}) \, d \bmath{\theta},$$

donde $s(\bmath{\theta})$ es una función de densidad de probabilidad sobre $\Theta$. En otras palabras, $I = E_s(f(\bmath{\theta})/s(\bmath{\theta}))$. La distribución $s(\bmath{\theta})$ se conoce como la distribución de muestreo por importancia y generalmente se elige de manera que sea fácil de simular.

Si ahora generamos una muestra $\bmath{\theta}_1,\ldots,\bmath{\theta}_N$ de $s(\bmath{\theta})$ entonces podemos aproximar la integral $I$ a través del estimador insesgado

$$\hat{I}_{MI} = \frac{1}{N} \, \sum_{i=1}^N \frac{f(\bmath{\theta}_i)}{s(\bmath{\theta}_i)}.$$

La varianza de este estimador está dada por

$$\mbox{Var}(\hat{I}_{MI}) = \frac{1}{N} \, \mbox{Var}_s(f(\b... ...} \, \{ E_s(f(\bmath{\theta})^2/s(\bmath{\theta})^2) - I^2 \}.$$

La precisión de $\hat{I}_{MI}$ depende tanto del tamanõ de muestra, $N$, como de la distribución de muestreo por importancia, $s(\bmath{\theta})$. De hecho, si $s(\bmath{\theta})$ se elige proporcional a $f(\bmath{\theta})$ entonces Var $(\hat{I}_{MI}) = 0$ sin importar el tamaño de muestra. Por supuesto en la práctica dicha elección no es posible, pero esta idea sugiere que $s(\bmath{\theta})$ debe tener una forma similar a la de $f(\bmath{\theta})$, excepto tal vez en regiones donde los valores de $f(\bmath{\theta})$ sean despreciables.

Tenemos entonces que, para una integral $I$ dada, existe una infinidad de estimadores insesgados, en principio con precisiones distintas. Un aspecto importante del método de Monte Carlo se refiere al diseño de técnicas de reducción de varianza para dichos estimadores. Una de las técnicas más sencillas consiste precisamente en elegir una distribución de muestreo por importancia adecuada. Generalmente se requiere que $s(\bmath{\theta})$ satisfaga las siguientes condiciones:

(a) debe ser fácil de simular;

(b) debe tener una forma similar a la de $f(\bmath{\theta})$, la función que se desea integrar;

(c) debe tener las colas más pesadas que $f(\bmath{\theta})$, pues de otra forma la varianza de $\hat{I}_{MI}$ podría llegar a ser muy grande o incluso infinita.

En la práctica es común trabajar en términos de alguna reparametrización $\bmath{\varphi} = \bmath{\varphi}(\bmath{\theta})$, de manera que la integral esté definida sobre todo $\Rex^d$. Es este caso, el uso de la distribución $t$ de Student (con pocos grados de libertad) como distribución de muestreo por importancia es bastante frecuente.

Como una aplicación interesante de este método, consideremos el problema de calcular el valor esperado de una función $g(\bmath{\theta})$ respecto a la distribución final $p(\bmath{\theta} \vert \bmath{x})$, i.e.

$$E(g(\bmath{\theta}) \vert \bmath{x}) = \frac{ \int g(\bmath... ...th{\theta} } { \int p_x(\bmath{\theta}) \, d \bmath{\theta} }.$$ (10)

Sea $s(\bmath{\theta})$ una distribución de muestreo por importancia. Estimando por separado cada una de las integrales en (10) tenemos

$$\begin{eqnarray*} \int g(\bmath{\theta}) \, p_x(\bmath{\theta}) \, d \bmath{\the... ..._{i=1}^N \frac{ p_x(\bmath{\theta}_i) }{ s(\bmath{\theta}_i) }, \end{eqnarray*}$$

donde $\bmath{\theta}_1,\ldots,\bmath{\theta}_N$ es una muestra de $s(\bmath{\theta})$. El estimador resultante puede escribirse entonces como

$$\hat{E}_{MI}(g(\bmath{\theta}) \vert \bmath{x}) = \sum_{i=1}^N w(\bmath{\theta}_i) \, g(\bmath{\theta}_i),$$ (11)

donde los pesos $\{ w(\bmath{\theta}_i) \}$ están dados por

$$w(\bmath{\theta}_i) = \frac{ \left\{ p_x(\bmath{\theta}_i) ... ...a}_j) \right\} } \; \; \; \; \; \; \; \; \; \; (i=1,\ldots,N).$$

Este estimador es sólo asintóticamente insesgado. En problemas de este tipo es común elegir $s(\bmath{\theta})$ de acuerdo con la forma de $p_x(\bmath{\theta})$, es decir, de la distribución final, pero con las colas más pesadas. Si las colas de $s(\bmath{\theta})$ son más ligeras que las de $g(\bmath{\theta}) \, p_x(\bmath{\theta})$ entonces $\{ g(\bmath{\theta}) \, p_x(\bmath{\theta}) \} / s(\bmath{\theta})$ puede llegar a ser muy grande, lo cual haría que el estimador (11) fuera inestable. Por esta razón es conveniente que los pesos sean todos del mismo orden. De hecho, si $s(\bmath{\theta}) \propto p_x(\bmath{\theta})$ entonces $w(\bmath{\theta}_i) = 1/N$ para toda $i=1,\ldots,N$.

Lo anterior sugiere una manera de verificar la convergencia del estimador (11): comparar los pesos $\{ w(\bmath{\theta}_i) \}$ con una distribución uniforme discreta sobre $\{ 1,\ldots,N \}$, digamos a través de la divergencia logarítmica de Kullback-Leibler, i.e.

$$- \frac{1}{N} \sum_{i=1}^N \log( N \, w(\bmath{\theta}_i) ).$$

La divergencia logarítmica mide la discrepancia entre dos distribuciones. Es no negativa y se anula si y sólo si las distribuciones comparadas son iguales.