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

Answer 2

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);

Answer 3

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.

So to rewrite your code:

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.

Answer 4

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);