Dalam berbagai bidang ilmu, seperti statistika, teknik, ekonomi, manajemen operasional, hingga ilmu komputer, proses pengambilan keputusan sering kali dilakukan di bawah kondisi yang mengandung ketidakpastian. Ketidakpastian tersebut muncul karena nilai suatu variabel tidak selalu dapat diprediksi secara pasti, melainkan mengikuti suatu distribusi probabilitas tertentu. Oleh karena itu, diperlukan metode yang mampu memberikan pendekatan numerik terhadap karakteristik suatu distribusi ketika penyelesaian analitik menjadi sulit atau bahkan tidak memungkinkan.
Salah satu karakteristik penting dalam distribusi probabilitas adalah nilai ekspektasi (expected value). Nilai ekspektasi menggambarkan rata-rata jangka panjang yang diharapkan diperoleh apabila suatu percobaan dilakukan berulang kali dalam jumlah yang sangat besar. Pada distribusi yang sederhana, nilai ekspektasi dapat dihitung secara langsung menggunakan rumus matematis berdasarkan peluang setiap kejadian. Namun, pada permasalahan nyata yang lebih kompleks, distribusi probabilitas sering kali tidak diketahui secara eksplisit atau melibatkan perhitungan yang rumit sehingga pendekatan simulasi menjadi alternatif yang efektif.
Metode yang paling banyak digunakan untuk mengatasi permasalahan tersebut adalah Simulasi Monte Carlo. Metode ini merupakan teknik numerik berbasis pembangkitan bilangan acak yang digunakan untuk memperkirakan nilai suatu parameter statistik melalui proses pengambilan sampel secara berulang. Dengan memanfaatkan prinsip probabilitas, hasil simulasi akan semakin mendekati nilai sebenarnya ketika jumlah pengulangan simulasi terus ditingkatkan.
Secara teoritis, fenomena tersebut dijelaskan oleh Hukum Bilangan Besar (Law of Large Numbers) (Ross, 2020) yang menyatakan bahwa rata-rata sampel akan semakin mendekati nilai ekspektasi populasi seiring bertambahnya ukuran sampel. Selain itu, Teorema Limit Pusat (Central Limit Theorem) (Devore, 2016) menjelaskan bahwa distribusi rata-rata sampel akan cenderung mengikuti distribusi normal ketika jumlah observasi cukup besar, sehingga estimasi yang dihasilkan menjadi lebih stabil dan memiliki tingkat kesalahan yang semakin kecil.
Pada praktikum ini digunakan data distribusi permintaan diskrit yang terdiri atas lima tingkat permintaan, yaitu 50, 60, 70, 80, dan 90 unit dengan probabilitas yang telah diketahui. Berdasarkan distribusi tersebut, nilai ekspektasi teoritis dihitung secara manual sebagai acuan dalam mengevaluasi hasil simulasi. Selanjutnya dilakukan modifikasi program Monte Carlo dengan menggunakan tiga ukuran simulasi yang berbeda, yaitu 1.000, 5.000, dan 20.000 iterasi.
Perbedaan jumlah simulasi tersebut memungkinkan dilakukan evaluasi terhadap bagaimana peningkatan ukuran simulasi memengaruhi ketepatan estimasi nilai ekspektasi, tingkat kesalahan estimasi, kestabilan hasil simulasi, serta proses konvergensi menuju nilai teoritis. Dengan demikian, praktikum ini tidak hanya menunjukkan implementasi algoritma Monte Carlo, tetapi juga memberikan pemahaman empiris mengenai hubungan antara ukuran simulasi dengan kualitas hasil estimasi yang diperoleh.
Selain analisis deskriptif dan visual, dilakukan pula pengujian statistik inferensial berupa One Sample t-test, ANOVA, dan uji homogenitas varians untuk menguji secara formal apakah perbedaan rata-rata dan varians antar simulasi bersifat signifikan secara statistik.
Praktikum ini bertujuan untuk:
Tahap persiapan merupakan fondasi utama sebelum proses simulasi Monte Carlo dijalankan. Pada tahap ini dilakukan pemanggilan berbagai paket (library) yang mendukung proses komputasi, manipulasi data, visualisasi hasil, serta penyajian tabel. Selain itu, ditentukan pula parameter awal yang akan digunakan secara konsisten selama proses simulasi sehingga seluruh hasil dapat direproduksi (reproducible) dan dianalisis secara objektif.
Penggunaan paket yang tepat tidak hanya mempermudah proses pemrograman, tetapi juga meningkatkan kualitas visualisasi serta efisiensi pengolahan data. Oleh karena itu, setiap paket yang digunakan memiliki fungsi spesifik sesuai kebutuhan analisis.
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(scales)
library(knitr)
library(boot)
Penjelasan Paket:
Pada praktikum ini digunakan distribusi permintaan diskrit yang telah disediakan pada modul pembelajaran. Distribusi tersebut terdiri atas lima tingkat permintaan dengan probabilitas yang telah diketahui berdasarkan frekuensi historis.
Sebagai bentuk modifikasi terhadap program awal, jumlah simulasi diubah menjadi tiga skenario, yaitu:
Skenario Jumlah Simulasi (N) Simulasi I 1.000 Simulasi II 5.000 Simulasi III 20.000
Nilai ekspektasi teoritis yang diperoleh dari distribusi probabilitas akan dijadikan benchmark untuk mengevaluasi ketepatan hasil simulasi pada setiap skenario.
Agar hasil simulasi dapat direproduksi pada waktu yang berbeda maupun pada perangkat yang berbeda, digunakan seed yang telah ditentukan sebelum proses pembangkitan bilangan acak.
parameter_simulasi <- data.frame(
Simulasi = c("N = 1000",
"N = 5000",
"N = 20000"),
N = c(1000,
5000,
20000),
Seed = c(101,
202,
303)
)
parameter_simulasi
## Simulasi N Seed
## 1 N = 1000 1000 101
## 2 N = 5000 5000 202
## 3 N = 20000 20000 303
Pemilihan seed yang berbeda pada setiap skenario bertujuan untuk menghasilkan urutan bilangan acak yang independen, namun tetap memberikan hasil simulasi yang dapat direproduksi. Dengan demikian, setiap proses analisis dapat diulang dan diverifikasi tanpa mengubah karakteristik distribusi probabilitas yang digunakan.
Secara umum, tahapan simulasi Monte Carlo pada praktikum ini dilakukan melalui beberapa langkah berikut.
Tahap pembangunan data merupakan langkah awal dalam implementasi Simulasi Monte Carlo. Pada tahap ini, data historis permintaan diubah menjadi distribusi probabilitas yang selanjutnya digunakan sebagai dasar pembangkitan bilangan acak selama proses simulasi.
Distribusi probabilitas memiliki peran penting karena menentukan peluang kemunculan setiap tingkat permintaan. Dengan kata lain, kualitas hasil simulasi sangat bergantung pada representasi distribusi probabilitas yang dibangun dari data awal. Oleh sebab itu, sebelum simulasi dilakukan, perlu dihitung probabilitas masing-masing nilai permintaan beserta probabilitas kumulatifnya untuk membentuk interval bilangan acak.
Selain itu, pada tahap ini juga dilakukan perhitungan nilai ekspektasi teoritis (expected value) secara manual. Nilai tersebut berfungsi sebagai nilai acuan (benchmark) dalam mengevaluasi tingkat akurasi estimasi yang diperoleh melalui Simulasi Monte Carlo.
Data permintaan yang digunakan merupakan distribusi diskrit dengan lima kemungkinan tingkat permintaan. Setiap tingkat permintaan memiliki frekuensi kemunculan yang berbeda sehingga menghasilkan probabilitas yang berbeda pula.
tabel_permintaan <- data.frame(
permintaan = c(50, 60, 70, 80, 90),
frekuensi = c(10, 20, 40, 20, 10)
)
# Total frekuensi
total_frekuensi <- sum(tabel_permintaan$frekuensi)
# Probabilitas
tabel_permintaan$probabilitas <-
tabel_permintaan$frekuensi / total_frekuensi
# Probabilitas kumulatif
tabel_permintaan$prob_kumulatif <-
cumsum(tabel_permintaan$probabilitas)
# Interval bilangan acak
tabel_permintaan$batas_bawah <- c(
1,
head(tabel_permintaan$prob_kumulatif, -1) * 100 + 1
)
tabel_permintaan$batas_atas <-
tabel_permintaan$prob_kumulatif * 100
knitr::kable(
tabel_permintaan,
digits = 2,
caption = "Distribusi Permintaan"
)
| permintaan | frekuensi | probabilitas | prob_kumulatif | batas_bawah | batas_atas |
|---|---|---|---|---|---|
| 50 | 10 | 0.1 | 0.1 | 1 | 10 |
| 60 | 20 | 0.2 | 0.3 | 11 | 30 |
| 70 | 40 | 0.4 | 0.7 | 31 | 70 |
| 80 | 20 | 0.2 | 0.9 | 71 | 90 |
| 90 | 10 | 0.1 | 1.0 | 91 | 100 |
Tabel Distribusi Permintaan
| Permintaan | Frekuensi | Probabilitas | Probabilitas Kumulatif | Interval Bilangan Acak |
|---|---|---|---|---|
| 50 | 10 | 0,10 | 0,10 | 1–10 |
| 60 | 20 | 0,20 | 0,30 | 11–30 |
| 70 | 40 | 0,40 | 0,70 | 31–70 |
| 80 | 20 | 0,20 | 0,90 | 71–90 |
| 90 | 10 | 0,10 | 1,00 | 91–100 |
Berdasarkan Tabel 3.1 dapat diketahui bahwa nilai permintaan sebesar 70 unit memiliki probabilitas terbesar, yaitu 40%, sehingga nilai tersebut mempunyai peluang paling tinggi untuk muncul selama proses simulasi.
Sementara itu, permintaan sebesar 50 unit dan 90 unit masing-masing hanya memiliki probabilitas 10%, sehingga frekuensi kemunculannya diperkirakan lebih sedikit dibandingkan nilai permintaan lainnya.
Interval bilangan acak dibangun berdasarkan probabilitas kumulatif sehingga setiap angka acak dari 1 hingga 100 dapat dipetakan secara langsung ke salah satu tingkat permintaan. Sebagai contoh, apabila bilangan acak yang dihasilkan adalah 25, maka angka tersebut berada pada interval 11–30, sehingga nilai permintaan yang dipilih adalah 60 unit.
Pendekatan ini merupakan prinsip dasar metode Monte Carlo, yaitu menggantikan distribusi probabilitas dengan representasi bilangan acak yang memiliki peluang kemunculan sesuai distribusi aslinya.
Nilai ekspektasi (Expected Value) merupakan rata-rata teoritis suatu variabel acak yang diperoleh berdasarkan distribusi probabilitasnya. Pada distribusi diskrit, nilai ekspektasi dihitung menggunakan persamaan berikut.
\[ E(X)=\sum_{i=1}^{n} x_iP(x_i) \]
dengan:
Berdasarkan distribusi pada Tabel 3.1, diperoleh perhitungan sebagai berikut.
\[ E(X)=50(0.10)+60(0.20)+70(0.40)+80(0.20)+90(0.10)=70 \]
ekspektasi_teoritis <-
sum(
tabel_permintaan$permintaan *
tabel_permintaan$probabilitas
)
ekspektasi_teoritis
## [1] 70
Perhitungan manual maupun menggunakan R menghasilkan nilai ekspektasi sebesar 70 unit. Nilai ini menunjukkan bahwa apabila proses permintaan diamati dalam jangka panjang dengan distribusi probabilitas yang tetap, maka rata-rata permintaan akan mendekati 70 unit.
Nilai ekspektasi teoritis tersebut akan dijadikan tolok ukur utama dalam mengevaluasi hasil Simulasi Monte Carlo. Semakin kecil selisih antara rata-rata hasil simulasi dengan nilai teoritis, semakin tinggi tingkat akurasi estimasi yang dihasilkan oleh simulasi.
Untuk memperoleh gambaran yang lebih intuitif mengenai karakteristik distribusi probabilitas, dilakukan visualisasi menggunakan diagram batang. Visualisasi ini membantu menunjukkan tingkat dominasi masing-masing nilai permintaan sebelum proses simulasi dijalankan.
ggplot(
tabel_permintaan,
aes(
x = factor(permintaan),
y = probabilitas
)
) +
geom_col(width = 0.65) +
geom_text(
aes(
label = scales::percent(probabilitas)
),
vjust = -0.4,
size = 4
) +
labs(
title = "Distribusi Probabilitas Permintaan",
x = "Permintaan",
y = "Probabilitas"
) +
theme_minimal(base_size = 13)
Visualisasi distribusi menunjukkan bahwa pola probabilitas membentuk distribusi yang simetris dengan puncak pada permintaan 70 unit. Nilai tersebut memiliki probabilitas tertinggi sehingga diperkirakan akan paling sering muncul selama proses simulasi Monte Carlo.
Sebaliknya, nilai permintaan 50 dan 90 unit memiliki probabilitas terendah, sehingga frekuensi kemunculannya relatif lebih sedikit. Pola distribusi yang simetris ini juga mengindikasikan bahwa nilai ekspektasi berada tepat di pusat distribusi, yaitu sebesar 70 unit, sehingga secara teoritis rata-rata hasil simulasi akan semakin mendekati nilai tersebut ketika jumlah iterasi diperbesar.
Setelah distribusi probabilitas permintaan berhasil dibangun, langkah berikutnya adalah mengimplementasikan algoritma Simulasi Monte Carlo. Metode ini bekerja dengan menghasilkan sejumlah bilangan acak yang kemudian dipetakan ke dalam interval probabilitas yang telah ditentukan sebelumnya. Pendekatan ini sesuai dengan prinsip dasar Simulasi Monte Carlo yang menggunakan pembangkitan bilangan acak untuk merepresentasikan distribusi probabilitas suatu variabel acak (Fishman, 1996). Setiap bilangan acak akan merepresentasikan satu kejadian permintaan sesuai distribusi probabilitas awal.
Pada praktikum ini, fungsi simulasi dirancang agar lebih fleksibel dibandingkan program dasar pada modul. Pengguna hanya perlu menentukan jumlah simulasi (n), kemudian fungsi secara otomatis menghasilkan data permintaan hasil simulasi beserta berbagai ukuran statistik yang diperlukan untuk proses analisis.
Selain menghasilkan nilai rata-rata hasil simulasi, fungsi ini juga menghitung simpangan baku (standard deviation), standard error, interval kepercayaan, serta galat estimasi terhadap nilai ekspektasi teoritis. Informasi tersebut akan digunakan pada bab selanjutnya untuk mengevaluasi kualitas estimasi Monte Carlo pada berbagai ukuran simulasi.
Secara umum, proses Simulasi Monte Carlo pada penelitian ini dilakukan melalui tahapan sebagai berikut.
simulasi_mc <- function(n, tabel, seed){
set.seed(seed)
# Membuat bilangan acak
bilangan_acak <- sample(
1:100,
size = n,
replace = TRUE
)
# Menentukan indeks berdasarkan interval probabilitas
# round() digunakan untuk menghindari masalah floating-point
# pada batas_bawah dan batas_atas (mis. 11.000000000000002)
batas_bawah_bulat <- round(tabel$batas_bawah)
batas_atas_bulat <- round(tabel$batas_atas)
get_idx <- function(x) {
which(x >= batas_bawah_bulat & x <= batas_atas_bulat)
}
idx <- sapply(bilangan_acak, get_idx)
# Mengubah menjadi permintaan
permintaan <- tabel$permintaan[idx]
# Running mean
running_mean <- cumsum(permintaan) / seq_len(n)
# Mengembalikan seluruh hasil
return(
list(
data = data.frame(
Iterasi = seq_len(n),
Bilangan_Acak = bilangan_acak,
Permintaan = permintaan,
Running_Mean = running_mean
),
mean = mean(permintaan),
sd = sd(permintaan),
se = sd(permintaan) / sqrt(n),
ci_lower = mean(permintaan) -
1.96 * sd(permintaan) / sqrt(n),
ci_upper = mean(permintaan) +
1.96 * sd(permintaan) / sqrt(n)
)
)
}
Penjelasan Fungsi
Fungsi simulasi_mc() dikembangkan agar seluruh proses simulasi dapat dilakukan hanya melalui satu pemanggilan fungsi. Pengguna cukup menentukan jumlah simulasi dan seed, sedangkan seluruh proses pembangkitan bilangan acak, pemetaan interval probabilitas, hingga perhitungan statistik dilakukan secara otomatis.
Berbeda dengan kode dasar pada modul yang hanya menghasilkan rata-rata permintaan, fungsi ini juga menghasilkan data simulasi lengkap yang dapat digunakan untuk berbagai analisis lanjutan, seperti visualisasi distribusi hasil simulasi, analisis konvergensi (running mean), perhitungan standard error, serta pembentukan interval kepercayaan.
Pendekatan ini membuat kode menjadi lebih modular, mudah dipelihara, dan dapat digunakan kembali (reusable) pada berbagai skenario simulasi tanpa perlu melakukan perubahan terhadap struktur program utama.
Sebelum digunakan pada seluruh skenario simulasi, fungsi terlebih dahulu diuji menggunakan jumlah simulasi yang relatif kecil untuk memastikan bahwa proses pemetaan bilangan acak ke distribusi permintaan telah berjalan dengan benar.
uji <- simulasi_mc(
n = 10,
tabel = tabel_permintaan,
seed = 123
)
head(uji$data)
## Iterasi Bilangan_Acak Permintaan Running_Mean
## 1 1 31 70 70.00000
## 2 2 79 80 75.00000
## 3 3 51 70 73.33333
## 4 4 14 60 70.00000
## 5 5 67 70 70.00000
## 6 6 42 70 70.00000
Hasil pengujian menunjukkan bahwa setiap bilangan acak berhasil dipetakan ke nilai permintaan sesuai interval probabilitas yang telah dibangun pada Bab 3. Nilai Running Mean juga berubah secara dinamis mengikuti setiap iterasi, yang nantinya akan digunakan untuk mengamati proses konvergensi estimasi menuju nilai ekspektasi teoritis.
Setiap pemanggilan fungsi simulasi_mc() menghasilkan beberapa informasi penting yang akan digunakan pada analisis selanjutnya, yaitu:
Setelah distribusi probabilitas berhasil dibangun dan fungsi Simulasi Monte Carlo selesai dikembangkan, tahap berikutnya adalah menjalankan simulasi pada tiga ukuran sampel yang berbeda, yaitu 1.000, 5.000, dan 20.000 iterasi.
Variasi jumlah simulasi ini bertujuan untuk mengevaluasi pengaruh ukuran sampel terhadap kualitas estimasi nilai ekspektasi. Secara teoritis, semakin banyak iterasi yang dilakukan, maka rata-rata hasil simulasi akan semakin mendekati nilai ekspektasi teoritis. Fenomena tersebut merupakan implementasi langsung dari Law of Large Numbers (LLN) yang menyatakan bahwa rata-rata sampel akan mengalami konvergensi menuju rata-rata populasi ketika ukuran sampel terus bertambah.
Pada setiap skenario simulasi dihitung beberapa ukuran statistik yang digunakan sebagai indikator kualitas estimasi, meliputi:
Nilai-nilai tersebut kemudian dibandingkan dengan nilai ekspektasi teoritis sebesar 70 unit.
sim_1000 <- simulasi_mc(
parameter_simulasi$N[1],
tabel_permintaan,
parameter_simulasi$Seed[1]
)
sim_5000 <- simulasi_mc(
parameter_simulasi$N[2],
tabel_permintaan,
parameter_simulasi$Seed[2]
)
sim_20000 <- simulasi_mc(
parameter_simulasi$N[3],
tabel_permintaan,
parameter_simulasi$Seed[3]
)
Tahap ini menghasilkan ringkasan statistik dari setiap simulasi sehingga mempermudah proses evaluasi terhadap kualitas estimasi yang diperoleh.
Untuk mengetahui tingkat ketepatan hasil simulasi terhadap nilai teoritis, dilakukan perhitungan dua ukuran galat, yaitu:
1. Absolute Error
\[ AE = |\bar X - E(X)| \] Absolute Error menunjukkan besarnya selisih rata-rata simulasi terhadap nilai ekspektasi teoritis.
2. Relative Error
\[ RE=\frac{|\bar X-E(X)|}{E(X)}\times100\% \] Relative Error menunjukkan besarnya persentase kesalahan estimasi dibandingkan nilai teoritis.
ekspektasi <- ekspektasi_teoritis
error_summary <- data.frame(
Simulasi = c(
"1000",
"5000",
"20000"
),
Mean = c(
sim_1000$mean,
sim_5000$mean,
sim_20000$mean
)
)
error_summary$Absolute_Error <-
abs(
error_summary$Mean -
ekspektasi
)
error_summary$Relative_Error <-
100 *
error_summary$Absolute_Error /
ekspektasi
hasil_perbandingan <- data.frame(
Jumlah_Simulasi = c(1000, 5000, 20000),
Mean = c(
sim_1000$mean,
sim_5000$mean,
sim_20000$mean
),
SD = c(
sim_1000$sd,
sim_5000$sd,
sim_20000$sd
),
Standard_Error = c(
sim_1000$se,
sim_5000$se,
sim_20000$se
),
CI_Bawah = c(
sim_1000$ci_lower,
sim_5000$ci_lower,
sim_20000$ci_lower
),
CI_Atas = c(
sim_1000$ci_upper,
sim_5000$ci_upper,
sim_20000$ci_upper
),
Ekspektasi_Teoritis = rep(
ekspektasi_teoritis,
3
)
)
hasil_perbandingan$Absolute_Error <-
abs(
hasil_perbandingan$Mean -
hasil_perbandingan$Ekspektasi_Teoritis
)
hasil_perbandingan$Relative_Error <-
100 *
hasil_perbandingan$Absolute_Error /
hasil_perbandingan$Ekspektasi_Teoritis
Pada simulasi 1.000 iterasi, diperoleh rata-rata sebesar 70,0700 dengan Absolute Error sebesar 0,0700. Ketika jumlah simulasi ditingkatkan menjadi 5.000 iterasi, rata-rata simulasi menjadi 69,9280 dengan Absolute Error 0,0720, yang sedikit lebih besar dibandingkan simulasi 1.000 iterasi. Sementara itu, simulasi 20.000 iterasi menghasilkan rata-rata 70,0235 dengan Absolute Error 0,0235, yang merupakan galat terkecil di antara ketiga skenario. Perbedaan kecil pada simulasi 1.000 dan 5.000 iterasi merupakan konsekuensi alami dari proses pembangkitan bilangan acak sehingga nilai galat tidak selalu menurun secara monoton pada setiap percobaan.
Meskipun Absolute Error sempat sedikit meningkat dari 1.000 ke 5.000 iterasi, nilai Standard Error terus menurun secara konsisten dari 0,3391 menjadi 0,1550, kemudian 0,0778. Penurunan Standard Error menunjukkan bahwa estimasi rata-rata menjadi semakin presisi seiring bertambahnya jumlah simulasi. Dengan demikian, peningkatan jumlah iterasi tetap menghasilkan estimasi yang lebih stabil dan memiliki tingkat ketidakpastian yang lebih rendah, sekaligus pada akhirnya juga memberikan Absolute Error terkecil pada ukuran simulasi terbesar (20.000 iterasi).
Tahap selanjutnya adalah menyusun tabel perbandingan yang merangkum seluruh hasil simulasi Monte Carlo pada tiga ukuran simulasi, yaitu 1.000, 5.000, dan 20.000 iterasi. Penyajian dalam bentuk tabel bertujuan untuk mempermudah proses evaluasi terhadap perubahan kualitas estimasi yang dihasilkan ketika jumlah simulasi ditingkatkan.
Selain membandingkan rata-rata hasil simulasi dengan nilai ekspektasi teoritis, tabel ini juga menyajikan beberapa ukuran statistik penting seperti simpangan baku (Standard Deviation), galat baku (Standard Error), interval kepercayaan 95%, galat absolut (Absolute Error), dan galat relatif (Relative Error). Kombinasi indikator tersebut memberikan gambaran yang lebih komprehensif mengenai tingkat akurasi, kestabilan, serta ketelitian estimasi Monte Carlo.
knitr::kable(
hasil_perbandingan,
digits = 4,
caption = "Perbandingan Hasil Simulasi Monte Carlo"
)
| Jumlah_Simulasi | Mean | SD | Standard_Error | CI_Bawah | CI_Atas | Ekspektasi_Teoritis | Absolute_Error | Relative_Error |
|---|---|---|---|---|---|---|---|---|
| 1000 | 70.0700 | 10.7243 | 0.3391 | 69.4053 | 70.7347 | 70 | 0.0700 | 0.1000 |
| 5000 | 69.9280 | 10.9608 | 0.1550 | 69.6242 | 70.2318 | 70 | 0.0720 | 0.1029 |
| 20000 | 70.0235 | 10.9977 | 0.0778 | 69.8711 | 70.1759 | 70 | 0.0235 | 0.0336 |
Tabel perbandingan menunjukkan bahwa seluruh hasil simulasi memberikan estimasi rata-rata yang sangat dekat dengan nilai ekspektasi teoritis sebesar 70 unit. Hal ini mengindikasikan bahwa algoritma Monte Carlo berhasil merepresentasikan distribusi probabilitas permintaan dengan baik.
Nilai Standard Deviation pada ketiga simulasi relatif sama, yaitu sekitar 11 (10,7243 hingga 10,9977), karena seluruh simulasi berasal dari distribusi probabilitas yang sama. Sebaliknya, Standard Error mengalami penurunan yang konsisten ketika jumlah simulasi diperbesar, yaitu dari 0,3391 pada 1.000 iterasi menjadi 0,1550 pada 5.000 iterasi, dan 0,0778 pada 20.000 iterasi. Penurunan ini menunjukkan bahwa estimasi rata-rata menjadi semakin presisi.
Nilai Confidence Interval juga menjadi semakin sempit pada ukuran simulasi yang lebih besar. Hal tersebut menunjukkan bahwa tingkat ketidakpastian estimasi semakin kecil sehingga rata-rata hasil simulasi semakin dapat diandalkan.
Absolute Error sempat meningkat tipis dari 1.000 ke 5.000 iterasi (dari 0,0700 menjadi 0,0720), namun kemudian mengecil signifikan pada 20.000 iterasi (0,0235) — nilai Absolute Error terkecil di antara ketiga skenario. Pola ini menunjukkan bahwa meskipun variasi acak dapat menyebabkan fluktuasi kecil antarhasil simulasi pada ukuran sampel sedang, peningkatan jumlah simulasi yang signifikan (20.000 iterasi) tetap menghasilkan estimasi yang paling mendekati nilai ekspektasi teoritis, sejalan dengan prediksi Law of Large Numbers.
Penurunan Standard Error sejalan dengan teori Monte Carlo yang menyatakan bahwa besar galat estimasi berkurang sebanding dengan 1/√N.
Berdasarkan hasil yang diperoleh, beberapa kecenderungan umum yang diharapkan muncul adalah sebagai berikut.
| Jumlah Simulasi | Standard Error |
|---|---|
| 1000 | Terbesar |
| 5000 | Lebih kecil |
| 20000 | Terkecil |
Tabel ringkasan tersebut menunjukkan adanya hubungan yang konsisten antara peningkatan jumlah simulasi dengan peningkatan kualitas estimasi Monte Carlo. Secara umum, peningkatan jumlah simulasi meningkatkan peluang estimasi mendekati nilai teoritis walaupun pada satu realisasi simulasi tertentu galat tidak selalu menurun secara monoton.
One Sample t-test dilakukan untuk menguji apakah rata-rata hasil simulasi berbeda secara signifikan dari nilai ekspektasi teoritis sebesar 70 unit. Apabila diperoleh p-value lebih besar dari 0,05 maka tidak terdapat bukti yang cukup untuk menyatakan bahwa rata-rata hasil simulasi berbeda dari nilai teoritis. Dengan demikian, hasil simulasi dapat dianggap konsisten dengan nilai ekspektasi teoritis.
t1000 <- t.test(
sim_1000$data$Permintaan,
mu = ekspektasi_teoritis
)
t5000 <- t.test(
sim_5000$data$Permintaan,
mu = ekspektasi_teoritis
)
t20000 <- t.test(
sim_20000$data$Permintaan,
mu = ekspektasi_teoritis
)
hasil_ttest <- data.frame(
Simulasi = c(1000,5000,20000),
Mean = c(
sim_1000$mean,
sim_5000$mean,
sim_20000$mean
),
p_value = c(
t1000$p.value,
t5000$p.value,
t20000$p.value
)
)
knitr::kable(
hasil_ttest,
digits = 5,
caption="Hasil One Sample t-test"
)
| Simulasi | Mean | p_value |
|---|---|---|
| 1000 | 70.0700 | 0.83651 |
| 5000 | 69.9280 | 0.64232 |
| 20000 | 70.0235 | 0.76251 |
Interpretasi One Sample t-test
Berdasarkan hasil One Sample t-test pada Tabel 6.4, p-value pada simulasi 1.000, 5.000, dan 20.000 iterasi masing-masing sebesar 0,83651, 0,64232, dan 0,76251 — seluruhnya jauh lebih besar dari 0,05. Hal ini menunjukkan bahwa tidak terdapat bukti yang cukup untuk menyatakan bahwa rata-rata hasil simulasi berbeda secara signifikan dari nilai ekspektasi teoritis sebesar 70 unit. Dengan demikian, secara statistik, estimasi Monte Carlo pada ketiga skenario tetap konsisten dengan nilai ekspektasi teoritis, meskipun terdapat sedikit selisih numerik akibat variasi bilangan acak.
ANOVA digunakan untuk menguji apakah rata-rata permintaan hasil simulasi berbeda antar ukuran simulasi. Jika p-value lebih besar dari 0,05 maka tidak terdapat perbedaan rata-rata yang signifikan antara simulasi 1.000, 5.000, dan 20.000 iterasi.
anova_data <-
dplyr::bind_rows(
sim_1000$data %>%
dplyr::mutate(Simulasi="1000"),
sim_5000$data %>%
dplyr::mutate(Simulasi="5000"),
sim_20000$data %>%
dplyr::mutate(Simulasi="20000")
)
anova_model <-
aov(
Permintaan~Simulasi,
data=anova_data
)
summary(anova_model)
## Df Sum Sq Mean Sq F value Pr(>F)
## Simulasi 2 41 20.31 0.168 0.845
## Residuals 25997 3134358 120.57
Interpretasi Uji ANOVA
Hasil ANOVA menunjukkan nilai p-value yang lebih besar dari 0,05, sehingga tidak terdapat perbedaan rata-rata permintaan yang signifikan antara simulasi 1.000, 5.000, dan 20.000 iterasi. Hasil ini sejalan dengan temuan pada Bab 6.1–6.3, di mana seluruh rata-rata hasil simulasi berada sangat dekat dengan nilai ekspektasi teoritis. Tidak adanya perbedaan signifikan ini juga memperkuat bahwa peningkatan jumlah simulasi tidak mengubah karakteristik distribusi permintaan, melainkan hanya meningkatkan presisi (Standard Error) dari estimasi rata-rata tersebut.
bartlett_test <- bartlett.test(
Permintaan ~ Simulasi,
data = anova_data
)
bartlett_test
##
## Bartlett test of homogeneity of variances
##
## data: Permintaan by Simulasi
## Bartlett's K-squared = 1.2271, df = 2, p-value = 0.5414
Interpretasi Uji Homogenitas Varians
Hasil uji Bartlett menunjukkan nilai p-value lebih besar dari 0,05, sehingga tidak terdapat bukti yang cukup untuk menyatakan bahwa varians permintaan berbeda secara signifikan antara simulasi 1.000, 5.000, dan 20.000 iterasi. Hasil ini sejalan dengan nilai Standard Deviation pada Tabel 6.1 yang relatif sama pada ketiga simulasi, karena seluruh data dibangkitkan dari distribusi probabilitas yang identik — perbedaan hanya terletak pada jumlah iterasi, bukan pada karakteristik distribusinya.
Visualisasi merupakan salah satu tahap penting dalam analisis Simulasi Monte Carlo karena mampu menyajikan pola distribusi data, tingkat konvergensi, serta perubahan kualitas estimasi secara lebih intuitif dibandingkan penyajian numerik semata. Melalui grafik, hubungan antara jumlah simulasi dengan tingkat akurasi estimasi dapat diamati secara langsung sehingga mempermudah proses interpretasi.
Pada praktikum ini digunakan empat jenis visualisasi individual yang saling melengkapi, ditambah satu panel gabungan yang merangkum keempatnya. Masing-masing visualisasi dirancang untuk mengevaluasi aspek yang berbeda, mulai dari distribusi hasil simulasi, proses konvergensi rata-rata, tingkat kesalahan estimasi, hingga perbandingan Standard Error berdasarkan pendekatan analitik dan bootstrap. Seluruh visualisasi tersebut kemudian dirangkum dalam sebuah panel terpadu untuk memberikan gambaran menyeluruh mengenai performa Simulasi Monte Carlo.
Histogram digunakan untuk mengevaluasi apakah distribusi hasil simulasi telah mengikuti pola distribusi probabilitas yang dibangun dari data awal. Selain itu, visualisasi ini memperlihatkan bagaimana peningkatan jumlah simulasi memengaruhi kestabilan distribusi empiris terhadap distribusi teoritis.
library(dplyr)
library(ggplot2)
gabungan_data <-
bind_rows(
sim_1000$data %>%
mutate(Simulasi="N = 1000"),
sim_5000$data %>%
mutate(Simulasi="N = 5000"),
sim_20000$data %>%
mutate(Simulasi="N = 20000")
)
plot1 <-
ggplot(
gabungan_data,
aes(x=factor(Permintaan))
)+
geom_bar(fill="#2C7FB8")+
facet_wrap(~Simulasi)+
labs(
title="Distribusi Permintaan Hasil Simulasi",
x="Permintaan",
y="Frekuensi"
)+
theme_minimal()
plot1
Interpretasi Distribusi Permintaan
Histogram menunjukkan bahwa seluruh simulasi berhasil menghasilkan pola distribusi yang menyerupai distribusi probabilitas awal. Pada simulasi dengan 1.000 iterasi, masih terlihat fluktuasi frekuensi akibat pengaruh variasi acak. Ketika jumlah simulasi meningkat menjadi 5.000 dan 20.000, bentuk histogram semakin konsisten dengan distribusi teoritis. Nilai permintaan sebesar 70 unit tetap menjadi kategori yang paling sering muncul karena memiliki probabilitas terbesar (40%), sedangkan nilai 50 dan 90 unit muncul lebih sedikit sesuai probabilitasnya yang hanya sebesar 10%.
Hasil ini menunjukkan bahwa peningkatan jumlah simulasi mampu memperbaiki representasi distribusi empiris sehingga semakin mendekati distribusi probabilitas sebenarnya.
Running Mean digunakan untuk mengamati proses konvergensi rata-rata hasil simulasi menuju nilai ekspektasi teoritis. Grafik ini memperlihatkan perubahan estimasi rata-rata pada setiap iterasi simulasi.
plot2 <-
ggplot(
gabungan_data,
aes(
x=Iterasi,
y=Running_Mean
)
)+
geom_line(color="steelblue")+
facet_wrap(~Simulasi,
scales="free_x")+
geom_hline(
yintercept=70,
linetype="dashed",
color="red"
)+
labs(
title="Konvergensi Running Mean",
x="Iterasi",
y="Running Mean"
)+
theme_minimal()
plot2
Interpretasi Plot Konvergensi Running Mean
Pada awal simulasi, nilai rata-rata masih berfluktuasi cukup besar karena jumlah observasi yang masih sedikit sehingga setiap data baru memberikan pengaruh yang signifikan terhadap rata-rata. Namun, setelah jumlah iterasi terus bertambah, fluktuasi tersebut semakin berkurang dan kurva mulai mendekati garis horizontal yang menunjukkan nilai ekspektasi teoritis sebesar 70 unit.
Simulasi dengan 5.000 dan 20.000 iterasi menunjukkan konvergensi yang lebih stabil dibandingkan 1.000 iterasi.
Grafik ini digunakan untuk membandingkan rata-rata hasil simulasi dengan nilai ekspektasi teoritis sekaligus memperlihatkan besarnya galat estimasi pada masing-masing ukuran simulasi.
plot3 <-
ggplot(
hasil_perbandingan,
aes(
x=factor(Jumlah_Simulasi),
y=Mean
)
)+
geom_col(fill="#3182BD")+
geom_hline(
yintercept=70,
linetype="dashed",
color="red"
)+
geom_text(
aes(label=round(Mean,3)),
vjust=-0.4
)+
theme_minimal()
plot3
InterpretasI Estimasi Nilai Ekspektasi Hasil Simulasi
Grafik memperlihatkan bahwa rata-rata hasil simulasi pada ketiga ukuran simulasi berada sangat dekat dengan garis horizontal yang menunjukkan nilai ekspektasi teoritis sebesar 70 unit. Hal ini menunjukkan bahwa metode Monte Carlo mampu menghasilkan estimasi yang mendekati nilai ekspektasi teoritis.
Pada percobaan ini, simulasi 20.000 iterasi justru menghasilkan rata-rata yang paling dekat dengan nilai teoritis (Absolute Error 0,0235), diikuti oleh simulasi 1.000 iterasi (0,0700), sedangkan simulasi 5.000 iterasi memiliki Absolute Error sedikit lebih besar (0,0720) dibandingkan keduanya. Simulasi dengan 20.000 iterasi juga memberikan estimasi yang paling presisi karena memiliki Standard Error paling kecil. Dengan demikian, pada percobaan ini peningkatan jumlah simulasi menghasilkan perbaikan yang konsisten baik dari segi kedekatan rata-rata maupun presisi estimasi, sejalan dengan prediksi Law of Large Numbers.
Selain menggunakan pendekatan analitik, dilakukan pula estimasi Standard Error menggunakan metode Bootstrap. Perbandingan kedua pendekatan ini bertujuan untuk mengevaluasi konsistensi ukuran ketelitian estimasi. Bootstrap dilakukan sebanyak 1000 kali resampling untuk memperoleh estimasi Standard Error secara empiris tanpa bergantung pada rumus analitik.
library(boot)
boot_mean <- function(data, index){
mean(data[index])
}
boot1000 <- boot(
data = sim_1000$data$Permintaan,
statistic = boot_mean,
R = 1000
)
boot5000 <- boot(
data = sim_5000$data$Permintaan,
statistic = boot_mean,
R = 1000
)
boot20000 <- boot(
data = sim_20000$data$Permintaan,
statistic = boot_mean,
R = 1000
)
se_plot <- data.frame(
Simulasi = factor(
c("1000","5000","20000",
"1000","5000","20000"),
levels = c("1000","5000","20000")
),
Metode = c(
"Analitik","Analitik","Analitik",
"Bootstrap","Bootstrap","Bootstrap"
),
SE = c(
sim_1000$se,
sim_5000$se,
sim_20000$se,
sd(boot1000$t),
sd(boot5000$t),
sd(boot20000$t)
)
)
plot4 <- ggplot(
se_plot,
aes(x = Simulasi,
y = SE,
fill = Metode)
) +
geom_col(position = "dodge") +
labs(
title = "SE Analitik vs Bootstrap",
x = "Jumlah Simulasi",
y = "Standard Error"
) +
theme_minimal()
plot4
Interpretasi
Hasil menunjukkan bahwa Standard Error yang dihitung menggunakan pendekatan analitik dan bootstrap memiliki nilai yang relatif serupa. Kedekatan nilai Standard Error analitik dan bootstrap menunjukkan bahwa pendekatan analitik memberikan estimasi ketelitian yang konsisten dengan pendekatan resampling. Perbedaan yang kecil di antara keduanya merupakan konsekuensi alami dari proses bootstrap yang menggunakan pengambilan sampel ulang secara acak. Secara keseluruhan, hasil ini memperkuat bahwa estimasi yang diperoleh cukup stabil dan andal.
Untuk memberikan gambaran yang lebih komprehensif, seluruh grafik utama digabungkan menjadi satu panel visual menggunakan paket patchwork. Panel ini memudahkan pembaca dalam membandingkan distribusi data, proses konvergensi, kualitas estimasi, serta perubahan Standard Error secara bersamaan.
library(patchwork)
(plot1 | plot2) /
(plot3 | plot4)
Interpretasi
Panel visual memperlihatkan bahwa secara umum estimasi menjadi lebih stabil dan ditunjukkan melalui penurunan Standard Error, meskipun Absolute Error pada satu percobaan tidak selalu menurun secara monoton.
Secara keseluruhan, visualisasi mendukung hasil analisis numerik pada bab sebelumnya dan memperkuat bahwa Simulasi Monte Carlo bekerja semakin efektif ketika jumlah iterasi diperbesar.
Bab ini bertujuan untuk menginterpretasikan seluruh hasil simulasi Monte Carlo yang telah diperoleh pada bab sebelumnya. Analisis dilakukan dengan membandingkan hasil simulasi terhadap nilai ekspektasi teoritis, mengevaluasi perubahan tingkat kesalahan estimasi, mengamati proses konvergensi, serta menghubungkan hasil empiris dengan teori probabilitas dan statistika yang menjadi dasar metode Monte Carlo.
Berbeda dengan bab sebelumnya yang berfokus pada penyajian hasil numerik dan visualisasi, pada bab ini pembahasan diarahkan untuk menjelaskan mengapa hasil tersebut terjadi serta bagaimana hubungan antara ukuran simulasi dengan kualitas estimasi yang dihasilkan.
Berdasarkan hasil simulasi, seluruh estimasi rata-rata berhasil mendekati nilai ekspektasi teoritis sebesar 70 unit. Selisih antara rata-rata simulasi dengan nilai teoritis sangat kecil sehingga menunjukkan bahwa metode Monte Carlo mampu menghasilkan estimasi yang akurat.
Peningkatan jumlah simulasi memberikan dampak yang jelas terhadap tingkat presisi estimasi. Hal ini ditunjukkan oleh penurunan Standard Error dari 0,3391 pada 1.000 iterasi menjadi 0,1550 pada 5.000 iterasi, kemudian 0,0778 pada 20.000 iterasi. Selain itu, interval kepercayaan juga semakin sempit sehingga estimasi rata-rata menjadi semakin stabil.
Pada percobaan ini, nilai Absolute Error terkecil diperoleh pada simulasi 20.000 iterasi (0,0235), sedangkan simulasi 5.000 iterasi menghasilkan galat yang sedikit lebih besar dibandingkan simulasi 1.000 iterasi (0,0720 berbanding 0,0700). Kondisi tersebut merupakan karakteristik alami simulasi Monte Carlo karena setiap percobaan menggunakan bilangan acak yang berbeda. Pada kasus ini, hasil yang diperoleh justru selaras dengan ekspektasi teori: Standard Error, Absolute Error, dan kestabilan Running Mean seluruhnya menunjukkan perbaikan yang konsisten seiring bertambahnya jumlah simulasi.
Hasil simulasi yang diperoleh pada praktikum ini sejalan dengan teori dasar Simulasi Monte Carlo. Metode Monte Carlo pada dasarnya memanfaatkan pembangkitan bilangan acak untuk menghasilkan sampel yang mengikuti distribusi probabilitas tertentu. Karena setiap sampel dibentuk secara acak, hasil simulasi pada ukuran sampel kecil masih mengandung variasi yang cukup besar. Variasi tersebut merupakan konsekuensi alami dari proses pengambilan sampel acak.
Ketika jumlah simulasi diperbesar, pengaruh variasi acak terhadap rata-rata sampel menjadi semakin kecil. Fenomena ini dijelaskan oleh Hukum Bilangan Besar (Law of Large Numbers) yang menyatakan bahwa rata-rata sampel akan mengalami konvergensi menuju nilai harapan populasi apabila ukuran sampel terus bertambah. Oleh karena itu, peningkatan jumlah simulasi dari 1.000 menjadi 20.000 iterasi menghasilkan estimasi yang semakin mendekati nilai ekspektasi teoritis.
Selain Hukum Bilangan Besar, hasil simulasi juga sesuai dengan Teorema Limit Pusat (Central Limit Theorem). Teorema ini menjelaskan bahwa distribusi rata-rata sampel akan semakin mendekati distribusi normal ketika ukuran sampel cukup besar.
Dari sudut pandang statistik komputasi, peningkatan jumlah simulasi juga menyebabkan Monte Carlo Error mengalami penurunannya mengikuti orde 1/√N. Secara teoritis, besar kesalahan Monte Carlo berbanding terbalik dengan akar kuadrat jumlah simulasi, atau dapat dinyatakan sebagai:
\[ \text{Standard Error = Monte Carlo Error} \propto \frac{1}{\sqrt{N}} \] Persamaan tersebut menunjukkan bahwa peningkatan jumlah simulasi akan mengurangi tingkat kesalahan estimasi. Namun, penurunan kesalahan tidak berlangsung secara linear. Sebagai contoh, untuk mengurangi kesalahan hingga setengahnya, jumlah simulasi harus ditingkatkan sekitar empat kali lipat.
Dalam konteks praktikum ini, 20.000 iterasi memberikan estimasi yang paling presisi karena memiliki Standard Error terkecil, dan pada percobaan ini rata-ratanya juga merupakan yang paling dekat dengan nilai teoritis (Absolute Error 0,0235) dibandingkan kedua skenario lainnya.
Hasil pengujian One Sample t-test menunjukkan bahwa seluruh rata-rata hasil simulasi tidak berbeda signifikan dari nilai ekspektasi teoritis (p-value > 0,05 pada ketiga skenario), sehingga memperkuat secara formal bahwa metode Monte Carlo menghasilkan estimasi yang konsisten dengan nilai teoritis.
Hasil ANOVA juga menunjukkan tidak adanya perbedaan rata-rata yang signifikan antara ketiga ukuran simulasi (p-value > 0,05), yang mengindikasikan bahwa peningkatan jumlah iterasi tidak mengubah nilai harapan distribusi permintaan, melainkan hanya memperbaiki presisi estimasi (ditunjukkan oleh penurunan Standard Error pada Bab 6).
Selanjutnya, uji homogenitas varians menggunakan Bartlett’s test menunjukkan bahwa varians ketiga simulasi juga tidak berbeda signifikan (p-value > 0,05), sesuai ekspektasi karena seluruh data dibangkitkan dari distribusi probabilitas yang identik. Temuan ini memperkuat asumsi dasar ANOVA bahwa varians antar kelompok bersifat homogen, sehingga hasil uji ANOVA pada bagian sebelumnya dapat dipercaya validitasnya.
Berdasarkan hasil implementasi dan analisis yang telah dilakukan, dapat disimpulkan bahwa metode Simulasi Monte Carlo mampu mengestimasi nilai ekspektasi distribusi permintaan dengan tingkat akurasi yang tinggi. Seluruh hasil simulasi menghasilkan rata-rata yang sangat dekat dengan nilai ekspektasi teoritis sebesar 70 unit, sehingga menunjukkan bahwa metode Monte Carlo mampu menghasilkan estimasi nilai ekspektasi yang konsisten dengan distribusi probabilitas yang digunakan.
Peningkatan jumlah simulasi dari 1.000 menjadi 5.000 dan 20.000 iterasi memberikan pengaruh yang signifikan terhadap kualitas estimasi. Semakin besar jumlah simulasi, semakin kecil nilai Standard Error dan semakin sempit interval kepercayaan, sehingga estimasi rata-rata menjadi semakin presisi dan memiliki tingkat ketidakpastian yang lebih rendah. Walaupun Absolute Error tidak selalu menurun secara monoton pada setiap percobaan akibat variasi acak, secara keseluruhan hasil simulasi tetap sangat dekat dengan nilai ekspektasi teoritis.
Hasil praktikum ini juga memberikan bukti empiris yang konsisten dengan Hukum Bilangan Besar (Law of Large Numbers), yaitu rata-rata hasil simulasi menunjukkan proses konvergensi menuju nilai ekspektasi teoritis seiring bertambahnya jumlah iterasi. Selain itu, hasil yang diperoleh juga sejalan dengan Teorema Limit Pusat (Central Limit Theorem), yang menjelaskan bahwa distribusi rata-rata sampel menjadi semakin stabil ketika ukuran sampel diperbesar.
Hasil pengujian statistik inferensial (One Sample t-test, ANOVA, dan uji homogenitas varians) memperkuat temuan deskriptif di atas. Ketiga simulasi tidak menunjukkan perbedaan rata-rata yang signifikan terhadap nilai ekspektasi teoritis maupun antar sesama simulasi, dan varians antar kelompok juga homogen. Hal ini membuktikan secara statistik bahwa peningkatan jumlah simulasi tidak mengubah nilai harapan maupun keragaman data, melainkan hanya meningkatkan presisi estimasi melalui penurunan Standard Error.
Secara keseluruhan, Simulasi Monte Carlo merupakan metode yang efektif untuk mengestimasi parameter statistik, terutama ketika penyelesaian analitik sulit dilakukan. Meskipun peningkatan jumlah simulasi memerlukan waktu komputasi yang lebih besar, peningkatan tersebut menghasilkan estimasi yang lebih presisi dan andal sehingga metode ini layak digunakan pada berbagai permasalahan statistika, optimasi, maupun pengambilan keputusan di bawah kondisi ketidakpastian.
Sanusi, R. N. M. (2025). Simulasi Monte Carlo. Modul Mata Kuliah Pemodelan Statistika dan Simulasi.
Fishman, G. S. (1996). Monte Carlo: Concepts, Algorithms, and Applications. Springer.
Glasserman, P. (2004). Monte Carlo Methods in Financial Engineering. Springer.
Robert, C. P., & Casella, G. (2004). Monte Carlo Statistical Methods (2nd ed.). Springer.
Rubinstein, R. Y., & Kroese, D. P. (2017). Simulation and the Monte Carlo Method (3rd ed.). Wiley.
Ross, S. M. (2020). Introduction to Probability Models (12th ed.). Academic Press.
Efron, B., & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall/CRC.
Law, A. M. (2015). Simulation Modeling and Analysis (5th ed.). McGraw-Hill Education.
Devore, J. L. (2016). Probability and Statistics for Engineering and the Sciences (9th ed.). Cengage Learning.
Kroese, D. P., Taimre, T., & Botev, Z. I. (2011). Handbook of Monte Carlo Methods. Wiley.