Simulasi Monte Carlo

Topik ini menerapkan Generator NumPy pada simulasi Monte Carlo, yaitu teknik estimasi nilai numerik melalui pengambilan sampel acak yang berulang. Dua simulasi dibahas, yaitu estimasi nilai pi dari titik acak serta estimasi peluang dan nilai harapan dari pelemparan koin dan dadu, disertai perilaku errornya saat jumlah sampel diperbesar.

Ide dasar simulasi Monte Carlo#

Simulasi Monte Carlo mengestimasi besaran yang sulit dihitung langsung dengan cara mengambil sampel acak dalam jumlah besar, lalu meringkas sampel tersebut menjadi satu angka. Peluang sebuah kejadian diestimasi dari frekuensi kemunculannya, sedangkan nilai harapan diestimasi dari rata-rata sampel. Ketepatan estimasi bergantung pada jumlah sampel: semakin banyak sampel, ringkasannya semakin dekat ke nilai sejatinya. Keunggulan teknik ini muncul pada masalah yang sulit dihitung dengan rumus tertutup, misalnya integral berdimensi tinggi atau sistem dengan banyak kemungkinan, karena biayanya tumbuh jauh lebih lambat dibanding pemeriksaan seluruh kemungkinan. Harga yang dibayarkan adalah fluktuasi, sehingga setiap keluaran simulasi harus dibaca sebagai estimasi, bukan sebagai nilai pasti.

Contoh yang digunakan di sini adalah estimasi pi. Titik acak dibangkitkan pada persegi satu satuan, lalu dihitung proporsi titik yang jatuh di dalam seperempat lingkaran berjari-jari satu. Perbandingan luas keduanya menyatakan hubungan berikut.

Proporsi titik yang berada di dalam seperempat lingkaran, dikalikan empat, karenanya menjadi estimasi untuk pi.

Mengestimasi pi dengan titik acak#

Blok berikut membangkitkan sejuta pasang koordinat acak pada persegi $[0, 1) \times [0, 1)$. Sebuah titik berada di dalam seperempat lingkaran yang berpusat di titik asal dengan jari-jari satu apabila $x^2 + y^2 \le 1$. Pendekatan ini disebut hit-or-miss karena setiap titik hanya menerima dua kemungkinan, yaitu di dalam atau di luar lingkaran. Estimasi pi diperoleh dari proporsi titik yang memenuhi syarat tersebut.

import numpy as np

rng = np.random.default_rng(seed=0)

n = 1_000_000
x = rng.random(n)
y = rng.random(n)
di_dalam = np.count_nonzero(x**2 + y**2 <= 1.0)
estimasi = 4.0 * di_dalam / n

print(f"titik di dalam seperempat lingkaran : {di_dalam}")
print(f"estimasi pi                         : {estimasi:.6f}")
print(f"nilai np.pi                         : {np.pi:.6f}")
print(f"selisih mutlak                      : {abs(estimasi - np.pi):.6f}")

titik di dalam seperempat lingkaran : 785342
estimasi pi                         : 3.141368
nilai np.pi                         : 3.141593
selisih mutlak                      : 0.000225

Dari sejuta titik, 785.342 titik jatuh di dalam seperempat lingkaran, sehingga estimasinya 3.141368 dengan selisih 0.000225 terhadap np.pi. Seluruh proses berjalan tanpa perulangan Python yang eksplisit: rng.random(n) membangkitkan seluruh koordinat sebagai array, ekspresi x**2 + y**2 <= 1.0 menghasilkan mask boolean, dan np.count_nonzero menghitung titik yang lulus seleksi. Angka pada keluaran tetap sama setiap blok dijalankan karena Generatornya berseed. Nilai estimasi ini bukan nilai pi yang pasti, melainkan satu tarikan dari estimasi yang berfluktuasi mengikuti sampel yang ditarik.

Konvergensi seiring bertambahnya sampel#

Blok berikut mengulang estimasi yang sama untuk beberapa ukuran sampel, dari seratus hingga sepuluh juta titik, dengan satu Generator yang berlanjut antar ukuran.

import numpy as np

rng = np.random.default_rng(seed=0)

for n in [100, 10_000, 1_000_000, 10_000_000]:
    x = rng.random(n)
    y = rng.random(n)
    estimasi = 4.0 * np.count_nonzero(x**2 + y**2 <= 1.0) / n
    print(f"n = {n:>9}: estimasi = {estimasi:.6f}, selisih = {abs(estimasi - np.pi):.6f}")

n =       100: estimasi = 2.720000, selisih = 0.421593
n =     10000: estimasi = 3.144800, selisih = 0.003207
n =   1000000: estimasi = 3.141408, selisih = 0.000185
n =  10000000: estimasi = 3.142120, selisih = 0.000528

Error menurun drastis dari 0.421593 pada seratus titik menjadi 0.003207 pada sepuluh ribu titik. Setelah itu penurunannya melambat: error 0.000185 pada sejuta titik justru naik menjadi 0.000528 pada sepuluh juta titik. Perilaku ini wajar karena error Monte Carlo menurun secara rata-rata sebanding dengan $1/\sqrt{N}$, bukan menurun secara monoton pada setiap eksperimen tunggal. Angka pada keluaran adalah satu contoh tarikan, sehingga eksperimen dengan seed lain menghasilkan angka yang berbeda di setiap baris, tetapi pola menurunnya tetap sama. Mempersepuluh kali jumlah sampel hanya mengecilkan error rata-rata sekitar tiga kali, sehingga biaya penambahan sampel semakin berat dibanding manfaatnya.

Simulasi pelemparan koin dan dadu#

Prinsip yang sama berlaku pada simulasi sederhana. Proporsi muka pada pelemparan koin mengestimasi peluang muka, sedangkan rata-rata hasil pelemparan dadu mengestimasi nilai harapannya. Pada koin, angka 0 dan 1 mewakili dua sisi, sehingga rata-rata array yang berisi 0 dan 1 langsung menjadi proporsi muka. Pada dadu, integers(1, 7) memberikan keenam muka dengan peluang yang sama, sehingga rata-ratanya mendekati nilai harapan teoretis. Blok berikut menjalankan keduanya dengan sejuta pengulangan.

import numpy as np

rng = np.random.default_rng(seed=0)

lemparan = rng.integers(low=0, high=2, size=1_000_000)
proporsi_muka = np.mean(lemparan)

dadu = rng.integers(low=1, high=7, size=1_000_000)
rata_rata = np.mean(dadu)

print(f"proporsi muka  : {proporsi_muka:.6f} (nilai teoretis 0.5)")
print(f"rata-rata dadu : {rata_rata:.6f} (nilai teoretis 3.5)")

proporsi muka  : 0.500418 (nilai teoretis 0.5)
rata-rata dadu : 3.500266 (nilai teoretis 3.5)

Proporsi muka tercatat 0.500418, berdekatan dengan peluang teoretis 0.5. Rata-rata dadu tercatat 3.500266, berdekatan dengan nilai harapan $(1 + 2 + 3 + 4 + 5 + 6) / 6 = 3.5$. Kedua hasil menunjukkan pola yang sama dengan estimasi pi: angka simulasi mendekati nilai teoretis, tetapi tidak pernah persis sama. Perbedaan kecil yang tersisa mengecil pada ukuran sampel yang lebih besar, dengan laju penurunan yang sama seperti error estimasi pi.

Kesalahan umum#

Menganggap satu keluaran simulasi sebagai nilai pasti. Estimasi Monte Carlo selalu disertai fluktuasi, sehingga dua eksperimen dengan ukuran sampel sama dapat menghasilkan error yang berbeda. Ukuran sampel perlu dilaporkan bersama estimasinya karena pembaca memerlukannya untuk menilai seberapa besar fluktuasi yang wajar. Bandingkan hasil dengan nilai teoretis apabila nilai tersebut diketahui.

Membandingkan hasil float dengan kesamaan sempurna. Estimasi pi hampir tidak pernah persis sama dengan np.pi, sehingga perbandingan == bernilai False meskipun estimasinya sudah sangat dekat. Gunakan selisih mutlak dengan toleransi, sebagaimana blok berikut. Dengan seed yang tercatat, hasil yang sama dapat diulang kapan pun perbandingannya perlu diperiksa ulang.

import numpy as np

rng = np.random.default_rng(seed=0)

n = 1_000_000
x = rng.random(n)
y = rng.random(n)
estimasi = 4.0 * np.count_nonzero(x**2 + y**2 <= 1.0) / n

print(estimasi == np.pi)
print(abs(estimasi - np.pi) < 0.001)

False
True

Menjalankan ulang simulasi tanpa seed saat hasil perlu dibandingkan. Tanpa seed, setiap eksekusi menghasilkan titik yang berbeda, sehingga perubahan error antar eksekusi tidak dapat diatribusikan ke perubahan kode. Isi seed sebelum melakukan perbandingan antar versi program. Kebiasaan ini juga berlaku saat keluaran simulasi dimasukkan ke laporan yang dibangkitkan secara otomatis.

Ringkasan#

• Monte Carlo mengestimasi peluang dan nilai harapan dari frekuensi serta rata-rata sampel acak.
• Estimasi pi diperoleh dari empat kali proporsi titik acak yang jatuh di dalam seperempat lingkaran.
• Error menurun secara rata-rata sebanding dengan $1/\sqrt{N}$, bukan monoton pada setiap eksperimen tunggal.
• Laporkan estimasi dengan ukuran sampel dan toleransi selisihnya, bukan dengan perbandingan kesamaan float.

Referensi#

• Random sampling (numpy.random) — NumPy — Ikhtisar submodul numpy.random yang menjadi dasar seluruh contoh.
• numpy.random.Generator — NumPy — Dokumentasi metode Generator yang digunakan untuk membangkitkan sampel.