所以我有这个概率分布
X = {0 概率 7/8}
{1/60 概率 1/8}
詹姆斯的车每年坏掉 N 次,其中 N ~ Pois(2),X 是维修费用,Y 是詹姆斯一年内造成的总费用。
我想计算 E[Y] 和 V(Y),这应该给我 E[X]=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,但方差不正确。那么这个模拟有什么问题呢?
最佳答案
没有足够的时间阅读您的代码,但我认为错误来自乘法。
下面是一个非常粗略的实现,首先编写一个函数来模拟成本,给定 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
您需要相当多的模拟,但上面应该是您可以开始优化以使其更接近的代码。
关于r - 使用蒙特卡罗模拟计算方差的期望值,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/63967198/