用蒙特卡罗模拟计算方差期望值

Question

所以我有这个概率分布

X= {0 _=_

(1/60概率1/8}

他的车每年坏N次，其中N~ Pois(2)和X，修理费用和Y是James在一年内造成的总成本。

我要计算EY和V(Y)，这应该给我EX=15和V(Y) = 1800。

我有一个蒙特卡洛模拟：

expon_dis <- rexp(200, 1/60)

result_matrix2 <- rep(0, 200)
expected_matrix <- rep(0, runs)

for (u in 1:runs){
  expon_dis <- rexp(200, 1/60)
  N <- rpois(200, 2)
  for (l in 1:200){
    result_matrix2[l] <- (expon_dis[l] * (1/8)) * (N[l])
  }
  expected_matrix[u] <- mean(result_matrix2)
}

这个代码给出了15的期望值，但是方差是不正确的。那么这个模拟有什么问题呢？

Answer 1

没有足够的时间来阅读你的代码，但是我认为这个错误伴随着乘法。

下面是一个非常粗略的实现，首先编写一个函数来模拟成本，给定x个故障数：

sim_cost = function(x){
cost = rexp(x,1/60)
prob = sample(c(0,1/60),x,prob=c(7/8,1/8),replace=TRUE)
sum(cost[prob>0])
}

然后生成每年的故障数：

set.seed(111)
N <- rpois(500000, 2)

多年来迭代，如果0，我们返回0：

set.seed(111)
sim = sapply(N,function(i)if(i==0){0}else{sum(sim_cost(i))})

mean(sim)
[1] 14.98248
var(sim)
[1] 1797.549

您需要相当多的模拟，但上面应该是一个代码，您可以开始优化，以使它更接近。