¿Cómo calcular el muestreo de importancia?

Question

Estoy tratando de implementar el muestreo de importancia de esta integral $$\mathfrak{I} = \int\limits_{-\infty }^{\infty }{\sqrt{\left| \frac{\theta }{1-\theta } \right|}}f(\theta )\text{d}\theta$$ donde $f(\theta )\propto {{(1+{{\theta }^{2}}/5)}^{-3}}$ es una distribución t con df=5.

Ya tomé una muestra de la distribución anterior y me dijeron que utilizara las muestras de la 1 y la 2 siguientes:

1. Función de importancia igual a $$0.5\{{{g}_{1}}(\theta )+{{g}_{2}}(\theta )\},$$ donde $${{g}_{1}}(\theta )=\frac{1}{\pi }\frac{1}{1+{{\theta }^{2}}}$$ y
   $${{g}_{2}}(\theta )=\frac{1}{4\sqrt{\left|1-\theta \right| }} \quad\text{on}\quad\space [0,2]$$

2. Función de importancia igual a $$g'(\theta)\propto \frac{1}{\sqrt{\left| 1-\theta \right|}}\exp (-\left| 1-\theta \right|)$$

¿Cómo debo hacerlo? He buscado por ahí, y vagamente entiende el concepto. ¿Podría alguien explicar con más detalle lo que debo hacer? Parece que tengo todas las herramientas que necesito.

Answer 1

En caso de que su dificultad sea con la simulación por sí mismo Aquí está mi código R para comparar las simulaciones de $f$ (llano), $g$ igual a $$\frac{1}{2} \frac{1}{\pi} \frac{1}{1+x^2} + \frac{1}{2}\frac{1}{4}\frac{\mathbb{I}_{[0,2]}(x)}{\sqrt{|1-x|}}$$ (mezcla de distribuciones de Cauchy y de potencia) y $m$ igual a $$\frac{1}{2}\frac{1}{\Gamma(1/2)}\frac{1}{\sqrt{|1-x|}}\exp\{-|1-x|\}$$ (Gamma doblada).

Simulación de $f$ es sencillo

> sam1=matrix(rt(10^6,df=5),ncol=100)
> fam1=h(sam1)

donde

> h
function(x){ 
sqrt(abs(x/{1-x}))}

Simulación de $g$ requiere la simulación de la parte de la raíz cuadrada. Si se integra $1/4{\sqrt{|1-x|}}$ en $[0,2]$ , se obtiene $1-\sqrt{1-x}$ en $[0,1]$ o $\sqrt{x-1}$ en $[1,2]$ lo que significa que esta distribución puede representarse como $$1\pm \mathcal{U}(0,1)^2.$$ (A continuación, fuerzo que ambas submuestras tengan el mismo tamaño $5\cdot10^5$ , que es un truco de Rao-Blackwellización para reducir la varianza sin impacto en la expectativa).

> sam22=1+sample(c(-1,1),5*10^5,rep=TRUE)*runif(5*10^5)^2
> sam21=rcauchy(5*10^5)
> sam2=matrix(sample(c(sam21,sam22)),ncol=100)
> fam2=h(sam2)*dt(sam2,df=5)/g(sam2)

donde

g=función(x){.5*dcauchy(x)+.125*((x>0)*(x<2))/sqrt(abs(1-x))}

Simulación de $m$ se desprende de la representación plegada:

> sam3=matrix(1+sample(c(-1,1),10^6,rep=TRUE)*rgamma(10^6,.5),ncol=100)
> fam3=h(sam3)*dt(sam3,df=5)/(.5*dgamma(abs(1-sam3),.5))

La comparación de los tres métodos de simulación se ilustra en el siguiente boxplot (que deberíamos utilizar en la próxima edición de Métodos estadísticos de Monte Carlo !)