Bagian 3: Menjalankan Simulasi

Sampel kecil pada Bagian 2 membuktikan bahwa mesinnya bekerja; kini ukuran sampel dinaikkan ke n = 100.000 dan ketiga estimator dieksekusi penuh. Seluruh eksperimen memakai seed 42 sehingga setiap angka di bawah ini direproduksi persis. Seandainya proyek ini disimpan sebagai folder kode, ketiga estimator akan kita letakkan sebagai fungsi di src/model.py, sementara pemanggilnya dengan berbagai ukuran sampel di src/simulasi.py — model dan eksperimen dapat diuji terpisah. Di bab ini keduanya ditulis inline agar setiap blok tetap berdiri sendiri.

Estimasi pi pada 100.000 Titik #

Model lingkaran pada Bagian 1 kini diberi 100.000 titik acak. Berbeda dengan contoh ringkas pada Chapter Simulasi Monte Carlo yang memakai seperempat lingkaran di persegi satuan, proyek ini memakai lingkaran penuh berjari-jari 1 di persegi sisi 2; rasio yang dikejar sama, π/4. Koordinat dihasilkan rng.uniform(-1.0, 1.0) — bentuk uniform(low, high) yang dibahas pada Chapter Uniform Distribution:

import numpy as np

rng = np.random.default_rng(42)
n = 100_000
x = rng.uniform(-1.0, 1.0, size=n)
y = rng.uniform(-1.0, 1.0, size=n)
di_lingkaran = np.count_nonzero(x**2 + y**2 <= 1.0)
estimasi = 4 * di_lingkaran / n
print("Titik di dalam lingkaran :", di_lingkaran)
print("Proporsi                 :", di_lingkaran / n)
print(f"Estimasi pi (4 x proporsi) : {estimasi:.4f}")
print(f"Nilai np.pi                : {np.pi:.4f}")

Output yang Dihasilkan:

Titik di dalam lingkaran : 78577
Proporsi                 : 0.78577
Estimasi pi (4 x proporsi) : 3.1431
Nilai np.pi                : 3.1416

Mari kita bedah:

1. Baris 5-6 (rng.uniform(-1.0, 1.0, size=n)): dua array berisi 100.000 koordinat uniform pada [-1, 1), mewakili 100.000 titik (x, y) yang tersebar merata di persegi sisi 2.
2. Baris 7 (np.count_nonzero(x**2 + y**2 <= 1.0)): uji jarak ke pusat — titik berada di dalam lingkaran bila x² + y² ≤ 1; mask boolean yang lulus uji dihitung, hasilnya 78.577 titik.
3. Baris 8 (estimasi = 4 * di_lingkaran / n): proporsi 0.78577 dikalikan 4, menghasilkan estimasi pi = 3.1431 — meleset 0.0015 dari 3.1416, atau sekitar 0.05%.
4. Baris 11-12: estimator dan nilai teoretis dicetak dengan format empat desimal agar perbandingannya berdampingan.

Peluang Empiris Jumlah Dadu 7 #

Eksperimen kedua meniru pola Bagian 2 pada skala penuh: 100.000 lemparan dua dadu, lalu proporsi lemparan berjumlah 7 dibandingkan dengan teori 6/36:

import numpy as np

rng = np.random.default_rng(42)
n = 100_000
dadu = rng.integers(1, 7, size=(n, 2))
jumlah = dadu.sum(axis=1)
frekuensi_7 = np.count_nonzero(jumlah == 7)
proporsi = frekuensi_7 / n
print("Frekuensi jumlah = 7 :", frekuensi_7)
print(f"Proporsi empiris     : {proporsi:.4f}")
print(f"Teori 6/36           : {6 / 36:.4f}")

Output yang Dihasilkan:

Frekuensi jumlah = 7 : 16827
Proporsi empiris     : 0.1683
Teori 6/36           : 0.1667

Penjelasan tiap bagian:

1. Baris 5 (size=(n, 2)): satu panggilan membangkitkan 100.000 lemparan berpasangan; setiap baris adalah satu lemparan dua dadu.
2. Baris 6-7: sum(axis=1) merangkum tiap baris menjadi jumlah dua dadu, dan count_nonzero menghitung yang berjumlah 7 — tepat 16.827 lemparan.
3. Baris 9-11: proporsi empiris 0.1683 berbanding teori 0.1667; selisih +0.0016, dan presisinya diukur pada Bagian 4.

Integral sin(x) lewat Rata-rata Sampel #

Eksperimen ketiga memakai identitas Bagian 1: ∫ f(x) dx ≈ (b − a) × mean(f(x)) untuk x uniform pada [a, b]. Karena rng.random() menyuplai bilangan uniform pada [0, 1), nilai x pada [0, π) diperoleh dengan transformasi sederhana rng.random(n) * np.pi:

import numpy as np

rng = np.random.default_rng(42)
n = 100_000
x = rng.random(n) * np.pi
rata_rata = np.sin(x).mean()
estimasi = np.pi * rata_rata
print(f"Rata-rata sin(x)  : {rata_rata:.6f}")
print(f"Estimasi integral : {estimasi:.4f}")
print("Nilai eksak teori : 2.0000")

Output yang Dihasilkan:

Rata-rata sin(x)  : 0.637035
Estimasi integral : 2.0013
Nilai eksak teori : 2.0000

Mari kita telaah:

1. Baris 5 (x = rng.random(n) * np.pi): 100.000 bilangan uniform [0, 1) digeser skalanya menjadi [0, π); transformasi linear ini menjaga kemerataan sebarannya.
2. Baris 6 (np.sin(x).mean()): sin dievaluasi pada seluruh titik sekaligus (vectorized), lalu dirata-rata — hasilnya 0.637035, mendekati nilai analitik 2/π ≈ 0.6366198.
3. Baris 7: mengalikan dengan (b − a) = π menghasilkan 2.0013 — meleset 0.0013 dari nilai eksak 2.0.

Konvergensi pada Empat Ukuran Sampel #

Satu angka besar belum menceritakan perilaku estimator; kita ulangi estimasi pi pada empat ukuran sampel — 100, 1.000, 10.000, dan 100.000 titik — dengan fungsi yang membentuk Generator baru ber-seed 42 untuk setiap ukuran:

import numpy as np

def estimasi_pi(n, seed=42):
    rng = np.random.default_rng(seed)
    x = rng.uniform(-1.0, 1.0, size=n)
    y = rng.uniform(-1.0, 1.0, size=n)
    return 4 * np.count_nonzero(x**2 + y**2 <= 1.0) / n

for n in [100, 1_000, 10_000, 100_000]:
    print(f"n = {n:>7} : estimasi pi = {estimasi_pi(n):.4f}")

Output yang Dihasilkan:

n =     100 : estimasi pi = 3.1200
n =    1000 : estimasi pi = 3.0800
n =   10000 : estimasi pi = 3.1256
n =  100000 : estimasi pi = 3.1431

Bedah baris demi baris:

1. Baris 3 (def estimasi_pi(n, seed=42)): seluruh langkah estimasi dikemas sebagai fungsi; Generator dibentuk ulang di dalam fungsi, sehingga setiap ukuran sampel memakai deret acak tersendiri.
2. Baris 5-6: urutan pembangkitan identik dengan blok pertama — x dahulu, lalu y; karena itu baris n = 100000 mencetak 3.1431, persis hasil eksperimen pertama.
3. Baris 9-10: loop memanggil fungsi pada empat skala; format {n:>7} menjaga kolom output tetap sejajar.

Dibandingkan dengan pi = 3.141593, galat tiap skala:

Mengukur Waktu Eksekusi #

Sebelum menutup bagian ini, wajar bertanya berapa biaya 100.000 titik tersebut. Pengukuran waktu memakai time.perf_counter:

# notest
import numpy as np
import time

rng = np.random.default_rng(42)
n = 100_000
mulai = time.perf_counter()
x = rng.uniform(-1.0, 1.0, size=n)
y = rng.uniform(-1.0, 1.0, size=n)
estimasi = 4 * np.count_nonzero(x**2 + y**2 <= 1.0) / n
durasi = time.perf_counter() - mulai
print(f"Estimasi pi = {estimasi:.4f} selesai dalam {durasi * 1000:.1f} ms")

Pada mesin uji proyek ini, satu siklus estimasi pi dengan 100.000 titik berlangsung dalam orde dua milidetik. Angka waktu berbeda di tiap komputer, karena itu blok di atas tidak dipakai untuk membandingkan output.

Ketiga eksperimen kini punya angka; pertanyaan berikutnya: seberapa dekat dan seberapa yakin? Itulah pekerjaan Bagian 4.