蒙特卡洛风格评估积分 MATLAB

Question

我知道我们可以使用蒙特卡洛方法来近似 pi，方法是在右上角“抛出”点并计算其中有多少在圆圈内等。

我想对 每个函数 f 都这样做，所以我在矩形 [a,b] x [0; 中“抛出”随机点 max(f)] 并且我正在测试我的 random_point_y 是否低于 f(random_point_x)，然后我将总数除以 f 下方的点数。
这是代码：

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

例如我有 f = e^(-x^2) 在 -1 和 1 之间：

但我的结果 1.3414,1.3373 获得 500.000 分。 确切的结果是 1.49365

我错过了什么？

Answer 1

你有两个小错误：

• 应该是count/n，而不是n/count。使用正确的count/n 将给出曲线下方点的比例。
• 要获得曲线下方的面积，请将该比例乘以矩形的面积，(b-a)*max_value。

所以，请使用count/n * (b-a)*max_value。

---

除此之外，通过一些矢量化，您的代码会更快、更清晰：

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