1
votes

I am trying to implement Monte Carlo integration with importance sampling. I've created a trivial example - I wish to integrate h(x), which is has a student t-distribution (mu =1, sigma =1, df=100), but us scaled-up 4 times - I want to integrate over the interval [-2,2] - f(x), the pdf of my h(x), is then t(1,1,100) - my proposal distribution is g(x), and normal (0,1)

I can't get this work... I am confused with how to implement Importance Sampling and using a probability-density-function for h(x) and a proposal distribution g(x). I'm sure my implementation is wrong. I was hoping someone could please help me?

xtemp<-rnorm(100000) 
x<-xtemp[which(xtemp>=-2 & xtemp<=2)] 
hx<-dt(x,100)*4
fx<-dt(x,100) 
gx<-dnorm(x) 

IntMC<-sum(hx*fx/gx)/length(hx) 
IntAn <-(pt(2,100)-pt(-2,100))*4
1
What is IS? pdf? Get rid of acronims, please. - Stepan Yakovenko
I'm going to go out on a limb and guess that IS means Importance Sampling, and pdf is Probability Density Function. - Matthew Lundberg
edits done, pdf = probability density function, IS = importance sampling - user1375871
I'd recommend taking a little time to read Numerical Recipes in C. There are pdf versions available for download. You don't have to use the c-code, as the algorithms themselves are presented in excellent detail. - Carl Witthoft

1 Answers

2
votes

Since you are sampling from a truncated normal distribution, you should not use the probability density function of the normal distribution (dnorm in your example) but of the truncated normal distribution (e.g., dtnorm from the package msm) for the computation of the weights.

Try the following, it gives you the expected result:

library(msm)  # provides the pdf of the truncated normal distribution

xtemp <- rnorm(100000) 
x <- xtemp[which(xtemp>=-2 & xtemp<=2)] 

hx <- dt(x,100)*4
gx <- dtnorm(x, lower=-2, upper=2) 

IntMC <- mean(hx/gx) 
IntAn <- (pt(2,100)-pt(-2,100))*4