Monte carlo integration not working?

Question

I wish to integrate (1/y)*(2/(1+(log(y))^2)) from 0 to 1. Wolfram alpha tells me this should be pi. But when I do monte carlo integration in R, I keep getting 3.00 and 2.99 after trying over 10 times. This is what I have done:

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

I copied the exact function into wolfram alpha to check that the integral should be pi

I tried to check if my y is properly distributed by checking it's mean and plotting a historgram, and it seems to be ok. Could there be something wrong with my computer?

Edit: Maybe someone else could copy my code and run it themselves, to confirm that it isn't my computer acting up.

Answer 1

Ok, first let's start with simple transformation, log(x) -> x, making integral

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

where S is integration sign.

So function 1/(1+x^2) is falling monotonically and reasonable fast. We need some reasonable PDF to sample points in [0...infinity] interval, such that most of the region where original function is significant is covered. We will use exponential distribution with some free parameter which we will use to optimize sampling.

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

So, we have k*e^-kx as properly normalized PDF in the range of [0...infinity]. Function to integrate is (2/(1+x^2))*exp(k*x)/k. We know that sampling from exponential is basically -log(U(0,1)), so code to do that 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)

Result is equal to 3.145954, for seed 22345 result is equal to 3.135632, for seed 32345 result is equal to 3.146081.

UPDATE

Going back to original function [0...1] is quite simple

UPDATE II

changed per prof.Bolker suggestion