r - 蒙特卡洛积分不起作用？

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...infinity] 区间内的点进行采样，以便覆盖原始函数显着的大部分区域。我们将使用带有一些自由参数的指数分布来优化采样。

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 .

更新

回到原来的函数 [0...1] 很简单

更新二

根据 Bolker 教授的建议进行更改