蒙特卡罗积分无法工作？

Question

蒙特卡罗积分无法工作？

5

我希望集成(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是否正确分布，看起来没问题。我的电脑可能有问题吗？

编辑：也许其他人可以复制我的代码并自行运行，以确认这不是我的电脑出了问题。

- user124249

1

说到10次，你是指10个样本吗？那太少了！应该做大约100万个样本。同时，还要考虑重要性采样（而不是均匀随机数分布），这有助于减少误差。 - a_guest

我做了一百万个样本，没错。我的意思是我运行了相同的代码10次。 - user124249

这不是使用mc方法的最佳积分示例，因为在零点处有渐近线 - 你会期望它低估。 - user20650

1

@Roland 当然可以。 - Severin Pappadeux

1

@Roland，请阅读答案。 - Severin Pappadeux

显示剩余8条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Severin Pappadeux · Accepted Answer

好的，首先让我们从简单的转换开始，log(x) -> x，将其积分。

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

其中S是积分符号。

因此，函数1/(1+x^2)单调递减且下降得相当快。我们需要一些合理的概率密度函数来在[0...无穷大]区间内采样点，以便涵盖原始函数显著的大部分区域。我们将使用指数分布，并设置一些自由参数来优化采样。

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

因此，在[0...无穷大]范围内，我们将k*e^-kx作为适当归一化的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]非常简单。

更新 II

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