Monte Carlo Stil einen integrierten MATLAB

Question

Ich wusste, zu bewerten, dass wir pi mit Monte-Carlo-Methode von „Werfen“ Punkt auf der oberen rechten Ecke nähern könnte und zählen, wie viele von ihnen innerhalb des Kreises sind etc ..

ich möchte, dass jede Funktion f für tun, also bin ich zufällige Punkte im Rechteck [a, b] x [0; max (f)] "werfen" und teste ich, ob mein random_point_y ist niedriger als f (random_point_x) und dann dividiere ich den Gesamtbetrag durch die Anzahl der Punkte unter f.
Hier ist der 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 

ich zum Beispiel für f hatte = e^(- x^2) zwischen -1 und 1:

Aber ich Folge haben 1,3414, 1.3373 für 500.000 Punkte. Das genaue Ergebnis ist 1.49365

Was fehlt mir?

Answer 1

Sie haben zwei kleine Fehler:

• Es count/n werden sollte, nicht n/count. Mit dem richtigen count/n wird der Anteil der Punkte unter der Kurve geben.
• Um den Bereich unter der Kurve zu erhalten, multiplizieren Sie diesen Anteil mit der Fläche des Rechtecks, (b-a)*max_value.

Verwenden Sie also count/n * (b-a)*max_value.

---

Abgesehen davon, würde der Code schneller und klarer mit einigen Vektorisierung:

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