问题描述
所以我有这个概率分布
X = {0概率7/8}
{1/60概率1/8}
James的汽车每年要发生N次故障,其中N〜Pois(2),X是维修成本,Y是James一年造成的总成本。
我想计算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
您需要进行大量仿真,但是上面应该是可以开始优化以使其更接近的代码。