Monte Carlo integration not working?

I want to integrate (1/y)*(2/(1+(log(y))^2)) from 0 to 1. Wolfram alpha tells me that it must be pi. But when I do monte carlo integration in R, I keep getting 3.00 and 2.99 after trying more than 10 times. This is what I did:

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

I copied the exact function to wolfram alpha to verify that the integral must be pi

I tried to check if my y is correctly distributed by checking its value and building a historical diagram, and it seems that everything is in order. Is there something wrong with my computer?

Edit: Perhaps someone else can copy my code and run it on their own to make sure that it does not work on my computer.

+5
source share
1 answer

Ok, first we start with a simple transformation, log(x) -> x , making the integral

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

where S is the sign of integration.

Thus, the function 1 / (1 + x ^ 2) decreases monotonously and reasonably quickly. We need some reasonable PDF to select points in the interval [0 ... infinity], so that most of the area in which the original function is significant is covered. We will use an exponential distribution with some free parameter, which we will use to optimize the selection.

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

So, we have k * e -kx as a properly normalized PDF in the range [0 ... infinity]. The function for integration is (2/(1+x^2))*exp(k*x)/k . We know that the selection from the exponent is mainly -log(U(0,1)) , so the code for this is very simple

 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) 

The result is 3.145954 , for the result of seed 22345 is 3.135632 , for the result of seed 32345 is 3.146081 .

UPDATE

Returning to the original function [0 ... 1] is quite simple

UPDATE II

changed to prof.Bolker offer

+6
source

Source: https://habr.com/ru/post/1240654/


All Articles