Monte Carlo 计算积分 MATLAB 的样式

Question

我知道我们可以通过右上角的 "throwing" 点使用 Monte Carlo 方法近似 pi 并计算其中有多少个在圈子等..

我想对 每个函数 f 都这样做，所以我在 "throwing" 矩形 [=] 中的随机点 23=][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 之间：

但我有 500.000 点的结果 1.3414,1.3373。确切的结果是 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

Monte Carlo 计算积分 MATLAB 的样式

Monte Carlo style to evaluate an integral MATLAB

matlab

montecarlo

approximation