5
$\begingroup$

I am trying implement importance sampling of this integral $$ \mathfrak{I} = \int\limits_{-\infty }^{\infty }{\sqrt{\left| \frac{\theta }{1-\theta } \right|}}f(\theta )\text{d}\theta $$ where $f(\theta )\propto {{(1+{{\theta }^{2}}/5)}^{-3}}$ is a t-distribution with df=5.

I already sampled from the above distribution and was told to use samples from 1 and 2 below:

  1. Importance function equal to $$ 0.5\{{{g}_{1}}(\theta )+{{g}_{2}}(\theta )\}, $$ where $$ {{g}_{1}}(\theta )=\frac{1}{\pi }\frac{1}{1+{{\theta }^{2}}} $$ and
    $$ {{g}_{2}}(\theta )=\frac{1}{4\sqrt{\left|1-\theta \right| }} \quad\text{on}\quad\space [0,2] $$

  2. Importance function equal to $$ g'(\theta)\propto \frac{1}{\sqrt{\left| 1-\theta \right|}}\exp (-\left| 1-\theta \right|) $$

How should I do this? I have searched around, and vaguely understands concept. Could someone explain more in detail what I should do? Seems like I have all the tools I need.

$\endgroup$
  • 1
    $\begingroup$ The first suggestion is wrong: $g_2$ is not defined on the real line. See Example 3.13 of our book Monte Carlo Statistical Methods where this problem is fully processed! $\endgroup$ – Xi'an Feb 13 '12 at 21:04
  • $\begingroup$ @Xi'an: wow, you have keen eyes, I missed out typing $for \theta \in [0,2]$. Actually, I have a general question, in this formula $$h(x)=\int{h(x)\frac{f(x)}{g(x)}g(x)dx=E\left[ \frac{h(X)f(X)}{g(X)} \right]}$$ $$\simeq \frac{1}{n}\sum\limits_{j=1}^{n}{\frac{f\left( {{X}_{j}} \right)}{g\left( {{X}_{j}} \right)}h\left( {{X}_{j}} \right)}$$ Where do you find your $h\left( {{x}_{j}} \right)$ ? $\endgroup$ – user1061210 Feb 13 '12 at 22:39
  • 1
    $\begingroup$ 1. $\theta\in [0,2]$ still does not work in $\sqrt{1-\theta}$! $\endgroup$ – Xi'an Feb 14 '12 at 6:40
  • 1
    $\begingroup$ I see, I didn't notice it is a folded gamma. I did this $$F(x)=\left\{ \begin{array}{*{35}{l}} \int\limits_{0}^{1}{\frac{1}{\sqrt{1-x}}\exp \left( -(1-x) \right)dx} & \text{if x}<1 \\ \int\limits_{1}^{\infty }{\frac{1}{\sqrt{x-1}}\exp \left( -(x-1) \right)}dx & \text{if x}>1 \\ \end{array} \right.$$ $$\int\limits_{0}^{1}{\frac{1}{\sqrt{1-x}}{{e}^{-(1-x)}}}dx=-\int\limits_{1}^{0}{2{{e}^{-{{u}^{2}}}}}=\sqrt{\pi }erf(1)=1.49$$ $$\int\limits_{1}^{\infty }{\frac{1}{\sqrt{x-1}}\exp \left( -(x-1) \right)}dx=\int\limits_{0}^{\infty }{2{{e}^{-{{u}^{2}}}}}du=\sqrt{\pi }=1.77$$ 1.49+1.77=3.27 $\endgroup$ – user1061210 Feb 19 '12 at 22:52
  • 1
    $\begingroup$ I've also been reading your other book"Introducing Monte Carlo Methods with R", pretty good examples, I find it very practical for learning. $\endgroup$ – user1061210 Feb 20 '12 at 17:25
5
$\begingroup$

In case your difficulty is with the simulation per se, here is my R code to compare simulations from $f$ (plain), $g$ equal to $$ \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|}} $$ (mixture of Cauchy and power distributions) and $m$ equal to $$ \frac{1}{2}\frac{1}{\Gamma(1/2)}\frac{1}{\sqrt{|1-x|}}\exp\{-|1-x|\} $$ (folded Gamma).

Simulating from $f$ is straightforward

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

where

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

Simulating from $g$ requires simulating from the square-root part. If you integrate out $1/4{\sqrt{|1-x|}}$ over $[0,2]$, you get either $1-\sqrt{1-x}$ over $[0,1]$ or $\sqrt{x-1}$ over $[1,2]$, which means that this distribution can be represented as $$ 1\pm \mathcal{U}(0,1)^2. $$ (In the following, I force both subsamples to have the same size $5\cdot10^5$, which is a Rao-Blackwellisation trick to reduce the variance with no impact on the expectation.)

> 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)

where

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

Simulating from $m$ follows from the folded representation:

> 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))

The comparison of the three simulation methods is illustrated in the following boxplot (that we should use in the next edition of Monte Carlo Statistical Methods!)

Boxplot of the variations of three estimators of \mathfrak{I}

$\endgroup$
  • $\begingroup$ I will run this in R, if I got questions, I'll leave a msg, thanks again! $\endgroup$ – user1061210 Feb 20 '12 at 17:26
  • 1
    $\begingroup$ Here the definition of sam3 takes 5 times as many samples as the other tests. This makes the folded Gamma look better than the Cauchy-power mixture, while the reverse is true. The folded Gamma estimator has infinite variance, while the Cauchy-power mixture does not. $\endgroup$ – deinst Feb 25 '12 at 16:46
  • $\begingroup$ Indeed, the mixture does better now. About the infinite variance, we noticed this (major) drawback in the book. However, it does not seem to impact the practical implementation (in that we do not see strange realisations in the range of the estimates)... $\endgroup$ – Xi'an Feb 25 '12 at 18:11
  • 1
    $\begingroup$ The problems with the variance are akin to the game where you win $2^n/n$ dollars with probability $2^{-n}$. The expected payout may be infinite, but I am not about to sell my house to play. With the folded gamma there will be large excursions, but they will be extremely rare. In practice using $\frac{\sum hf/g}{\sum f/g}$ with the folded gamma performs worse, although in theory it has finite variance. $\endgroup$ – deinst Feb 25 '12 at 19:23
  • 1
    $\begingroup$ @Xi'an Sorry to bother you again, but your definition of h contain the t distribution, when it should just have the $\sqrt{|x/(1-x)|}$. This causes all the estimates to be off because you are multiplying by f^2 instead of f. The approximations should be around 1.015, while they are around 0.25. $\endgroup$ – deinst Feb 27 '12 at 3:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.