Monte Carlo Integration funktioniert nicht?

Question

Ich möchte (1/y)*(2/(1+(log(y))^2)) von 0 bis 1 integrieren. Wolfram alpha sagt mir, das sollte pi sein. Aber wenn ich monte carlo Integration in R mache, bekomme ich 3,00 und 2,99, nachdem ich 10 mal probiert habe. Dies ist, was ich getan habe:

y=runif(10^6) 
f=(1/y)*(2/(1+(log(y))^2)) 
mean(f) 

ich die genaue Funktion in Wolfram Alpha kopiert zu überprüfen, ob das Integral pi sein sollte

Ich versuchte, mein y zu überprüfen, ob richtig Mittelwert es verteilt wird, indem überprüft und Plotten ein historgram, und es scheint in Ordnung zu sein. Könnte etwas mit meinem Computer nicht stimmen?

Edit: Vielleicht könnte jemand anderes meinen Code kopieren und selbst ausführen, um zu bestätigen, dass es nicht mein Computer ist.

Answer 1

Ok, lassen Sie uns zuerst mit einfachen Transformation beginnen, log(x) -> x, so dass Integral

I = S 2/(1+x^2) dx, x in [0...infinity] 

wo S Integration Zeichen.

Also Funktion 1/(1 + x^2) fällt monoton und vernünftig schnell. Wir brauchen ein vernünftiges PDF, um Punkte im Intervall [0 ... unendlich] zu samplen, so dass der Großteil der Region, in der die ursprüngliche Funktion signifikant ist, abgedeckt ist. Wir werden die Exponentialverteilung mit einigen freien Parametern verwenden, die wir verwenden werden, um die Abtastung zu optimieren.

I = S 2/(1+x^2)*exp(k*x)/k k*exp(-k*x) dx, x in [0...infinity] 

So haben wir k * e ^-kx als richtig normalisiert PDF im Bereich von [0 ... infinity]. Zu integrierende Funktion ist (2/(1+x^2))*exp(k*x)/k. Wir wissen, dass Probenahme aus exponentiellen grundsätzlich -log(U(0,1)) ist, so Code zu tun, ist sehr einfach

k <- 0.05 

# exponential distribution sampling from uniform vector 
Fx <- function(x) { 
    -log(x)/k 
} 

# integrand 
Fy <- function(x) { 
    (2.0/(1.0 + x*x))*exp(k*x)/k 
} 

set.seed(12345) 

n <- 10^6L 
s <- runif(n) 

# one could use rexp() as well instead of Fx 
# x <- rexp(n, k) 
x <- Fx(s) 

f <- Fy(x) 

q <- mean(f) 

print(q) 

Ergebnis zu 3.145954 gleich ist, für Saatgut 22345 Ergebnis ist gleich 3.135632, für Saatgut 32345 Ergebnis 3.146081 gleich ist.

UPDATE

Gehen wir zurück zu den ursprünglichen Funktion [0 ... 1] ist ganz einfach

UPDATE II

pro prof.Bolker Vorschlag geändert