[Math] Monte Carlo double integral over a non-rectangular region (Matlab)

Question

I want to evaluated the following integral using Monte Carlo method: $$\int_{0}^{1}\int_{0}^{y}x^2y\ dxdy $$
What I tried using Matlab:

output=0;
a=0;
b=@(y) y;
c=0;
d=1;

f=@(x,y) x^2*y;
N=50000;

for i=1:N
  y=c+(d-c)*rand();
  x=a+(b(y)-a)*rand();
  output=output+f(x,y);
end

output=0.5*output/N; % 0.5 because it's the area of the region of integration
fprintf('Value : %9.8f\n',output);

However this code didn't give me the value I was looking for ( the exact answer is $1/15$) and I have no idea why this won't work. Anyone knows what the problem is and what to fix? Cheers

Answer 1

Here is a function that will give you a vector of points uniformly distributed over the required triangular region:

function [rv]=RandTri(m)
% function to generate an output array of M x,y coordinates uniformly
% distributed over the upper triangle within the unit square

  rv=rand(2,m);
  rv1t=rv(1,:)+(rv(2,:)<rv(1,:)).*(-rv(1,:)+rv(2,:));
  rv(2,:)=rv(2,:)+(rv(2,:)<rv(1,:)).*(-rv(2,:)+rv(1,:));
  rv(1,:)=rv1t;
endfunction

And here is the plot generated by:

>> xy=RandTri(100);
>> plot(xy(1,:),xy(2,:),'o')

Now taking the average of your function over this sample provides an estimate of $$ I=\int_{y=0}^1\int_{x=0}^y x^2y\; p(x,y)\;dx\; dy=2\int_{y=0}^1\int_{x=0}^y x^2y\;dx\;dy $$