Monte Carlo Integration in MATLAB, help?

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

Answer 1

[John D'Errico]

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

[Ron Burgundy]

Perfect answer. My initial sampling scheme made sense to me logically, but clearly was not uniform. Thank you very much John, and nice trick with the 'sort'.

[John D'Errico]

Until I plotted the points you were generating, it was not obvious what the problem was. What you were doing seemed to make sense at first glance. Once they were plotted, it became clear of course. The point being that you should always plot everything. I have found this to be perhaps the most useful piece of advice I can ever offer.

Thre are other ways to generate a uniform sample in a triangular region, but that was an easy one.