7
$\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.

Improve this question
$\endgroup$
11
  • 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, 2012 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, 2012 at 22:39
  • 1
    $\begingroup$ 1. $\theta\in [0,2]$ still does not work in $\sqrt{1-\theta}$! $\endgroup$
    – Xi'an
    Feb 14, 2012 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, 2012 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, 2012 at 17:25

1 Answer 1

Reset to default
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}

Improve this answer
$\endgroup$
5
  • $\begingroup$ I will run this in R, if I got questions, I'll leave a msg, thanks again! $\endgroup$
    – user1061210
    Feb 20, 2012 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, 2012 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, 2012 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, 2012 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, 2012 at 3:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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