R:使用重要性采样的蒙特卡罗积分

Question

我有一个积分要评估

      "x^(-0.5)" ; x in [0.01,1] 

我正在使用重要性采样 MC :
该理论说必须使用近似的 PDF 来计算期望值(几乎肯定会收敛到积分的平均值)

在绘制给定的积分和指数 PDF 后，仅基于图，我选择了
rexp和 dexp生成 PDF - 我的代码看起来像这样 -

#Without Importance Sampling
set.seed(1909)
X <- runif(1000,0.01,1)
Y <- X^(-0.5)
c( mean(Y), var(Y) )

#Importance sampling Monte Carlo
w <- function(x) dunif(x, 0.01, 1)/dexp(x,rate=1.5)
f <- function(x) x^(-0.5)
X= rexp(1000,rate=1.5)
Y=w(X)*f(X)
c( mean(Y), var(Y) )

有人可以确认我的思路是否正确吗？
如果错了，我应该如何处理这个问题？
请阐明 - 我已经理解了该理论，但事实证明，实现对我来说是有问题的。

对于不那么简单的积分，

1.) f(x) = [1+sinh(2x)ln(x)]^-1
我选择了 正常 PDF = g(x) (均值 = 0.5 和 SD = 5)仅在观察绘图后作为近似值。我写了一个类似于它的代码，但它说 NaN 是在重要性采样的情况下产生的。 (这在理想情况下意味着未定义的函数，但我不知道如何解决这个问题)。

2.) f(x,y) = exp(-x^4 - y^4)

我该如何选择 g(x,y) 对于上述功能？

Answer 1

一般来说，您的方法似乎是正确的，但您必须对要集成的域更加小心。在您的原始示例中，大约 20% 的值 rexp(1000, 1.5)高于 1. 函数 dexp(x, rate=1.5)不是区间 [0,1] 上的密度函数。你必须除以 pexp(1, rate=1.5) .因此，对于重要性采样示例，我将执行以下操作:

#Importance sampling Monte Carlo
w <- function(x) dunif(x, 0.01, 1)/dexp(x,rate=1.5) * pexp(1, rate=1.5)
f <- function(x) x^(-0.5)
X <- rexp(1000,rate=1.5)
X <- X[X<=1]
Y <- w(X)*f(X)
c(mean(Y), var(Y))

在你的第二个例子中，同样的事情导致了问题。您得到负 X，因此得到 log(X) 的 NA 值。此外，您的正常函数应以 0.5 为中心，方差较小。这是我的方法:

#Without Importance Sampling
set.seed(1909)
X <- runif(1000,0.01,1)
Y <- (1+sinh(2*X)*log(X))^(-1)
c(mean(Y), var(Y))

#Importance sampling Monte Carlo
w <- function(x) dunif(x, 0.01, 1)/dnorm(x, mean=0.5, sd=0.25) * (1-2*pnorm(0, mean=0.5, sd=0.25))
f <- function(x) (1+sinh(2*x)*log(x))^(-1)
X <- rnorm(1000, mean=0.5, sd=0.25)
Y1 <- w(X)
Y2 <- f(X)
Y <- Y1*Y2
Y <- Y[!(is.na(Y2)&Y1==0)]
c(mean(Y), var(Y))

在你的第二个例子中，我不太明白 y是。它只是一个常数吗？那么也许威 bool 分布可能会起作用。

编辑:关于您在评论中的其他问题。
(1) 任何概率密度函数都应该积分为 1。因此 dexp(x, rate=1.5)不是区间 [0,1] 上的密度函数，它只积分到 pexp(1, rate=1.5) .然而，函数

dexp01 <- function(x, rate){
  dexp(x, rate=rate)/pexp(1, rate=rate)
}

实际上集成为 1:

integrate(dexp, 0, 1, rate=1.5)
integrate(dexp01, 0, 1, rate=1.5)

这就是包含概率分布函数的基本原理。如果您有不同的间隔，例如[0.3,8]，你必须相应地调整函数:

dexp0.3_8 <- function(x, rate){
  dexp(x, rate=rate)/(pexp(8, rate=rate)-pexp(0.3, rate=rate))
}
integrate(dexp0.3_8, 0.3, 8, rate=1.5)

(2) 在这里，我选择方差使得 rnorm(1000, .5, .25) 中大约 95% 的值位于 0 到 1 的区间内(在此区间之外有许多值肯定会增加方差)。但是，我不确定这是分布函数的最佳选择。重要性函数的选择是一个我不太熟悉的问题。你可以问 CrossValidated .你的下一个问题也是如此。