蒙特卡洛風格來評估整體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

回答

相關問題