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

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

• I guess that's because the probability law for your $x$ and $y$ is not uniform over the triangle, and should be. May 21, 2016 at 14:36

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$$

• Thanks for your answer! That should do the job , however I figured out what to do to evaluated double integrals over non-rectangular regions from the definition of this method. Basically instead of summing $f(x,y)$ we have to sum $f(x,y)*(y-0)$ in every iteration and after the cycle is done we just divide this sum by $N$ and it should work! May 21, 2016 at 16:09

First of all I want to make it clear that the answer Conrad Turner gave is the actual answer to the question I made and I'm posting this just to give further clarification on the topic.

I found out what to do in these cases where we have a non-rectangular region of integration and it follows from the definition of the Monte Carlo method: $$\int_{a}^{b}f(x)dx\approx(b-a)\frac{1}{N}\sum_{i}f(x_{i})$$ where $f(x_{i})$ is a random evaluation of the function in the interval $[a,b]$ and $N$ the number of points to consider.

Extending this definition to two dimension and generalizing to non-rectangular regions we have: $$\int_{c}^{d}\int_{a(y)}^{b(y)}f(x,y)dxdy\approx\int_c^d(b(y)-a(y))\frac{1}{N}\sum_if(x_i,y)dy\approx(b(y)-a(y))(d-c)\frac{N}{N^2}\sum_{i,j}f(x_i,y_j)$$

Using this definition the code should look like this:

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)*b(y);
end

output=output/N;
fprintf('Value : %9.8f\n',output);


When using Matlab try and avoid for-loops as they often are extremely much slower than vectorized code. So instead generate two vectors, make the function f able to handle vectors and then do a sum on the samples-index. The difference in speed can easily be at least hundreds of times.

f=@(x,y) x.^2.*y;    % dot is entrywise or Hadamard or Schur
% operation on the in-vectors x and y
N=50000;
y=c+(d-c)*rand(N,1);     % row vector of N elements of random numbers.
x=a+(b(y)-a)*rand(N,1);  % another one.
output=sum(f(x,y));      % apply the vector function and sum it up.


As noted, the way you pick $(x,y)$ is not uniform over the triangle $T:=\{(x,y)\in[0,1]^2,y\geq x\}$. To see that, note that the probability for $y$ to be less than $1/2$ is $1/2$ in your algorithm, where it should be $$\frac{|\{(x,y)\in T, y\leq 1/2\}|}{|T|}=\frac{1}{4}.$$ Let's be more specific, and explicit the law of $(x,y)$ in your algorithm. For a given $y$, $x$ is picked uniformly on $[0,y]$. Thus, our law $\mu$ has density $d\mu(x,y)=\frac{1}{y}dxdy$. The Monte-Carlo integration method gives, according to the law of large number, that $$\int_T f(x,y)d\mu(x,y)=\mu(T)\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^Nf(x_i,y_i)$$ where $((x_i,y_i))_{i\in\mathbb N^*}$ are random variables independent and identically distributed with probability distribution $\mu$. Now, it is easy to check that $\mu(T)=1$, so your algorithm actually gives $$0.5*\sum_{i=1}^{50000}f(x_i,y_i)\approx\frac{1}{2}\int_T x^2yd\mu(x,y)=\frac{1}{2}\int_T x^2dxdy=\frac{1}{24}.$$

To correct this, you need to find a way of picking $(x,y)$ uniformly on $T$. Conrad's way, for example, works fine, or you can simply adjust your algorithm this way by integrating $f\chi_T$ on $[0,1]^2$ instead, where it is easier to draw random variables ($\chi$ being the characteristic function):

output=0;
a=0;
b=1;
c=0;
d=1;

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

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

output=(b-a)*(d-c)*output/N;
fprintf('Value : %9.8f\n',output);