Search code examples
matlabmontecarloapproximation

Monte Carlo style to evaluate an integral MATLAB

I knew that we could approximate pi using Monte Carlo method by "throwing" point on the top right corner and count how many of them are inside the circle etc..

I want to do that for every function f, so I am "throwing" random points in the rectangle [a,b] x [0;max(f)] and I'm testing if my random_point_y is lower than f(random_point_x) and then I divide the total amount by the number of point below f.
Here is the code : 

clear
close all
%Let's define our function f
clear
close all
f = @(x) exp(-x.^2);
a=-1; b=1;
range = [a:0.01:b];
f_range = f(range);

%Let's find the maximum value of f
max_value = f(1);
max_x = range(1);
for i = range
    if (f(i) > max_value) %If we have a new maximum
        max_value = f(i);
        max_x = i;
    end
end


n=5000;
count=0;
%Let's generate uniformly distributed points over [a,b] x [0;max_f]
x = (b-a)*rand(1,n) + a;
y = rand(1,n) * max_value;

for i=1:n
    if y(i)<f(x(i)) %If my point is below the function
        count = count + 1;
    end
end


%PLOT
hold on

%scatter(x,y,'.')
plot(range,f_range,'LineWidth',2)
axis([a-1 b+1 -1 max_value+1])
integral = (n/count)

hold off

for example I had for f = e^(-x^2) between -1 and 1 : monte carlo

But I have for result 1.3414,1.3373 for 500.000 points. The exact result is 1.49365

What am I missing ?

Solution

  • You have two small mistakes:

    • It should be count/n, not n/count. Using the correct count/n will give the proportion of points below the curve.
    • To get the area below the curve, multiply that proportion by the area of the rectangle, (b-a)*max_value.

    So, use count/n * (b-a)*max_value.

    Apart from this, your code would be faster and clearer with some vectorization:

    clear
close all
f = @(x) exp(-x.^2);
a=-1; b=1;
range = [a:0.01:b];

%Let's find the maximum value of f
max_value = max(f(range));

n=50000;
%Let's generate uniformly distributed points over [a,b] x [0;max_f]
x = (b-a)*rand(1,n) + a;
y = rand(1,n) * max_value;

count = sum(y < f(x));
result = count/n * (b-a)*max_value