Bagian 1: Membangun Model Probabilitas

Sebelum satu angka acak pun dibangkitkan, sebuah proyek Monte Carlo dimulai dari pekerjaan paling tenang: merumuskan pertanyaan menjadi model probabilitas. Model itu menjawab tiga hal — besaran acak apa yang kita bangkitkan, aturan apa yang menghubungkan angka acak dengan jawaban, dan berapa nilai teoretisnya bila dihitung dengan rumus. Proyek ini sengaja memilih tiga persoalan yang jawaban eksaknya sudah diketahui: estimasi pi, peluang jumlah dua dadu sama dengan 7, dan integral sin(x) pada [0, π]. Dengan nilai teoretis di tangan, galat setiap simulasi pada bagian-bagian berikutnya bisa diukur, bukan ditebak.

Peta ketiga model yang akan kita bangun:

Model 1: Peluang Jumlah Dua Dadu #

Ruang sampel dua dadu berisi 36 kombinasi yang masing-masing sama mungkinnya. Jumlah 7 bisa terwujud lewat enam kombinasi — (1,6), (2,5), (3,4), (4,3), (5,2), dan (6,1) — sehingga peluangnya 6/36 = 1/6. Kita bisa membiarkan NumPy menghitung enumerasi itu tanpa satu pun loop for, hanya dengan broadcasting:

import numpy as np

dadu1 = np.arange(1, 7)
dadu2 = np.arange(1, 7)
jumlah = dadu1[:, None] + dadu2[None, :]
print(jumlah)
print("Banyak kombinasi jumlah = 7 :", np.count_nonzero(jumlah == 7))
print("Total kemungkinan           :", jumlah.size)
print("P(jumlah = 7) = 6/36        =", 6 / 36)

Output yang Dihasilkan:

[[ 2  3  4  5  6  7]
 [ 3  4  5  6  7  8]
 [ 4  5  6  7  8  9]
 [ 5  6  7  8  9 10]
 [ 6  7  8  9 10 11]
 [ 7  8  9 10 11 12]]
Banyak kombinasi jumlah = 7 : 6
Total kemungkinan           : 36
P(jumlah = 7) = 6/36        = 0.16666666666666666

Mari kita bedah:

1. Baris 3-4 (np.arange(1, 7)): dua array 1-D berisi muka dadu 1 sampai 6; batas atas 7 tidak disertakan.
2. Baris 5 (dadu1[:, None] + dadu2[None, :]): [:, None] menjadikan dadu1 kolom dan [None, :] menjadikan dadu2 baris, sehingga penjumlahan broadcasting menghasilkan matriks 6×6; sel [i, j] menyimpan jumlah muka dadu1[i] + dadu2[j].
3. Baris 7 (np.count_nonzero(jumlah == 7)): perbandingan array menghasilkan mask boolean berukuran 6×6, dan count_nonzero menghitung sel True — tepat 6 kombinasi penyusun jumlah 7.
4. Baris 9: peluang eksak 6/36 dicetak sebagai desimal 0.16666666666666666; pecahan 1/6 memang tidak pernah habis dibagi, dan float64 menyimpan nilai terdekatnya.

Model 2: Rasio Luas Lingkaran dalam Persegi #

Estimasi pi lahir dari perbandingan dua luas. Lingkaran berjari-jari 1 beralaskan persegi bersisi 2: luas lingkaran adalah π·r² = π dan luas persegi adalah 4, sehingga rasio keduanya adalah π/4 ≈ 0.7854. Jika kita menebar titik acak merata di persegi itu, proporsi titik yang jatuh di dalam lingkaran akan mendekati rasio tersebut — dan mengalikannya dengan 4 mengembalikan pi:

import numpy as np

r = 1.0
luas_persegi = (2 * r) ** 2
luas_lingkaran = np.pi * r**2
print("Luas lingkaran berjari-jari 1   :", luas_lingkaran)
print("Luas persegi sisi 2             :", luas_persegi)
print("Rasio luas (target Monte Carlo) :", luas_lingkaran / luas_persegi)
print("Nilai np.pi / 4                 :", np.pi / 4)

Output yang Dihasilkan:

Luas lingkaran berjari-jari 1   : 3.141592653589793
Luas persegi sisi 2             : 4.0
Rasio luas (target Monte Carlo) : 0.7853981633974483
Nilai np.pi / 4                 : 0.7853981633974483

Penjelasan tiap baris:

1. Baris 3 (luas_persegi = (2 * r) ** 2): persegi yang memayungi lingkaran berjari-jari 1 punya sisi 2·r, maka luasnya 4.0.
2. Baris 4 (luas_lingkaran = np.pi * r**2): rumus luas lingkaran memakai konstanta np.pi; hasilnya 3.141592653589793, nilai yang nanti kita kejar lewat sampel acak.
3. Baris 8-9: rasio luas dan np.pi / 4 mencetak nilai yang sama persis — inilah besaran yang diestimasi simulasi pada Bagian 3, yaitu proporsi titik acak yang berada di dalam lingkaran.

Model 3: Integral sebagai Rata-rata Nilai Fungsi #

Model ketiga memanfaatkan sebuah identitas: bila x disebar uniform pada [a, b], rata-rata f(x) mendekati nilai rata-rata integral, sehingga ∫ f(x) dx ≈ (b − a) × mean(f(x)). Untuk f(x) = sin(x) pada [0, π], nilai eksaknya bisa dihitung analitik: hasilnya 2.0. Sebagai alat kontrol — belum dengan angka acak — kita uji identitas itu pada grid tetap yang sangat rapat:

import numpy as np

a, b = 0.0, np.pi
x = np.linspace(a, b, 1_000_001)
rata_rata_sin = np.sin(x).mean()
print("Rata-rata sin(x) pada grid rapat :", rata_rata_sin)
print("Estimasi integral (b-a) * mean   :", (b - a) * rata_rata_sin)
print("Nilai eksak analitik             :", 2.0)

Output yang Dihasilkan:

Rata-rata sin(x) pada grid rapat : 0.6366191357479217
Estimasi integral (b-a) * mean   : 1.999998000000354
Nilai eksak analitik             : 2.0

Mari kita baca hasilnya:

1. Baris 3-4 (np.linspace(a, b, 1_000_001)): baris 3 menetapkan batas integral a = 0 dan b = np.pi, lalu baris 4 menyusun grid berisi satu juta titik berjarak sama dari 0 sampai π; grid ini deterministik, bukan sampel acak.
2. Baris 5: rata-rata sin(x) pada grid menghasilkan 0.6366191357479217, mendekati nilai analitik 2/π ≈ 0.6366197723.
3. Baris 7: mengalikan dengan (b − a) = π menghasilkan 1.999998000000354 — meleset sekitar 0.000002 dari 2.0 karena grid berhingga tidak pernah mencapai kepadatan kontinu.

Ketiga model kini siap: masing-masing punya besaran acak, estimator, dan nilai teoretis sebagai patokan. Bagian berikutnya membangkitkan sampel acak yang membuat ketiga estimator itu bergerak.