Computing an integral with Monte Carlo method

Question

I have to write a MATLAB code to approximate the following integral using Monte Carlo integration: $$ \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}(ax^3+bx^2+cx+d)e^{-\frac{1}{2}x^2}dx $$ with $a=2,\ b=0,\ c=1,\ d=10$ and using $n=1000$ realizations. My attempt is the following

ApproxInt = mean(a * rand(1000, 1).^3+b*rand(1000, 1).^2+c * rand(1000,1)+d)

and it gives perhaps $11$ as approximated result of the integral. Is this correct?

Answer 1

Presumably the integral is

$$ \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}(ax^3+bx^2+cx+d)e^{-\frac{1}{2}x^2}dx=b+d $$

Note that $\frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty x^3 e^{-\frac{1}{2}x^2}dx =0$, $\frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty x^2 e^{-\frac{1}{2}x^2}dx =1$, $\frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty x e^{-\frac{1}{2}x^2}dx =0$, $\frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty e^{-\frac{1}{2}x^2}dx =1$ by the moments of the gaussian distribution.

The correct way to calculate this integral using $n$ points in MATLAB is:

x=randn(n,1);
integral=mean(a*x.^3+b*x.^2+c*x+d);

since you can write the integral is $E[a X^3 + b X^2 + c X + d]$ where $ X \sim N(0,1)$.

octave:14> a=2,b=0,c=1,d=10
a =  2
b = 0
c =  1
d =  10
octave:15> n=1000
n =  1000
octave:16> x=randn(n,1);
octave:17> mean(a*x.^3+b*x.^2+c*x+d)
ans =  9.8828

Answer 2

You should use the Box–Muller transform to transform a pair of uniformly distributed random numbers from the interval $(0,1)$ to the standard normal distribution. Then you can simply average over the polynomial.