蒙特卡罗整合不起作用？

Question

我希望(1/y)*(2/(1+(log(y))^2))从0到1 进行整合.Wolfram alpha告诉我这应该是pi.但是当我在R中进行蒙特卡罗整合时,我在尝试超过10次后仍然获得3.00和2.99.这就是我所做的:

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

我将确切的函数复制到wolfram alpha中以检查积分应该是pi

我试图通过检查它的平均值并绘制历史图来检查我的y是否正确分布,它似乎没问题.我的电脑有问题吗？

编辑:也许其他人可以复制我的代码并自己运行,以确认它不是我的电脑表现.

Answer 1

好的,首先让我们从简单的转换开始log(x) -> x,制作积分

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

S集成标志在哪里.

因此,函数1 /(1 + x ^ 2)单调且合理地快速下降.我们需要一些合理的PDF来对[0 ...无穷大]区间中的点进行采样,以便涵盖原始函数重要的大部分区域.我们将使用指数分布和一些我们将用于优化采样的自由参数.

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

因此,我们将k*e ^-kx设置为[0 ... infinity]范围内的正确标准化PDF.整合的功能是(2/(1+x^2))*exp(k*x)/k.我们知道从指数采样基本上是-log(U(0,1))这样的,所以这样做的代码非常简单

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)

结果等于3.145954,对于种子22345结果等于3.135632,对于种子32345结果等于3.146081.

UPDATE

回到原始函数[0 ... 1]非常简单

更新II

每个prof.Bolker建议改变