First we’ll setup the Monte Carlo’s table

This is a 100 days observation, so if we sum the frequency its equal to 100

demand <- c(50, 60, 70, 80, 90)
frequency  <- c(10, 20, 40, 20, 10)
probability <- frequency / sum(frequency)

cat("Probability Distribution:\n")
## Probability Distribution:
print(data.frame(demand, probability))
##   demand probability
## 1     50         0.1
## 2     60         0.2
## 3     70         0.4
## 4     80         0.2
## 5     90         0.1

Forecast for 5 Days & 20 Days

monte_carlo <- function(n_day, seed = 42) {
  set.seed(seed)
  result <- sample(demand, size = n_day, replace = TRUE, prob = probability)
  return(result)
}

sim5 <- monte_carlo(5)
cat("\n--- 5 Days Forecast---\n")
## 
## --- 5 Days Forecast---
cat("Demand everyday:", sim5, "\n")
## Demand everyday: 50 50 70 90 80
cat("Average demand:", mean(sim5), "\n")
## Average demand: 68
sim20 <- monte_carlo(20)
cat("\n--- 20 Days Forecast ---\n")
## 
## --- 20 Days Forecast ---
cat("Demand everyday:", sim20, "\n")
## Demand everyday: 50 50 70 90 80 60 80 70 80 80 60 80 50 70 60 50 50 70 60 60
cat("Average demand:", mean(sim20), "\n")
## Average demand: 66
ev <- sum(demand * probability)
cat("\nExpectation Value (E[X]):", ev, "\n")
## 
## Expectation Value (E[X]): 70

Monte Carlo is forecasts is based on the probability. If we forecast 100 days, the appearance of some value will equal to it’s probability. Let’s say Days with 50 demand has probability 0.1, so when we forecast 100 days, it will has the same output. 10 days with 50 demand, with random position.

Exponential & Normal Distribution

set.seed(2)

rate_param <- 1 / 70
demand_df <- rexp(10, rate = rate_param)

cat("10 Demands (Exponential):\n")
## 10 Demands (Exponential):
print(round(demand_df, 2))
##  [1] 130.57  28.33  10.27 121.15   6.27  46.68  75.21 105.81  92.00  10.96
mean_freq <- 20; sd_freq <- 5
data_frequency <- rnorm(10, mean = mean_freq, sd = sd_freq)
data_frequency <- abs(data_frequency)

cat("\n10 Data Frequency (Normal):\n")
## 
## 10 Data Frequency (Normal):
print(round(data_frequency, 2))
##  [1] 19.31 22.09 24.91 18.04 14.80 28.91  8.44 24.39 20.18 25.06

The exponential data shows that the time between demands is highly variable, with many short intervals and a few very long ones. The normal data indicates that the frequency is relatively stable, with most values clustered around the average of 20.

Monte Carlo simulation based on fitted distributions

It’s basically a Monte Carlo simulation that uses exponential demand and normal-based probabilities from recent simulation.

prob_sim <- data_frequency / sum(data_frequency)

simulation_from_exp <- function(n_day, seed = 42) {
  set.seed(seed)
  sample(demand_df, size = n_day, replace = TRUE, prob = prob_sim)
}

for (n in c(5, 20, 100, 1000)) {
  output <- simulation_from_exp(n)
  cat(sprintf("\n--- Predict %d Days ---\n", n))
  if (n <= 20) cat("Detail:", round(output, 2), "\n")
  cat(sprintf("Average: %.4f | Min: %.4f | Max: %.4f\n",
              mean(output), min(output), max(output)))
}
## 
## --- Predict 5 Days ---
## Detail: 6.27 6.27 10.27 121.15 92 
## Average: 47.1897 | Min: 6.2668 | Max: 121.1497
## 
## --- Predict 20 Days ---
## Detail: 6.27 6.27 10.27 121.15 92 28.33 130.57 46.68 92 92 105.81 130.57 6.27 10.96 105.81 6.27 75.21 46.68 105.81 28.33 
## Average: 62.3633 | Min: 6.2668 | Max: 130.5747
## 
## --- Predict 100 Days ---
## Average: 58.7043 | Min: 6.2668 | Max: 130.5747
## 
## --- Predict 1000 Days ---
## Average: 58.8888 | Min: 6.2668 | Max: 130.5747
cat("\nExponential Expectation Value (1/rate):", 1/rate_param, "\n")
## 
## Exponential Expectation Value (1/rate): 70

As the number of simulations increases, the average demand stabilizes around 59–62. This is slightly lower than the expected value of 70, suggesting the simulation does not fully match the theoretical distribution.