R:使用重要性采样的蒙特卡洛集成

我有一个整体来评估

      "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)？

最佳答案

通常，您的方法似乎是正确的，但是对于要集成的域，您必须格外小心。在您的原始示例中，rexp(1000, 1.5)的值大约有20％大于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。下一个问题也一样。

关于R:使用重要性采样的蒙特卡洛集成，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/22060675/