通过 R 进行数值积分的问题

Question

我有以下功能

f(x)∝|x| exp(-1/2 |x| )+1/(1+(x-40)^4 ),xϵR

我想通过辛普森方法（数值积分）、标准蒙特卡罗方法、接受拒绝抽样、重要性抽样、Metropolis-Hasting 算法、吉布斯抽样和贝叶斯算法找出 E(X) 和 E(X^3)使用 MCMC 的模型（我还没有决定）。

如何验证从不同方法获得的结果？ 我试图用数学方法求解 E(X) 但找不到任何接近的形式。此功能可以分为不同的部分，如

absolute(x)*双指数密度 + 另一个利用 X 的更高幂 (4) 的反函数。 由于绝对 (x) 和范围 [-Inf, Inf] 我们总是必须将它划分为 [-Inf, 0] 和 [0, Inf]。通过按部分积分，我能够将第一部分视为（无限范围内的绝对（x）+（x^2/2））+这部分的积分无法在数学上找到。

所以我利用下面的代码得到数值积分结果为

Library(stats)
integrand <- function(x) {x*(abs(x)* exp(-0.5*abs(x))+(1/(1+(x-40)^4)))}
integrate(integrand, lower = -Inf, upper = Inf)

因此结果是 E(X)= 88.85766，绝对误差

例如，我从这些方法中获得的结果并不相似

(i) 通过辛普森方法我得到 E(X) = 0.3222642 和 E(X^3)=677.0711..

simpson_v2 <- function(fun, a, b, n=100) {
    # numerical integral using Simpson's rule
    # assume a < b and n is an even positive integer
    if (a == -Inf & b == Inf) {
        f <- function(t) (fun((1-t)/t) + fun((t-1)/t))/t^2
        s <- simpson_v2(f, 0, 1, n)
    } else if (a == -Inf & b != Inf) {
        f <- function(t) fun(b-(1-t)/t)/t^2
        s <- simpson_v2(f, 0, 1, n)
    } else if (a != -Inf & b == Inf) {
        f <- function(t) fun(a+(1-t)/t)/t^2
        s <- simpson_v2(f, 0, 1, n)
    } else {
        h <- (b-a)/n
        x <- seq(a, b, by=h)
        y <- fun(x)
        y[is.nan(y)]=0
        s <- y[1] + y[n+1] + 2*sum(y[seq(2,n,by=2)]) + 4 *sum(y[seq(3,n-1, by=2)])
        s <- s*h/3
    }
    return(s)
}

EX  <- function(x) x*(abs(x)* exp(-0.5*abs(x))+(1/(1+(x-40)^4)))
simpson_v2(EX, -Inf, Inf, n=100)

EX3 <- function(x) (x^3)*(abs(x)* exp(-0.5*abs(x))+(1/(1+(x-40)^4)))
simpson_v2(EX3, -Inf, Inf, n=100) 

(ii) 重要性抽样 我的提案密度是正常的，平均值=0，标准差=4。我应用的重要性抽样过程总结如下

假设我不能从 f(x) 中采样，这是真的，因为它没有众所周知的形式，并且 R 中没有内置函数可用于采样。因此，我建议使用另一个对数凹尾分布 N(0, 4) 来采样，这样我就不会估计 E(x)，而是估计 E(x*f(x)/N(0,1))。我为此使用以下代码，它从 N(0,4) 中获取 100000 个样本

X <- rnorm(1e5, sd=4)
Y <- (X)*(abs(X)*exp(-0.5*abs(X))+(1/(1+(x-40)^4)))/(dnorm(X, sd=4))
mean(Y)

由于此代码需要从正态分布中随机抽样，因此每次我得到不同的答案，但它大约是 -0.1710694，与 0.3222642 几乎相似。我是从辛普森的方法中得到的。但是这些结果与积分（）有很大的不同 E(X)= 88.85766。请注意，integrate() 使用自适应求积方法。这种方法与辛普森一家和重要性抽样不同吗？在比较这些方法时，我应该期望结果有什么相似性

Answer 1

首先，EX 和 EX3 的定义是错误的，你错过了指数下的减号

好吧，这里有一些简化

• 如果你整合这部分x*(abs(x)* exp(-0.5*abs(x)))，从-infinity...infinity 结果将是0

• 如果你整合这部分x^3*(abs(x)* exp(-0.5*abs(x)))，从-infinity...infinity 结果将是0

• 积分 x^3/(1+(x-40)^4) 从 -infinity...infinity 将是无穷大，我敢冒险，你会得到无穷大的对数，请参阅 http://integrals.wolfram.com/index.jsp?expr=%28xxx% 29%2F%281+%2B+%28x-a%29%5E4%29&random=false

• 积分 x/(1+(x-40)^4) 看起来类似于反正切，尽管在线积分器提供了丑陋的输出 http://integrals.wolfram.com/index.jsp?expr=x%2F%281+%2B+%28x-a%29%5E4%29&random=false

更新

看起来你的EX 应该是40*\pi / \sqrt{2}

而EX3 不是无穷大，我这里可能错了

更新 2

是的，EX3 是有限的，应该是 a^2*EX + \pi*a*3/\sqrt{2}，其中a 等于 40

更新 3

如前所述，还需要规范化才能获得 EX 和 EX3 的真实值

N = 8 + \pi/\sqrt{2}

计算积分除以N 以获得正确的矩。