蒙特卡洛风格来评估整体MATLAB

我知道我们可以利用蒙特卡罗方法，通过“扔”点上右上角近似PI并计算其中有多少是圆等内部..蒙特卡洛风格来评估整体MATLAB

我想这样做，对于每函数f，所以我在矩形“扔” 随机点[A，b]×[0;最大（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

我错过了什么？

来源

2017-11-10 Romain B.

顺便说一句，你可以这样做：'一= -1;''B = 1;'' F = @（x）exp（-x。^ 2）;' 'n = 5000;' 'randnums = a +（ba）* rand（1，n）;' 'intg =（ba）* mean randnums））' –

是的，它的工作原理，但我真的想实施“解雇”。而且如果我设置'f = @（x）exp（-x。^ 2）;'和测试为'if x（i）^ 2 + y（i）^ 2 <= 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

来源

2017-11-10 22:27:32

@RomainB。那么，它将取代循环中的count = count + 1。比较'y

蒙特卡洛风格来评估整体MATLAB

回答

相关问题