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

【问题标题】：Monte-Carlo Simulation of seasonal and rare order quantities季节性和稀有订单数量的蒙特卡罗模拟
【发布时间】：2023-04-07 04:33:01
【问题描述】：

给定订单数量的季节性需求模式和可能订单数量的历史样本。如何执行 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 的形状应该与季节性模式匹配。

【问题讨论】：

标签： r montecarlo

【解决方案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")

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

【讨论】：