Monte Carlo Approximation of $\int_0^1\int_0^x x^2y dy dx$

Question

I am trying to evaluate the integral $I=\int_0^1\int_0^x x^2y \cdot dydx$ using the Monte Carlo approach. Therefore I generate $n$ random x values for $x_i\in (0, 1)$ and then generate $n$ random y values for $y_i\in (0, x_i)$. I attempted this in matlab as follows:

n=100000;
s=0;

for i=1:n 
x=rand(1,1); 
y=rand(1,1) * x;
s = s + x^2 * y; 
end 

I = s/n;

This gives $I\approx 0.125...$ however, the integral can be easily evaluated analytically to obtain $I= 0.1$. I'm really at a loss as for where the overestimation is coming from.

Answer 1

If you want to find $\displaystyle\int^b_a f(y)\; dy$ by Monte Carlo integration then you want $(b-a)$ times the average of the sample values of $f(y)$ with $y$ drawn uniformly from that interval

So $\displaystyle \int_0^1\int_0^x x^2y \; dy \;dx = \int_0^1 x^2 \left(\int_0^x y \; dy\right) \; dx$ could be modelled by adapting your code to

n=100000;
s=0;

for i=1:n 
x=rand(1,1); 
y=rand(1,1) * x;
s = s + (1-0) * x^2 * (x-0) * y; 
end 

I = s/n;

Answer 2

Your region of integration is a right triangle that is half of the square $[0,1]\times [0,1]$ below the line $y=x$. As such you should have $y \in [0,x]$. However, you should not take $y=rand*x$ as that will have points much more densely packed for small $x$. You should just throw random numbers for each of $x$ and $y$ then throw out any where $y \gt x$.