MATLAB: Monte Carlo Integration in MATLAB, help

integrationmathematicsMATLAB

Use Monte Carlo Integration to evaluate the integral of f(x,y)=x*(y^2), over x(0,2) and y(0,x/2).
My code is below, however it generates an answer of roughly 0.3333, which is incorrect because the exact value is 0.2667. Please help in correcting my code.
samplesize = 1000;
fxy = @(x,y) x.*y.^2; %integrand
x = 2*rand(1, samplesize); %uniform x ~(0,2)
y = (x./2).* rand(1, samplesize); %uniform y ~(0,x/2)
m = 2; %measure of domain
Ef = (1/samplesize)*sum(fxy(x,y)); %expected value
integral_value = m*Ef %estimation of integral

Best Answer

What is the problem?
samplesize = 1000;
fxy = @(x,y) x.*y.^2; %integrand
x = 2*rand(1, samplesize); %uniform x ~(0,2)
y = (x./2).* rand(1, samplesize); %uniform y ~(0,x/2)
plot(x,y,'.')
Do you, for some reason, expect Monte Carlo to be exact? Only kidding. ;-)
Your problem is you are not actually sampling uniformly over that triangular domain.
That is, when x is large, you have a lower density of points in y, than when x is small.
Your sampling scheme is completely wrong. I can see what you were thinking, but, still wrong. You sampled over the proper region. But it was not uniform, so the Monte Carlo failed. How can we fix this?
A simple solution is to sample uniformly over the triangle. For example, you might have done it using a trick like this:
xy = sort(rand(2,samplesize),1,'descend');
x = 2*xy(1,:);
y = xy(2,:);
plot(x,y,'.')
Next, the area of that domain is 1, NOT 2. So m=1 is correct.
m = 1; %measure of domain
Ef = (1/samplesize)*sum(fxy(x,y)); %expected value
integral_value = m*Ef %estimation of integral
integral_value =
      0.27315
With a larger value for samplesize, so 1000000, I get
integral_value =
      0.26646
which seems quite reasonable.
syms x y
>> int(int(x.*y^2,y,[0,x/2]),[0,2])
ans =
4/15
>> vpa(ans)
ans =
0.26666666666666666666666666666667
Related Question