4
$\begingroup$

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

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

4 Answers 4

4
$\begingroup$

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

enter image description here

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

$\endgroup$
1
  • $\begingroup$ 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! $\endgroup$
    – J. Barbosa
    May 21, 2016 at 16:09
2
$\begingroup$

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);
$\endgroup$
2
$\begingroup$

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.
$\endgroup$
2
$\begingroup$

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);
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.