r - 季节性和稀有订单数量的蒙特卡罗模拟

Question

给定订单数量的季节性需求模式和可能的订单数量的历史样本。如何执行 MCS？

# months 1:12
mon <- 1:12

# seasonal probability density of rare orders including the material "A"
prob <- sin( -(mon + 2.5) * 2 * pi / 12)
prob <- prob - min(prob) + 0.5
prob <- prob / (sum(prob))
plot(mon, prob, type = "b", main = "seasonality and density of probabilities", xlab = "months", ylab = "probability", ylim = c(0, .2))

# historical order quantities except 0 for "A"
quantity_A <- c(15, 3.4, 3.4) # there are 3 observations, in the other months quantiy should be 0
paste("expected average of simulated 12 months=", round(sum(quantity_A) / 12, 4))

想要：多个“历史上”（=数量可以是：15、3.4 或 0）模拟年（=12 个月），其数量包络类似于上面的季节性概率。

示例：0 0 0 0 15 0 3.4 0 0 0 0。因此，如果我多次运行模拟，形状quantity_A应该与季节性模式相匹配。

标签： rmontecarlo

Answer 1

最终这就是我所做的。它旨在模拟每 12 个月不超过 6 个非零订单数量的订单数量。

set.seed(76543)

# months 1:12
mon <- 1:12

# seasonal probability density of rare orders including the material "A"
prob <- sin( -(mon + 2.5) * 2 * pi / 12)
prob <- prob - min(prob) + 0.5
prob <- prob / (sum(prob))
plot(mon, prob, type = "b", main = "seasonality and density of probabilities", xlab = "months", ylab = "probability", ylim = c(0, .2))

# historical order quantities except 0 for "A"
quantity_A <- c(15, 3.4, 3.4) # there are 3 observations, in the other months quantiy should be 0
paste("expected average of simulated 12 months=", round(sum(quantity_A) / 12, 4))
# expected average of simulated 12 months= 1.8167

# sample 1 year = 12 months
sample_year <- function(probs, qtys) {
  res <- sample(qtys, 12, replace = TRUE) # sample historical order quantities ..
  res[(runif(12) / length(qtys)) > probs] <- 0 # .. but only accept, if greater than PDF
  res
}

sims <- t(replicate(100000, sample_year(prob, quantity_A)))

barplot(colSums(sims))
# there is definitely a peak in the high season

paste("simulated average = ", round(mean(rowMeans(sims)), 4))
# simulated average =  1.8166

colSums(sims)/sum(sims)

plot(cumsum(prob), cumsum(colSums(sims)/sum(sims)), main = "simulated vs target probabilities", type = "b")

但我绝对仍然有兴趣获得对我的问题的答复。