蒙特卡洛集成的答案不正确

Question

如何用R中的蒙特卡罗方法近似[x ^ 4 * sin（x）] / [exp（1）^（x / 5）]（0到+ Inf）的积分？

我试图做的是

set.seed(666)
func1 <- function(x)
    {
      (x^4 * sin(x))/exp(1)^(x/5)
    }

n <- 1000000

x <- rexp(n, 0.2)
f <- func1(x)
E <- mean(f)

但是E的结果是不对的

Answer 1

如果您要从指数进行采样，则不应在函数中再次使用它。

从代码

set.seed(32345)

func <- function(x) { (x^4 * sin(x)) }

n <- 10000000

x <- rexp(n, 0.2)
f <- func(x)
E <- mean(f)

我得到了答案

[1] 13.06643

UPDATE

它波动，并且波动很大。

至少首先从正确的答案开始，根据Mathematica等于4453125/371293 = 11.9936。

我改变了积分

I =∫dxexp（-x / 5）x4 sin（x）

使用替换y=x/5来

I = 55 C（5）∫dyexp（-y）y5-1 / c（5）sin（5 * y）

除了sin(5*y)之外的所有东西都是归一化的伽马分布，我们将用它来进行采样，而sin(5*y)将是我们计算平均值的函数。

并使用以下技巧和大量样本：我分析正值和负值的计算。如果你有相互抵消的价值波动的答案，它会有所帮助。我也分批进行了计算。 Gamma函数5只有4！（阶乘）

码

set.seed(32345)

N  <- 10000000 # number of samples per batch
NN <- 640      # number of batches

pos <- rep(0, NN) # positive values
neg <- rep(0, NN) # negative values

for(k in 1:NN) { # loop over batches
    y <- rgamma(N, shape=5, scale=1)
    f <- sin(5.0 * y)
    pnf <- ifelse(f > 0.0, f, 0.0)
    pos[k] <- mean(pnf)
    pnf <- ifelse(f < 0.0, -f, 0.0)
    neg[k] <- mean(pnf)
    print(k)
}

mean(pos)
sd(pos)/sqrt(NN)

mean(neg)
sd(neg)/sqrt(NN)

5*5*5*5*5*4*3*2*(mean(pos) - mean(neg))

产量

> mean(pos)
[1] 0.3183912
> sd(pos)/sqrt(NN)
[1] 4.749269e-06
> 
> mean(neg)
[1] 0.3182223
> sd(neg)/sqrt(NN)
[1] 5.087734e-06
> 
> 5*5*5*5*5*4*3*2*(mean(pos) - mean(neg))
[1] 12.67078

您可以看到我们确实计算了两个非常接近的值的差异，这就是为什么难以获得收敛的原因。在我的Xeon工作站上计算需要20多分钟。

并与不同的种子= 12345

> mean(pos)
[1] 0.3183917
> sd(pos)/sqrt(NN)
[1] 4.835424e-06
> 
> mean(neg)
[1] 0.3182268
> sd(neg)/sqrt(NN)
[1] 4.633129e-06
> 
> 5*5*5*5*5*4*3*2*(mean(pos) - mean(neg))
[1] 12.36735

Answer 2

在下面我故意不设置随机种子。

正如我在评论中提到的，在Stack Overflow上至少有两个关于Monte Carlo集成的问答：

两者都解释了如何获得蒙特卡洛估计，但忘记了估计的标准误差。事实证明，蒙特卡罗估计的函数收敛速度极慢。

众所周知，蒙特卡洛积分具有O(1 / sqrt(N))收敛速率，其中N是样本大小，O()是大O符号。然而，对于某些函数，大O后面的常数可能非常大，因此实际收敛速度可能要慢得多。

您的功能至少可以通过两种方式定义：

## direct definition
f <- function (x) x^4 * sin(x) * exp(-x/5)

## using gamma distribution; see ?rgamma
g <- function (x) sin(x) * 5 ^ 5 * gamma(5) * dgamma(x, 5, 1/5)

curve(f, from = 0, to = 100)
curve(g, add = TRUE, col = 2)

第一个问答解释了如何使用均匀分布的样本计算蒙特卡洛积分。你的函数f或g几乎为零，超过x = 200，因此[0, +Inf)上的积分在[0, 200]上有效。以下函数将返回集成及其标准错误：

MCI1 <- function (n) {
  x <- runif(n, 0, 200)
  y <- 200 * f(x)
  c(mean.default(y), sqrt(var(y) / n))
  }

另一种方法是通过重要性抽样，如第二个问答所述。这里的伽马分布用作提案分布（正如Ben Bolker建议的那样）。

MCI2 <- function (n) {
  x <- rgamma(n, 5, 0.2)
  y <- sin(x) * 75000
  c(mean.default(y), sqrt(var(y) / n))
  }

现在让我们检查收敛速度。

n <- seq(1000, by = 5000, length = 100)
tail(n)
#[1] 471000 476000 481000 486000 491000 496000

b1 <- sapply(n, MCI1)
b2 <- sapply(n, MCI2)

对于统一采样，我们有

par(mfrow = c(1, 2))
plot(b1[1, ], main = "estimate")
plot(b1[2, ], main = "standard error")

b1[, (ncol(b1) - 5):ncol(b1)]
#         [,1]     [,2]      [,3]      [,4]      [,5]      [,6]
#[1,] 115.1243 239.9631  55.57149 -325.8631 -140.3745  78.61126
#[2,] 181.0025 179.9988 178.99367  178.2152  177.2193 175.31446

对于伽马采样，我们有

par(mfrow = c(1, 2))
plot(b2[1, ], main = "estimate")
plot(b2[2, ], main = "standard error")

b2[, (ncol(b2) - 5):ncol(b2)]
#           [,1]       [,2]     [,3]      [,4]      [,5]      [,6]
#[1,] -100.70344 -150.71536 24.40841 -49.58032 169.85385 122.81731
#[2,]   77.22445   76.85013 76.53198  76.03692  75.69819  75.25755

无论采用何种方法，都要注意标准误差有多大（与估计本身相比），以及减少的速度有多慢。

使用数值积分要容易得多（集成单变量函数并不奇怪）：

integrate(f, 0, 200)
#11.99356 with absolute error < 0.0012

## trapezoidal rule
200 * mean.default(f(seq(0, 200, length = 10000)))
#[1] 11.99236

在梯形规则中，即使仅采用1e + 4个均匀间隔的采样点，积分也足够接近事实。

备注

如果我们在更受限制的域上进行集成，那么蒙特卡罗集成将会有更少的困难。从f或g的图中，我们看到这是一个振荡函数。实际上，它穿过x轴，周期为pi。让我们考虑在[lower, upper]上进行整合。

MCI3 <- function (n, lower, upper) {
  x <- runif(n, lower, upper)
  y <- (upper - lower) * f(x)
  c(mean.default(y), sqrt(var(y) / n))
  }

a1 <- sapply(n, MCI3, lower = 0, upper = pi)
a2 <- sapply(n, MCI3, lower = pi, upper = 2 * pi)
a3 <- sapply(n, MCI3, lower = 2 * pi, upper = 3 * pi)
a4 <- sapply(n, MCI3, lower = 3 * pi, upper = 4 * pi)

a1[, (ncol(a1) - 5):ncol(a1)]
#            [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
#[1,] 17.04658711 16.97935808 17.01094302 17.02117843 16.96935285 16.99552898
#[2,]  0.02407643  0.02390894  0.02379678  0.02368683  0.02354298  0.02342799

a2[, (ncol(a2) - 5):ncol(a2)]
#             [,1]         [,2]         [,3]         [,4]         [,5]
#[1,] -406.5646843 -404.9633321 -405.4300941 -405.4799659 -405.8337416
#[2,]    0.3476975    0.3463621    0.3442497    0.3425202    0.3409073
#             [,6]
#[1,] -405.8628741
#[2,]    0.3390045

a3[, (ncol(a3) - 5):ncol(a3)]
#            [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
#[1,] 1591.539911 1592.280780 1594.307951 1591.375340 1593.171500 1591.648529
#[2,]    1.197469    1.190251    1.183095    1.177079    1.172049    1.165667

a4[, (ncol(a4) - 5):ncol(a4)]
#             [,1]         [,2]         [,3]         [,4]         [,5]
#[1,] -3235.561677 -3239.147235 -3241.532097 -3238.421556 -3238.667702
#[2,]     2.336684     2.321283     2.311647     2.300856     2.286624
#             [,6]
#[1,] -3237.043068
#[2,]     2.279032

蒙特卡洛集成的答案不正确

问题描述投票：3回答：2

2个回答

最新问题

蒙特卡洛集成的答案不正确

问题描述 投票：3回答：2

2个回答

最新问题

问题描述投票：3回答：2