Quasi Monte Carlo or Gaussian Hermite quadrature for a statistical model with random effect

Question

In my likelihood function, I need to integrate a random effect out as follows

$$\int g(x,c)\exp(-c^2/2)dc .$$

Since the likelihood function is really complicated, I need the approximation to be fast and fairly accurate. I use Matlab to code it and have tried two approaches to approximate, but I can't figure out which one is the correct way.

Firstly, I have tried the quasi Monte Carlo integral. I generate $c$ as follows:

sim_max=50;
q=qrandstream('halton',1,'Skip',1e3,'Leap',1e2);
RandMat=qrand(q,sim_max);
c=norminv(RandMat,0,1);

and then I calculate the average.

Secondly, I am trying Matlab built-in function integral(fun,-10,10). Does this function implement Gaussian Hermite quadrature? Is the integration from -10 to 10 sufficient?

The problem is that the two methods give different results, and I don't know which one is more accurate.

Thank you for sharing you suggestions.

Answer 1

Your numerical results should converge

• as you use more and more sampling points
• as you sample a larger and larger fraction of the true integration range (presumably $\mathbb{R}$)

So you can plot the numerical result vs. number of sampling points used, and judge from this how many points you need. You can do this for different methods, and use that to judge which is the best for your problem.

You can also plot the numerical result vs. the width of your integration interval (20 in your case) and verify that your interval is wide enough for convergence to occur.