Monte carlo 集成不起作用?
Monte carlo integration not working?
我希望将 (1/y)*(2/(1+(log(y))^2))
从 0 积分到 1。Wolfram alpha 告诉我这应该是 pi。但是当我在 R 中进行 monte carlo 集成时,我在尝试 10 次以上后不断得到 3.00 和 2.99。这就是我所做的:
y=runif(10^6)
f=(1/y)*(2/(1+(log(y))^2))
mean(f)
我将确切的函数复制到 wolfram alpha 中以检查积分应该是 pi
我试图通过检查它的均值和绘制直方图来检查我的 y 是否正确分布,它似乎没问题。会不会是我的电脑出问题了?
编辑:也许其他人可以复制我的代码并 运行 自己复制,以确认不是我的电脑出了问题。
好,先从简单的变换开始,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] 非常简单
更新二
根据 prof.Bolker 建议更改
我希望将 (1/y)*(2/(1+(log(y))^2))
从 0 积分到 1。Wolfram alpha 告诉我这应该是 pi。但是当我在 R 中进行 monte carlo 集成时,我在尝试 10 次以上后不断得到 3.00 和 2.99。这就是我所做的:
y=runif(10^6)
f=(1/y)*(2/(1+(log(y))^2))
mean(f)
我将确切的函数复制到 wolfram alpha 中以检查积分应该是 pi
我试图通过检查它的均值和绘制直方图来检查我的 y 是否正确分布,它似乎没问题。会不会是我的电脑出问题了?
编辑:也许其他人可以复制我的代码并 运行 自己复制,以确认不是我的电脑出了问题。
好,先从简单的变换开始,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] 非常简单
更新二
根据 prof.Bolker 建议更改