我有以下功能
f(x)∝|x| exp(-1/2 |x| )+1/(1+(x-40)^4 ),xϵR
我想通过辛普森的方法（数值积分）、标准蒙特卡洛方法、接受拒绝抽样、重要性抽样、Metropolis-Hasting 算法、吉布斯抽样和贝叶斯模型找出 E(X) 和 E(X^3) 使用 MCMC （我还没决定）。
如何验证从不同方法获得的结果？我试图用数学方法求解 E(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，绝对误差 < 0.004
例如，我从这些方法中获得的结果并不相似
(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() 使用自适应求积方法。这种方法与辛普森一家和重要性抽样不同吗？在比较这些方法时，我应该期望结果有什么相似之处