Simulasi Monte Carlo
Desmiwa Azhar
2026-06-01
Pendahuluan
Dalam analisis probabilitas dan statistika, nilai ekspektasi (expected value) merepresentasikan rata-rata jangka panjang suatu variabel acak. Secara teoritis, nilai ini dapat dihitung langsung dari distribusi probabilitas yang diketahui. Simulasi Monte Carlo menawarkan pendekatan empiris sebagai alternatif: dengan membangkitkan sampel acak berulang kali, kita dapat mendekati nilai parameter tersebut tanpa harus menyelesaikannya secara analitik.
Tugas ini merupakan modifikasi dari kode program yang diberikan pada modul perkuliahan. Perubahan yang dilakukan adalah memperluas jumlah simulasi menjadi 1.000, 5.000, dan 20.000 iterasi, kemudian membandingkan hasilnya dengan nilai ekspektasi teoritis dan menganalisis pola konvergensinya.
Data dan Persiapan
Data Awal
Sebuah warung es teh mencatat frekuensi harian permintaan selama 100 hari:
tabel_permintaan <- data.frame(
permintaan = c(50, 60, 70, 80, 90),
frekuensi = c(10, 20, 40, 20, 10)
)
tabel_permintaan## permintaan frekuensi
## 1 50 10
## 2 60 20
## 3 70 40
## 4 80 20
## 5 90 10Langkah 1 – Distribusi Probabilitas
Probabilitas setiap nilai permintaan dihitung sebagai frekuensi relatif terhadap total observasi (100 hari).
total_frekuensi <- sum(tabel_permintaan$frekuensi)
tabel_permintaan$probabilitas <- tabel_permintaan$frekuensi / total_frekuensi
tabel_permintaan## permintaan frekuensi probabilitas
## 1 50 10 0.1
## 2 60 20 0.2
## 3 70 40 0.4
## 4 80 20 0.2
## 5 90 10 0.1Langkah 2 – Distribusi Probabilitas Kumulatif
tabel_permintaan$probabilitas_kumulatif <- cumsum(tabel_permintaan$probabilitas)
tabel_permintaan## permintaan frekuensi probabilitas probabilitas_kumulatif
## 1 50 10 0.1 0.1
## 2 60 20 0.2 0.3
## 3 70 40 0.4 0.7
## 4 80 20 0.2 0.9
## 5 90 10 0.1 1.0Langkah 3 – Interval Bilangan Acak
Rentang bilangan acak 1–100 dibagi sesuai lebar probabilitas masing-masing nilai.
tabel_permintaan$batas_bawah <- c(
1,
head(tabel_permintaan$probabilitas_kumulatif, -1) * 100 + 1
)
tabel_permintaan$batas_atas <- tabel_permintaan$probabilitas_kumulatif * 100
tabel_permintaan## permintaan frekuensi probabilitas probabilitas_kumulatif batas_bawah
## 1 50 10 0.1 0.1 1
## 2 60 20 0.2 0.3 11
## 3 70 40 0.4 0.7 31
## 4 80 20 0.2 0.9 71
## 5 90 10 0.1 1.0 91
## batas_atas
## 1 10
## 2 30
## 3 70
## 4 90
## 5 100Nilai Ekspektasi Teoritis (Perhitungan Manual)
Nilai ekspektasi dihitung sebagai jumlahan dari setiap nilai permintaan dikalikan probabilitasnya:
\[E[X] = \sum_{i} x_i \cdot P(x_i) = 50(0{,}10) + 60(0{,}20) + 70(0{,}40) + 80(0{,}20) + 90(0{,}10)\]
nilai_ekspektasi <- sum(tabel_permintaan$permintaan * tabel_permintaan$probabilitas)
cat("Ekspektasi teoritis E[X] =", nilai_ekspektasi, "gelas/hari\n")## Ekspektasi teoritis E[X] = 70 gelas/hari\[E[X] = 5 + 12 + 28 + 16 + 9 = \mathbf{70 \text{ gelas/hari}}\]
Fungsi Simulasi Monte Carlo
Fungsi berikut diambil langsung dari kode program modul perkuliahan, tanpa perubahan logika.
permintaan_simulasi <- function(n, tabel) {
# Langkah 4: Bangkitkan bilangan acak diskrit 1–100
bilangan_acak <- sample(1:100, n, replace = TRUE)
# Petakan setiap bilangan acak ke nilai permintaan
get_demand <- function(x) {
index <- which(x >= tabel$batas_bawah & x <= tabel$batas_atas)
if (length(index) == 0) return(NA)
return(tabel$permintaan[index])
}
prediksi_permintaan <- sapply(bilangan_acak, get_demand)
result <- data.frame(
bilangan_acak = bilangan_acak,
prediksi_permintaan = prediksi_permintaan
)
na.omit(result)
}Pelaksanaan Simulasi (Tugas Modifikasi)
set.seed(42) digunakan agar hasil simulasi dapat
direproduksi.
set.seed(42)
sim_1000 <- permintaan_simulasi(1000, tabel_permintaan)
sim_5000 <- permintaan_simulasi(5000, tabel_permintaan)
sim_20000 <- permintaan_simulasi(20000, tabel_permintaan)Pratinjau Hasil — Simulasi 1.000 Hari
cat("Jumlah observasi:", nrow(sim_1000), "\n")## Jumlah observasi: 990cat("Rata-rata permintaan:", round(mean(sim_1000$prediksi_permintaan), 4), "gelas/hari\n\n")## Rata-rata permintaan: 70.101 gelas/harihead(sim_1000, 10)## bilangan_acak prediksi_permintaan
## 1 49 70
## 2 65 70
## 3 25 60
## 4 74 80
## 5 100 90
## 6 18 60
## 7 49 70
## 8 47 70
## 9 24 60
## 10 71 80Pratinjau Hasil — Simulasi 5.000 Hari
cat("Jumlah observasi:", nrow(sim_5000), "\n")## Jumlah observasi: 4959cat("Rata-rata permintaan:", round(mean(sim_5000$prediksi_permintaan), 4), "gelas/hari\n\n")## Rata-rata permintaan: 69.8891 gelas/harihead(sim_5000, 10)## bilangan_acak prediksi_permintaan
## 1 22 60
## 2 46 70
## 3 7 50
## 4 47 70
## 5 55 70
## 6 38 70
## 7 43 70
## 8 74 80
## 9 94 90
## 10 49 70Pratinjau Hasil — Simulasi 20.000 Hari
cat("Jumlah observasi:", nrow(sim_20000), "\n")## Jumlah observasi: 19805cat("Rata-rata permintaan:", round(mean(sim_20000$prediksi_permintaan), 4), "gelas/hari\n\n")## Rata-rata permintaan: 70.0722 gelas/harihead(sim_20000, 10)## bilangan_acak prediksi_permintaan
## 1 100 90
## 2 56 70
## 3 2 50
## 4 84 80
## 5 32 70
## 6 96 90
## 7 35 70
## 8 94 90
## 9 61 70
## 11 2 50Perbandingan Hasil Simulasi dengan Ekspektasi Teoritis
perbandingan <- data.frame(
Metode = c(
"Ekspektasi Teoritis (Manual)",
"Simulasi 1.000 hari",
"Simulasi 5.000 hari",
"Simulasi 20.000 hari"
),
Rata_Rata = c(
nilai_ekspektasi,
mean(sim_1000$prediksi_permintaan),
mean(sim_5000$prediksi_permintaan),
mean(sim_20000$prediksi_permintaan)
)
)
perbandingan$Rata_Rata <- round(perbandingan$Rata_Rata, 4)
perbandingan$Selisih_Absolut <- round(abs(perbandingan$Rata_Rata - nilai_ekspektasi), 4)
perbandingan$Selisih_Persen <- round(
abs(perbandingan$Rata_Rata - nilai_ekspektasi) / nilai_ekspektasi * 100, 4
)
perbandingan## Metode Rata_Rata Selisih_Absolut Selisih_Persen
## 1 Ekspektasi Teoritis (Manual) 70.0000 0.0000 0.0000
## 2 Simulasi 1.000 hari 70.1010 0.1010 0.1443
## 3 Simulasi 5.000 hari 69.8891 0.1109 0.1584
## 4 Simulasi 20.000 hari 70.0722 0.0722 0.1031Visualisasi Konvergensi
Grafik Konvergensi Running Mean — Simulasi 20.000 Hari
Grafik ini menunjukkan bagaimana rata-rata simulasi bergerak mendekati ekspektasi teoritis seiring bertambahnya jumlah iterasi.
library(ggplot2)
n_aktual <- nrow(sim_20000)
running_avg <- cumsum(sim_20000$prediksi_permintaan) / seq_len(n_aktual)
df_konvergensi <- data.frame(
iterasi = seq_len(n_aktual),
running_avg = running_avg
)
ggplot(df_konvergensi, aes(x = iterasi, y = running_avg)) +
geom_line(color = "#2E86C1", linewidth = 0.5, alpha = 0.8) +
geom_hline(
yintercept = nilai_ekspektasi,
color = "#E74C3C",
linetype = "dashed",
linewidth = 0.8
) +
annotate(
"text",
x = n_aktual * 0.75,
y = nilai_ekspektasi + 0.6,
label = paste("Ekspektasi Teoritis =", nilai_ekspektasi),
color = "#E74C3C",
size = 3.5
) +
scale_x_continuous(labels = scales::comma) +
labs(
title = "Konvergensi Running Mean – Simulasi Monte Carlo (20.000 Iterasi)",
subtitle = "Garis merah putus-putus = nilai ekspektasi teoritis (E[X] = 70 gelas/hari)",
x = "Jumlah Iterasi",
y = "Rata-rata Kumulatif (gelas/hari)"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
Distribusi Permintaan Hasil Simulasi
# Gabungkan ketiga simulasi untuk perbandingan distribusi
df_dist <- data.frame(
permintaan = c(
sim_1000$prediksi_permintaan,
sim_5000$prediksi_permintaan,
sim_20000$prediksi_permintaan
),
simulasi = factor(
rep(
c("1.000 hari", "5.000 hari", "20.000 hari"),
times = c(nrow(sim_1000), nrow(sim_5000), nrow(sim_20000))
),
levels = c("1.000 hari", "5.000 hari", "20.000 hari")
)
)
# Hitung proporsi empiris
df_prop <- aggregate(
x = list(frekuensi = df_dist$permintaan),
by = list(simulasi = df_dist$simulasi, permintaan = df_dist$permintaan),
FUN = length
)
df_prop <- do.call(rbind, lapply(split(df_prop, df_prop$simulasi), function(x) {
x$proporsi <- x$frekuensi / sum(x$frekuensi)
x
}))
# Probabilitas teoritis sebagai referensi
df_teoritis <- data.frame(
permintaan = c(50, 60, 70, 80, 90),
proporsi = c(0.10, 0.20, 0.40, 0.20, 0.10),
simulasi = "Teoritis"
)
ggplot(df_prop, aes(x = factor(permintaan), y = proporsi, fill = simulasi)) +
geom_col(position = "dodge", alpha = 0.85) +
geom_point(
data = df_teoritis,
aes(x = factor(permintaan), y = proporsi),
color = "#E74C3C", size = 3, shape = 18,
inherit.aes = FALSE
) +
geom_line(
data = df_teoritis,
aes(x = factor(permintaan), y = proporsi, group = 1),
color = "#E74C3C", linewidth = 0.8, linetype = "dashed",
inherit.aes = FALSE
) +
scale_fill_manual(values = c("#5DADE2", "#1ABC9C", "#F39C12")) +
labs(
title = "Distribusi Permintaan Hasil Simulasi vs. Teoritis",
subtitle = "Garis merah = distribusi teoritis; kolom = hasil simulasi",
x = "Permintaan (gelas)",
y = "Proporsi",
fill = "Jumlah Simulasi"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
Analisis dan Pembahasan
Pengaruh Jumlah Simulasi terhadap Kedekatan Hasil
Berdasarkan tabel perbandingan yang telah dihitung, dapat diamati pola berikut:
## Jumlah Simulasi Rata-rata Simulasi Selisih Absolut
## 1 1.000 70.1010 0.1010
## 2 5.000 69.8891 0.1109
## 3 20.000 70.0722 0.0722Hasil simulasi menunjukkan bahwa semakin besar jumlah iterasi, rata-rata simulasi semakin mendekati nilai ekspektasi teoritis sebesar 70 gelas/hari. Grafik konvergensi memperlihatkan bahwa pada iterasi-iterasi awal, running mean masih berfluktuasi cukup jauh dari nilai teoritis, namun kemudian berangsur stabil dan bergerak semakin dekat ke angka 70 seiring bertambahnya jumlah simulasi.
Mengapa Peningkatan Jumlah Simulasi Menghasilkan Estimasi yang Lebih Stabil?
Fenomena ini dijelaskan oleh Hukum Bilangan Besar (Law of Large Numbers, LLN):
Jika \(X_1, X_2, \ldots, X_n\) merupakan variabel acak independen dan identik terdistribusi (i.i.d.) dengan ekspektasi \(E[X] = \mu\), maka rata-rata sampel \(\bar{X}_n\) akan konvergen ke \(\mu\) ketika \(n \to \infty\): \[\bar{X}_n = \frac{1}{n}\sum_{i=1}^{n} X_i \xrightarrow{p} \mu\]
Dalam konteks simulasi Monte Carlo ini, setiap iterasi adalah satu pengambilan sampel acak dari distribusi permintaan. Ada dua penjelasan statistik yang mendasari mengapa simulasi lebih besar memberikan hasil lebih akurat:
1. Varians Rata-rata Sampel Mengecil Seiring Bertambahnya \(n\)
Varians dari rata-rata sampel adalah \(\text{Var}(\bar{X}_n) = \dfrac{\sigma^2}{n}\), di mana \(\sigma^2\) adalah varians populasi. Artinya, ketika \(n\) diperbesar empat kali lipat (misalnya dari 1.000 ke 4.000), standar deviasi rata-rata sampel berkurang setengahnya. Inilah sebabnya simulasi 20.000 iterasi menghasilkan estimasi yang lebih konsisten dibandingkan simulasi 1.000 iterasi.
2. Laju Konvergensi Monte Carlo
Laju konvergensi metode Monte Carlo adalah \(O\!\left(\dfrac{1}{\sqrt{n}}\right)\). Ini berarti kesalahan estimasi berbanding terbalik dengan akar kuadrat jumlah iterasi. Secara praktis:
| Jumlah Simulasi (\(n\)) | Faktor Pengurangan Galat (\(\sqrt{n}\)) |
|---|---|
| 1.000 | ≈ 31,6 |
| 5.000 | ≈ 70,7 |
| 20.000 | ≈ 141,4 |
Artinya, simulasi 20.000 iterasi secara teoritis memiliki galat estimasi sekitar 4,5 kali lebih kecil dibandingkan simulasi 1.000 iterasi.
3. Cakupan Sampel yang Lebih Representatif
Dengan jumlah iterasi yang lebih besar, distribusi empiris hasil simulasi semakin mendekati distribusi teoritis. Hal ini terlihat jelas pada grafik distribusi di atas, di mana proporsi tiap nilai permintaan pada simulasi 20.000 hari jauh lebih dekat ke nilai teoritis (50: 10%, 60: 20%, 70: 40%, 80: 20%, 90: 10%) dibandingkan simulasi 1.000 hari.
Implikasi Praktis
Bagi pengambil keputusan bisnis, hasil ini memberikan wawasan penting: menggunakan lebih banyak data historis atau menjalankan lebih banyak skenario simulasi akan menghasilkan prediksi yang lebih andal. Dalam kasus warung es teh ini, estimasi permintaan rata-rata sebesar 70 gelas/hari dapat dijadikan acuan yang cukup akurat untuk perencanaan stok, dengan catatan bahwa estimasi tersebut hanya valid selama pola permintaan historis tetap berlaku.
Kesimpulan
Nilai ekspektasi teoritis dari distribusi permintaan es teh adalah E[X] = 70 gelas/hari, dihitung sebagai \(\sum x_i \cdot P(x_i)\).
Simulasi Monte Carlo dengan 1.000, 5.000, dan 20.000 iterasi menghasilkan rata-rata masing-masing yang mendekati nilai teoritis tersebut.
Peningkatan jumlah simulasi menurunkan selisih absolut antara hasil simulasi dan ekspektasi teoritis, sesuai dengan Hukum Bilangan Besar.
Laju konvergensi Monte Carlo sebesar \(O(1/\sqrt{n})\) menjelaskan mengapa penggandaan iterasi tidak serta-merta menggandakan akurasi — dibutuhkan peningkatan \(n\) sebesar 4× untuk menurunkan galat sebesar 2×.
Simulasi 20.000 iterasi memberikan estimasi yang paling stabil dan paling dekat ke nilai teoritis, sehingga dapat diandalkan untuk pengambilan keputusan.