14 Rastgele Sayılar ve Monte Carlo Yöntemleri
İki zar atıldığında toplamın 7 gelme olasılığı kaçtır? Sayma teknikleriyle (bkz. Olasılık Teorisi) eşit olasılıklı \(36\) sonucun \(6\)’sının toplamı 7 olduğunu görür ve \(1/6\) buluruz. Aynı soruyu bilgisayara deneyi yüz bin kez yaptırarak da sorabiliriz:
import numpy as np
rng = np.random.default_rng(1)
n = 100_000
d1 = rng.integers(1, 7, size=n) # birinci zar: 1, 2, ..., 6
d2 = rng.integers(1, 7, size=n) # ikinci zar
seven = d1 + d2 == 7 # toplam 7 mi?
print(d1[:8])
print(d2[:8])
print(seven.sum(), seven.mean(), 1 / 6)Çıktı:
[3 4 5 6 1 1 5 6]
[5 4 6 3 5 4 1 3]
16450 0.1645 0.16666666666666666rng.integers(1, 7, size=n) 1 ile 6 arasında \(n\) tane tamsayı üretir; üst sınır olan 7 dahil değildir. d1 + d2 == 7 bir boolean dizidir: sum doğru olan elemanları sayar, mean ise doğruların oranını verir, çünkü hesapta True 1, False 0 sayılır. Yüz bin denemenin \(16450\)’sinde toplam 7 geldi; göreli sıklık \(0{,}1645\), gerçek olasılık ise \(1/6 \approx 0{,}1667\)’dir.
Bir olasılığı, bir beklenen değeri ya da bir integrali, ilgili rastgele deneyi çok kez tekrarlayıp sonuçların ortalamasını alarak kestirme yöntemine Monte Carlo yöntemi denir; adını kumarhaneleriyle ünlü Monte Carlo şehrinden alır. Vektörizasyon ve Broadcasting bölümünde \(\pi\) sayısını rastgele noktalarla kestirirken bu yöntemi bir kez kullanmıştık. Kestirim tam değildir: kodu başka bir tohumla çalıştırınca başka bir sayı çıkar. Bu bölüm iki soruyu birlikte ele alıyor: rastgele sayıları bilgisayarda nasıl üretiriz ve onlarla yaptığımız kestirimlerin hatası ne kadardır? Önce NumPy’ın üretecini ve temel dağılımları tanıyacağız. Sonra büyük sayılar yasasını ve merkezi limit teoremini simülasyonla gözleyeceğiz; bu iki teorem Monte Carlo yönteminin neden çalıştığını ve hatasının neden \(1/\sqrt{N}\) hızıyla küçüldüğünü açıklar. Bölüm Monte Carlo integrali ve rastgele yürüyüşle bitiyor.
14.1 Sözde Rastgele Sayılar
Bir bilgisayar programı aynı girdiyle her zaman aynı işi yapar; bu yüzden gerçekten rastgele bir sayı üretemez. Kullandığımız “rastgele” sayılar, rastgele görünen ama tamamen belirlenmiş bir diziden gelir.
Tanım 14.1 (Sözde Rastgele Sayı Üreteci) Bir sözde rastgele sayı üreteci (pseudorandom number generator) sonlu bir \(S\) durum kümesinden, bir \(T : S \to S\) geçiş fonksiyonundan ve bir \(G : S \to [0, 1)\) çıktı fonksiyonundan oluşur. Başlangıç durumu \(s_0 \in S\) tohum (seed) ile belirlenir ve üreteç
\[ s_k = T(s_{k-1}), \qquad u_k = G(s_k), \qquad k = 1, 2, \ldots \]
sayılarını üretir.
Yani bir üreteç, deterministik bir dizidir: aynı tohum her seferinde aynı \(u_1, u_2, \ldots\) sayılarını verir. İyi bir üretecin sayıları, bağımsız ve \(U(0, 1)\) dağılımlı (bkz. Olasılık Teorisi) sayılardan istatistiksel testlerle ayırt edilemez. \(S\) sonlu olduğundan dizi bir yerden sonra kendini tekrarlar; tekrarlanan parçanın uzunluğuna üretecin periyodu denir.
Örnek 14.1 (Doğrusal Kongrüans Üreteci) En eski üreteçlerden biri olan doğrusal kongrüans üretecinde (linear congruential generator) durumlar \(\{0, 1, \ldots, m - 1\}\) kümesindedir, geçiş \(x_{k+1} = (a x_k + c) \bmod m\) kuralıyla yapılır ve çıktı \(u_k = x_k / m\)’dir. \(m = 16\), \(a = 5\), \(c = 3\) ve \(x_0 = 7\) için dizinin ilk 20 terimini üretin ve periyodunu bulun.
Çözüm
Kuralı bir döngüyle uyguluyoruz. Periyodu bulmak için ilk terimin listede ikinci kez göründüğü yeri arıyoruz. Karşılaştırma için aynı üreteci bir de \(c = 2\) ile çalıştırıyoruz.
def lcg(seed, a=5, c=3, m=16, n=20):
"""x_{k+1} = (a*x_k + c) mod m dizisinin ilk n terimi."""
x = seed
xs = []
for _ in range(n):
x = (a * x + c) % m
xs.append(x)
return xs
xs = lcg(7)
print(xs)
print([x / 16 for x in xs[:5]]) # [0, 1) aralığına
print("periyot:", xs[1:].index(xs[0]) + 1)
ys = lcg(7, c=2) # c çift: kötü seçim
print(ys)
print("periyot:", ys[1:].index(ys[0]) + 1)Çıktı:
[6, 1, 8, 11, 10, 5, 12, 15, 14, 9, 0, 3, 2, 13, 4, 7, 6, 1, 8, 11]
[0.375, 0.0625, 0.5, 0.6875, 0.625]
periyot: 16
[5, 11, 9, 15, 13, 3, 1, 7, 5, 11, 9, 15, 13, 3, 1, 7, 5, 11, 9, 15]
periyot: 8İlk 16 terim \(0, 1, \ldots, 15\) sayılarının her birini tam bir kez içerir ve 17. terimden itibaren dizi baştan tekrarlanır: periyot \(16\)’dır, yani \(m\) kadar olabilecek en uzun değerdir. Bu bir rastlantı değildir. Sayılar teorisinde bilinen bir sonuca (Hull–Dobell teoremi) göre \(c \ne 0\) iken periyot ancak ve ancak şu üç koşul sağlanınca \(m\) olur: \(c\) ile \(m\) aralarında asaldır, \(a - 1\) sayısı \(m\)’nin her asal bölenine bölünür ve \(4 \mid m\) ise \(a - 1\) sayısı \(4\)’e de bölünür. Burada \(\gcd(3, 16) = 1\), \(m\)’nin tek asal böleni \(2\)’dir ve \(a - 1 = 4\) olduğundan üç koşul da sağlanır. \(c = 2\) seçince ilk koşul bozulur ve dizi yalnız 8 sayı arasında döner. \(\blacksquare\)
Böyle küçük bir üretecin sayıları rastgele görünmez; gerçek üreteçler aynı fikri çok daha büyük durum kümeleri ve daha karmaşık \(T\) ve \(G\) fonksiyonlarıyla uygular. NumPy’ın varsayılan üreteci PCG64’tür: 128 bitlik bir durumla çalışır ve periyodu \(2^{128} \approx 3{,}4 \cdot 10^{38}\)’dir. Saniyede bir milyar sayı üretseniz bile bu periyodu evrenin yaşı boyunca tüketemezsiniz.
Tanım 14.2 (Generator Nesnesi) np.random.default_rng(seed) çağrısı, PCG64 üretecini seed tamsayısını tohum alarak başlatır ve bir Generator nesnesi döndürür. Bu nesnenin random, integers, normal gibi metotları üretecin durumunu ilerleterek istenen dağılımdan sayılar üretir; metotların size argümanı çıktının şeklidir. seed verilmezse tohum işletim sisteminden alınan öngörülemez bir sayıdır.
Yani bütün rastgelelik tek bir nesnede, genellikle rng adını verdiğimiz üreteçte toplanır. Her çağrı yeni sayılar verir, çünkü durum ilerler; aynı sayıları yeniden görmek için üreteci aynı tohumla yeniden kurmak yeterlidir:
import numpy as np
rng = np.random.default_rng(42)
print(rng.random(3))
print(rng.random(3)) # durum ilerledi: yeni sayılar
rng = np.random.default_rng(42) # aynı tohumla baştan
print(rng.random(3))
print(type(rng).__name__, type(rng.bit_generator).__name__)Çıktı:
[0.77395605 0.43887844 0.85859792]
[0.69736803 0.09417735 0.97562235]
[0.77395605 0.43887844 0.85859792]
Generator PCG64İkinci rng.random(3) çağrısı birincinin devamıdır ve başka sayılar verir. Üreteci 42 tohumuyla yeniden kurunca ilk üç sayı aynen geri gelir. Bu notlardaki her kod sabit bir tohum kullanır; böylece kodu çalıştıran herkes aynı çıktıyı alır ve bir sonucu tekrarlamak mümkün olur.
Üreteci bir döngünün içinde aynı tohumla yeniden kurmak, her turda aynı sayıları üretir:
import numpy as np
for trial in range(3):
rng = np.random.default_rng(0) # YANLIŞ: her turda aynı tohum
print(rng.integers(1, 7, size=5))
rng = np.random.default_rng(0) # DOĞRU: üreteç bir kez kurulur
for trial in range(3):
print(rng.integers(1, 7, size=5))Çıktı:
[6 4 4 2 2]
[6 4 4 2 2]
[6 4 4 2 2]
[6 4 4 2 2]
[1 1 1 2 5]
[4 6 4 4 6]Üreteç programın başında bir kez kurulur ve sonra hep o kullanılır. Eski kodlarda sık görülen np.random.seed(0) ve np.random.rand(5) yazılışı ise bütün programın paylaştığı gizli, ortak bir üreteci kullanır. NumPy bu arayüzü yalnız eski kodlar çalışmaya devam etsin diye tutar; yeni kodda default_rng kullanılır.
14.2 Temel Dağılımlardan Örneklem
Generator nesnesinin metotları, Olasılık Teorisi notlarında tanıdığımız dağılımların hepsinden örneklem üretir. En sık kullanılanlar şunlardır:
| Metot | Dağılım |
|---|---|
rng.random(size) |
\(U(0, 1)\), değerler \([0, 1)\) aralığında |
rng.uniform(a, b, size) |
\(U(a, b)\) |
rng.integers(lo, hi, size) |
\(\{lo, lo + 1, \ldots, hi - 1\}\) üzerinde kesikli düzgün |
rng.normal(mu, sigma, size) |
\(N(\mu, \sigma^2)\) |
rng.exponential(scale, size) |
\(\operatorname{Üstel}(\lambda)\), scale \(= 1/\lambda\) |
rng.binomial(n, p, size) |
\(B(n, p)\) |
rng.poisson(lam, size) |
\(\operatorname{Poisson}(\lambda)\) |
rng.geometric(p, size) |
\(\operatorname{Geo}(p)\) |
rng.choice(a, size, p=p) |
a’nın elemanları arasından p olasılıklarıyla seçim |
rng.permutation(x) |
x’in elemanlarının rastgele bir sıralaması |
import numpy as np
rng = np.random.default_rng(7)
print(rng.uniform(-1, 1, size=3))
print(rng.integers(1, 7, size=10)) # 7 dahil değil
print(rng.normal(10, 2, size=3)) # ortalama 10, std 2
print(rng.exponential(0.5, size=3)) # ölçek 1/lambda = 0.5
print(rng.binomial(10, 0.3, size=10))
print(rng.poisson(4, size=10))
print(rng.choice(["yazı", "tura"], size=5))
print(rng.permutation(8))
print(rng.random((2, 3))) # size bir şekil de olabilirÇıktı:
[0.25019093 0.7944276 0.55137138]
[6 2 1 2 2 6 6 1 3 5]
[ 9.01558696 8.7590502 10.9796841 ]
[0.15607281 0.44988508 0.53685033]
[3 3 7 4 3 6 2 2 3 1]
[2 5 3 0 4 5 3 3 6 3]
['yazı' 'tura' 'tura' 'yazı' 'yazı']
[7 6 0 3 5 1 2 4]
[[0.84507432 0.94494817 0.90391679]
[0.56971915 0.14545995 0.19246349]]rng.choice bir listeden eleman seçer; p verilmezse bütün elemanlar eşit olasılıklıdır, replace=False ile aynı eleman ikinci kez seçilmez. rng.permutation(8) ise \(0, 1, \ldots, 7\) sayılarını karıştırır. Bir diziyi yerinde karıştırmak için rng.shuffle(x) de vardır.
Üç parametre sık karışır. rng.integers(lo, hi) üst sınırı dahil etmez: zar için integers(1, 7) yazılır (ya da integers(1, 6, endpoint=True)). rng.normal(mu, sigma) ikinci parametre olarak varyansı değil standart sapmayı alır: \(N(1, 4)\) için normal(1, 2) yazılır. rng.exponential(scale) oran parametresi \(\lambda\)’yı değil ortalama \(1/\lambda\)’yı alır: \(\operatorname{Üstel}(2)\) için exponential(0.5) yazılır.
Bir örneklemin dağılımına uyup uymadığını görmenin ilk yolu, örneklem ortalamasını ve örneklem varyansını (bkz. Matematiksel İstatistik) kuramsal değerlerle karşılaştırmaktır. x.mean() ortalamayı, x.var(ddof=1) ise \(n - 1\)’e bölen örneklem varyansını verir; ddof yazılmazsa NumPy \(n\)’ye böler.
import numpy as np
rng = np.random.default_rng(2025)
n = 10**6
samples = { # örneklem, E(X), Var(X)
"B(10, 0.3)": (rng.binomial(10, 0.3, n), 3.0, 2.1),
"Poisson(4)": (rng.poisson(4, n), 4.0, 4.0),
"Üstel(2)": (rng.exponential(0.5, n), 0.5, 0.25),
"N(1, 4)": (rng.normal(1, 2, n), 1.0, 4.0),
"U(0, 3)": (rng.uniform(0, 3, n), 1.5, 0.75),
}
print(f"{'dağılım':11}{'ort.':>9}{'E(X)':>7}{'varyans':>10}{'Var':>7}")
for name, (x, mean, var) in samples.items():
print(f"{name:11}{x.mean():9.4f}{mean:7.2f}"
f"{x.var(ddof=1):10.4f}{var:7.2f}")Çıktı:
dağılım ort. E(X) varyans Var
B(10, 0.3) 3.0014 3.00 2.1064 2.10
Poisson(4) 3.9993 4.00 4.0023 4.00
Üstel(2) 0.4995 0.50 0.2493 0.25
N(1, 4) 1.0002 1.00 3.9952 4.00
U(0, 3) 1.5013 1.50 0.7509 0.75Kuramsal değerler Olasılık Teorisi’nin moment formüllerinden gelir: binom için \(np\) ve \(np(1-p)\) (teorem), Poisson için \(\lambda\) ve \(\lambda\) (teorem), üstel dağılım için \(1/\lambda\) ve \(1/\lambda^2\) (teorem), normal dağılım için \(\mu\) ve \(\sigma^2\) (teorem), \(U(a, b)\) için \((a + b)/2\) ve \((b - a)^2/12\) (teorem). Bir milyon elemanlı örneklemlerde farklar binde birkaç mertebesindedir. Farkın neden tam bu büyüklükte olduğunu merkezi limit teoremi açıklayacak.
14.3 Histogram ve Ters Dönüşüm Yöntemi
Ortalama ve varyans, bir dağılımın yalnız iki sayısal özetidir. Dağılımın biçimini görmek için örneklemin histogramına bakarız.
Yoğunluk ölçeğinde histogramı Matplotlib ile Grafik Çizimi bölümünde tanımlamıştık: \(x_1, \ldots, x_n\) örnekleminde \([a_{j-1}, a_j)\) kutusuna \(c_j\) eleman düşüyorsa histogram o kutuda
\[ h(x) = \frac{c_j}{n\,(a_j - a_{j-1})}, \qquad a_{j-1} \le x < a_j \]
değerini alan basamak fonksiyonudur.
Böylece her çubuğun alanı, o kutuya düşen elemanların oranıdır; bütün elemanlar kutulara düşüyorsa çubukların toplam alanı \(1\)’dir. Örneklem \(f\) yoğunluklu bir dağılımdan geliyorsa \(c_j / n\) oranı \(P(a_{j-1} \le X < a_j) = \int_{a_{j-1}}^{a_j} f(x)\,dx\) olasılığına yakındır. Bu yüzden yoğunluk histogramı, büyük \(n\) ve dar kutular için yoğunluk fonksiyonunun grafiğine benzer. NumPy’da sayımı np.histogram yapar:
import math
import numpy as np
rng = np.random.default_rng(3)
n = 10_000
x = rng.normal(0, 1, size=n)
counts, edges = np.histogram(x, bins=8, range=(-4, 4))
print(counts, counts.sum())
print(edges)
width = np.diff(edges) # kutu genişlikleri
dens = counts / (n * width) # yoğunluk histogramı
print(dens.round(4))
print((dens * width).sum()) # çubukların toplam alanı
# kuramsal kutu olasılıkları: Phi(a_j) - Phi(a_{j-1})
Phi = [0.5 * (1 + math.erf(a / math.sqrt(2))) for a in edges]
print(np.diff(Phi).round(4))Çıktı:
[ 8 216 1361 3401 3433 1364 206 10] 9999
[-4. -3. -2. -1. 0. 1. 2. 3. 4.]
[0.0008 0.0216 0.1361 0.3401 0.3433 0.1364 0.0206 0.001 ]
0.9999
[0.0013 0.0214 0.1359 0.3413 0.3413 0.1359 0.0214 0.0013]np.histogram(x, bins=8, range=(-4, 4)) aralığı 8 eşit kutuya böler ve iki dizi döndürür: kutulardaki eleman sayıları ve 9 kutu sınırı. Son kutu dışındaki kutular sağdan açıktır; son kutu \([3, 4]\) kapalıdır. Aralığın dışında kalan değerler sayılmaz: on bin elemandan biri \([-4, 4]\) dışına düştüğü için sayıların toplamı \(9999\), çubukların toplam alanı da \(0{,}9999\)’dur. Son satırda kuramsal kutu olasılıkları \(\Phi(a_j) - \Phi(a_{j-1})\) var; \(\Phi\) standart normal dağılım fonksiyonudur ve math.erf hata fonksiyonuyla \(\Phi(z) = \frac{1}{2}\big(1 + \operatorname{erf}(z/\sqrt{2})\big)\) biçiminde hesaplanır. Histogramın değerleri bu olasılıklarla iki üç basamak uyuşuyor. Matplotlib’de aynı histogramı ax.hist(x, bins=8, range=(-4, 4), density=True) çizer (bkz. Matplotlib ile Grafik Çizimi); yalnız density=True, alanı aralığın içine düşen eleman sayısına göre \(1\) yapar.
Düzgün dağılımlı sayılar üretebiliyorsak başka her dağılımdan da örneklem üretebiliriz. Bunu sağlayan olgu şudur.
Önerme 14.1 (Ters Dönüşüm Yöntemi) \(F\) bir dağılım fonksiyonu, \(0 < u < 1\) için \(F^{-1}(u) = \inf\{x \in \mathbb{R} : F(x) \ge u\}\) ve \(U \sim U(0, 1)\) olsun. O zaman \(X = F^{-1}(U)\) rastgele değişkeninin dağılım fonksiyonu \(F\)’dir.
İspatı için bkz. Olasılık Teorisi. \(F\) sürekli ve kesin artan olduğunda \(F^{-1}\) sıradan ters fonksiyondur ve \(u = F(x)\) denklemini \(x\) için çözerek bulunur. Önerme, rng.random ile ürettiğimiz sayıları istediğimiz dağılıma taşımanın bir yolunu verir.
- İstenen dağılımın \(F\) dağılım fonksiyonunu yazın.
- \(u = F(x)\) denklemini \(x\) için çözerek \(F^{-1}(u)\)’yu bulun.
u = rng.random(n)ile düzgün sayılar üretipx = Finv(u)hesaplayın.
Örnek 14.2 (Üstel Dağılımdan Ters Dönüşümle Örneklem) \(\lambda = 2\) için \(\operatorname{Üstel}(\lambda)\) dağılımından ters dönüşüm yöntemiyle \(10^5\) elemanlı bir örneklem üretin. Örneklemin ortalamasını, varyansını ve \(P(X > 1)\) için verdiği kestirimi kuramsal değerlerle karşılaştırın.
Aşağıdaki şekil yöntemin işleyişini ve ürettiği örneklemin histogramını gösteriyor.
default_rng(11) üretecinin ilk altı u sayısı dikey eksenden F(x) = 1 − e−2x eğrisine, oradan yatay eksendeki x = −ln(1 − u)/2 noktasına taşınır. Eğri sıfırın yakınında dik olduğundan u ekseninin geniş bir parçası küçük x'lerin dar bir aralığına taşınır; x'ler bu yüzden sıfırın yakınında sıklaşır. Sağda 100 000 değerin 30 kutulu yoğunluk histogramı ve Üstel(2) yoğunluğu 2e−2x.Çözüm
Üç adımı uyguluyoruz.
- \(\operatorname{Üstel}(\lambda)\) dağılımının yoğunluğu \(x > 0\) için \(\lambda e^{-\lambda x}\)’tir (bkz. Olasılık Teorisi); integralini alarak \(x \ge 0\) için \(F(x) = 1 - e^{-\lambda x}\) buluruz.
- \(u = 1 - e^{-\lambda x}\) denkleminden \(e^{-\lambda x} = 1 - u\), yani \(F^{-1}(u) = -\ln(1 - u)/\lambda\) olur.
- Kodda
udizisine bu formülü uyguluyoruz.
import numpy as np
lam = 2.0
n = 100_000
rng = np.random.default_rng(11)
u = rng.random(n)
x = -np.log(1 - u) / lam # x = F^{-1}(u)
print(u[:4].round(4))
print(x[:4].round(4))
print(x.mean(), x.var(ddof=1)) # teori: 0.5 ve 0.25
print((x > 1).mean(), np.exp(-lam)) # P(X > 1) = e^(-2)
counts, edges = np.histogram(x, bins=30, range=(0, 3))
dens = counts / (n * np.diff(edges))
mid = (edges[:-1] + edges[1:]) / 2
print(np.abs(dens - lam * np.exp(-lam * mid)).max())Çıktı:
[0.1286 0.4993 0.6015 0.0287]
[0.0688 0.3459 0.46 0.0146]
0.5001468651191905 0.2495862780343443
0.1359 0.1353352832366127
0.011336413931975042Kuramsal değerler \(E(X) = 1/\lambda = 0{,}5\), \(\operatorname{Var}(X) = 1/\lambda^2 = 0{,}25\) ve \(P(X > 1) = e^{-2} \approx 0{,}1353\)’tür; örneklem bunları \(0{,}5001\), \(0{,}2496\) ve \(0{,}1359\) olarak veriyor. Son satır, 30 kutulu yoğunluk histogramının yüksekliklerinin kutu ortalarındaki \(2e^{-2x}\) değerlerinden en çok \(0{,}011\) saptığını gösteriyor. \(1 - U\) de \(U(0, 1)\) dağılımlı olduğundan -np.log(u) / lam da aynı dağılımı verir. NumPy’ın rng.exponential metodu ise daha hızlı olan başka bir yöntem kullanır. \(\blacksquare\)
14.4 Olasılıkları ve Beklenen Değerleri Kestirmek
Giriş örneğinde bir olasılığı, deneyi tekrarlayıp göreli sıklığa bakarak kestirmiştik. Aynı fikir her beklenen değer için çalışır.
Tanım 14.3 (Monte Carlo Kestiricisi) \(X\) bir rastgele değişken ya da rastgele vektör, \(g\) bir fonksiyon ve \(\theta = E[g(X)]\) sonlu olsun. \(X_1, \ldots, X_N\), \(X\) ile aynı dağılımlı ve bağımsız ise
\[ \hat\theta_N = \frac{1}{N} \sum_{i=1}^{N} g(X_i) \]
sayısına \(\theta\)’nın \(N\) örneklemli Monte Carlo kestiricisi denir.
Yani beklenen değeri, deneyi \(N\) kez yapıp \(g\)’nin aldığı değerlerin ortalamasıyla kestiriyoruz. Bir \(A\) olayı için \(g\), \(A\) gerçekleşince \(1\), gerçekleşmeyince \(0\) değerini alan gösterge fonksiyonu olursa \(E[g(X)] = P(A)\) olur ve \(\hat\theta_N\), \(A\)’nın \(N\) denemedeki göreli sıklığıdır. NumPy’da \(X_i\)’ler bir dizinin satırları, \(g\) vektörel bir ifade, \(\hat\theta_N\) de bir mean çağrısıdır.
- \(N\) denemenin hepsini tek bir diziyle üretin: her satır bir deneme olsun, örneğin
size=(N, k). - Olayı her satır için bir boolean değere çevirin; karşılaştırmalar ve
axis=1ileany,all,sumbunun araçlarıdır. - Bu boolean dizinin
meandeğerini alın.
Örnek 14.3 (Doğum Günü Problemi) Bir yıl 365 gün olsun ve herkesin doğum günü bu günlerden birine eşit olasılıkla, birbirinden bağımsız olarak düşsün. 23 kişilik bir grupta en az iki kişinin doğum gününün aynı olma olasılığını \(10^5\) denemeyle kestirin ve tam değerle karşılaştırın.
Çözüm
Üç adımı uyguluyoruz. Her deneme 23 kişinin doğum günleridir, yani (N, 23) şekilli bir dizinin bir satırı. Bir satırda aynı gün iki kez geçiyor mu sorusunu yanıtlamak için satırı sıralarız: sıralı bir listede eşit elemanlar yan yana gelir, dolayısıyla ardışık farklardan biri \(0\) olur.
import numpy as np
rng = np.random.default_rng(23)
N, k = 100_000, 23
days = rng.integers(0, 365, size=(N, k)) # her satır bir grup
days.sort(axis=1) # satırları sırala
shared = (np.diff(days, axis=1) == 0).any(axis=1)
print(shared[:6])
print(shared.mean())
exact = 1 - np.prod((365 - np.arange(k)) / 365)
print(exact)Çıktı:
[ True True False True True True]
0.50956
0.5072972343239857days.sort(axis=1) her satırı yerinde sıralar, np.diff(days, axis=1) satır içindeki ardışık farkları alır, .any(axis=1) her satırda sıfır fark olup olmadığını söyler. Simülasyonun kestirimi \(0{,}5096\)’dır. Tam değer, hiç çakışma olmaması olasılığının tümleyenidir:
\[ P(\text{çakışma}) = 1 - \frac{365 \cdot 364 \cdots 343}{365^{23}} = 1 - \prod_{k=0}^{22} \frac{365 - k}{365} \approx 0{,}5073. \]
Bu formülün sayma ile çıkarılışı için bkz. Olasılık Teorisi. Kodun son satırı çarpımı np.prod ile hesaplıyor. Kestirim ile tam değer arasındaki fark \(0{,}002\) civarındadır; bu büyüklükteki farkların neden beklendiğini merkezi limit teoremi gösterecek. \(\blacksquare\)
14.5 Büyük Sayılar Yasası
Monte Carlo kestiricisinin doğru değere gittiğini garanti eden teorem büyük sayılar yasasıdır. Önce bir zarı on bin kez atıp ortalamanın nasıl değiştiğini izleyelim. Bir dizinin ilk \(n\) elemanının ortalamasına, \(n\) ilerledikçe güncellendiği için koşu ortalaması (running mean) denir.
import numpy as np
rng = np.random.default_rng(5)
n = 10_000
rolls = rng.integers(1, 7, size=(5, n)) # beş ayrı deney
k = np.arange(1, n + 1)
means = rolls.cumsum(axis=1) / k # satır satır X̄_1, ..., X̄_n
for m in [10, 100, 1000, 10_000]:
print(m, means[:, m - 1].round(3))Çıktı:
10 [3.6 3.3 3.2 3.3 2.9]
100 [3.3 3.59 3.27 3.7 3.47]
1000 [3.451 3.422 3.41 3.477 3.523]
10000 [3.5 3.49 3.497 3.528 3.512]rolls.cumsum(axis=1) her satırın kısmi toplamlarını, yani \(X_1 + \cdots + X_n\) değerlerini verir. Bunu broadcasting ile k = [1, 2, ..., n] dizisine bölünce bütün koşu ortalamaları tek ifadeyle hesaplanır. On atışta ortalamalar \(2{,}9\) ile \(3{,}6\) arasında dağınıktır; on bin atışta hepsi \(3{,}5\)’e yüzde birden daha yakındır. Bir zarın beklenen değeri \((1 + 2 + \cdots + 6)/6 = 3{,}5\)’tir.
Teorem 14.1 (Güçlü Büyük Sayılar Yasası) \(X_1, X_2, \ldots\) bağımsız, aynı dağılımlı ve \(E|X_1| < \infty\) olan rastgele değişkenler, \(\mu = E(X_1)\) olsun. O zaman olasılığı \(1\) olan bir olay üzerinde
\[ \bar X_n = \frac{X_1 + X_2 + \cdots + X_n}{n} \longrightarrow \mu \qquad (n \to \infty) \]
olur.
Aşağıdaki şeklin solunda yukarıdaki kodun beş zar deneyi, sağında ise teoremin hipotezini sağlamayan bir dağılımın beş yolu var.
Yani bir deneyi yeterince çok tekrarlarsak sonuçların ortalaması neredeyse kesinlikle beklenen değere yaklaşır. Monte Carlo kestiricisi \(\hat\theta_N\), \(g(X_i)\) değerlerinin ortalaması olduğundan \(E|g(X)| < \infty\) ise \(\hat\theta_N \to \theta\) olur. Teoremin zayıf biçimi ve ispatı için bkz. Matematiksel İstatistik, güçlü biçimi için aynı bölümdeki teorem.
Teoremdeki \(E|X_1| < \infty\) koşulu bir ayrıntı değildir. Yoğunluğu \(\frac{1}{\pi(1 + x^2)}\) olan standart Cauchy dağılımında \(\int |x| / (\pi(1 + x^2))\,dx\) integrali ıraksar, yani beklenen değer yoktur:
import numpy as np
rng = np.random.default_rng(6)
n = 10_000
x = rng.standard_cauchy(size=(5, n))
means = x.cumsum(axis=1) / np.arange(1, n + 1)
for m in [10, 100, 1000, 10_000]:
print(m, means[:, m - 1].round(2))
print(np.abs(x).max().round(1))Çıktı:
10 [ 1.87 -0.9 -1.58 -0.85 0.87]
100 [ 2.96 6.21 -1.11 -0.34 -0.69]
1000 [-0.48 3.14 0.08 -0.72 13.03]
10000 [ -1.41 0.62 1.78 -14.78 5.75]
154047.4Koşu ortalamaları hiçbir sayıya yerleşmez. Arada bir gelen çok büyük bir değer (beş yoldaki elli bin sayının mutlak değerce en büyüğü yüz elli binden büyük) ortalamayı birden sıçratır. Hatta \(n\) tane bağımsız standart Cauchy değişkenin ortalaması yine standart Cauchy dağılımlıdır: ortalama almak hiçbir şeyi iyileştirmez. Bir Monte Carlo hesabında ortalamalar yerleşmiyorsa, kestirilen beklenen değerin sonlu olup olmadığını kontrol edin.
14.6 Merkezi Limit Teoremi
Büyük sayılar yasası ortalamanın \(\mu\)’ye gittiğini söyler ama ne kadar hızlı gittiğini söylemez. Bu soruyu merkezi limit teoremi yanıtlar: \(\bar X_n - \mu\) farkı \(\sigma/\sqrt{n}\) mertebesindedir ve yaklaşık normal dağılımlıdır. Bunu çok çarpık bir dağılımla, \(\operatorname{Üstel}(1)\) ile deneyelim; bu dağılımda \(\mu = \sigma = 1\)’dir. Her \(n\) için \(10^5\) kez \(n\) elemanlı bir örneklem alıp ortalamasını \(Z = (\bar X_n - \mu)/(\sigma/\sqrt{n})\) biçiminde standartlaştırıyor ve \(Z\)’lerin histogramını standart normal yoğunlukla karşılaştırıyoruz:
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(8)
reps = 100_000 # her n için tekrar sayısı
t = np.linspace(-4, 4, 201)
phi = np.exp(-t**2 / 2) / np.sqrt(2 * np.pi)
fig, axs = plt.subplots(1, 3, figsize=(12, 3.5), sharey=True,
layout="constrained")
for ax, n in zip(axs, [1, 4, 32]):
x = rng.exponential(1.0, size=(reps, n)) # mu = sigma = 1
z = (x.mean(axis=1) - 1.0) * np.sqrt(n) # standartlaştır
ax.hist(z, bins=40, range=(-4, 4), density=True, alpha=0.6)
ax.plot(t, phi, "k")
ax.set_title(f"n = {n}")
p1, p0 = np.mean(z <= -1), np.mean(z <= 0)
print(f"n = {n:2d}: P(Z <= -1) ~ {p1:.4f}, P(Z <= 0) ~ {p0:.4f}")
fig.savefig("merkezi_limit.png", dpi=150)Çıktı:
n = 1: P(Z <= -1) ~ 0.0000, P(Z <= 0) ~ 0.6308
n = 4: P(Z <= -1) ~ 0.1427, P(Z <= 0) ~ 0.5657
n = 32: P(Z <= -1) ~ 0.1558, P(Z <= 0) ~ 0.5205rng.exponential(1.0, size=(reps, n)) her satırı bir örneklem olan bir dizi üretir ve x.mean(axis=1) satır ortalamalarını alır. Standart normal dağılımda \(\Phi(-1) \approx 0{,}1587\) ve \(\Phi(0) = 0{,}5\)’tir. \(n = 1\)’de \(Z = X - 1 \ge -1\) olduğundan \(P(Z \le -1) = 0\)’dır ve \(P(Z \le 0) = 1 - e^{-1} \approx 0{,}632\)’dir. \(n\) büyüdükçe iki olasılık da normal değerlerine yaklaşıyor.
Teorem 14.2 (Merkezi Limit Teoremi) \(X_1, X_2, \ldots\) bağımsız, aynı dağılımlı, beklenen değeri \(\mu\) ve varyansı \(0 < \sigma^2 < \infty\) olan rastgele değişkenler ve \(\Phi\) standart normal dağılım fonksiyonu olsun. O zaman her \(z \in \mathbb{R}\) için
\[ \lim_{n \to \infty} P\left( \frac{\bar X_n - \mu}{\sigma / \sqrt{n}} \le z \right) = \Phi(z) \]
olur.
Aşağıdaki şekil yukarıdaki kodun çizdiği üç histogramdır.
Yani ortalamanın \(\mu\)’den sapması, \(n\) büyüdükçe \(\sigma/\sqrt{n}\) ölçekli bir normal dağılıma benzer; başlangıçtaki dağılım ne kadar çarpık olursa olsun. Özel olarak
\[ P\left( |\bar X_n - \mu| \le 1{,}96\, \frac{\sigma}{\sqrt{n}} \right) \approx \Phi(1{,}96) - \Phi(-1{,}96) \approx 0{,}95 \]
olur. Moment üreten fonksiyonun var olduğu durumda ispat için bkz. Matematiksel İstatistik. Önceki bölümlerdeki farklar bu yüzden beklenen büyüklükteydi. Örneğin doğum günü probleminde bir denemenin sonucu olasılığı \(p \approx 0{,}507\) olan bir Bernoulli değişkenidir, yani \(\sigma = \sqrt{p(1 - p)} \approx 0{,}5\)’tir. \(N = 10^5\) için \(\sigma/\sqrt{N} \approx 0{,}0016\) olur ve bulduğumuz \(0{,}002\)’lik fark bu ölçekte sıradandır.
Örnek 14.4 (On İki Düzgün Sayıdan Yaklaşık Normal Sayı) \(U_1, \ldots, U_{12}\) bağımsız ve \(U(0, 1)\) dağılımlı olmak üzere \(Z = U_1 + \cdots + U_{12} - 6\) olsun. \(10^6\) tane \(Z\) üretin. \(z = 0{,}5\), \(1\) ve \(2\) için \(P(Z \le z)\) kestirimlerini \(\Phi(z)\) ile, \(P(|Z| > 3)\) kestirimini de standart normal dağılımın değeriyle karşılaştırın.
Çözüm
Önce \(Z\)’nin ortalama ve varyansını bulalım. \(E(U_i) = 1/2\) ve \(\operatorname{Var}(U_i) = 1/12\) olduğundan beklenen değerin lineerliği ile \(E(Z) = 12 \cdot \frac{1}{2} - 6 = 0\), bağımsız değişkenlerin varyansları toplandığı için de \(\operatorname{Var}(Z) = 12 \cdot \frac{1}{12} = 1\) olur. Merkezi limit teoremine göre \(Z\) yaklaşık \(N(0, 1)\) dağılımlıdır. \(\Phi\)’yi math.erf ile hesaplıyoruz.
import math
import numpy as np
def Phi(z):
"""Standart normal dağılım fonksiyonu."""
return 0.5 * (1 + math.erf(z / math.sqrt(2)))
rng = np.random.default_rng(12)
N = 10**6
Z = rng.random((N, 12)).sum(axis=1) - 6 # her satır 12 düzgün sayı
print(Z.mean(), Z.var(ddof=1))
for z in [0.5, 1.0, 2.0]:
print(z, np.mean(Z <= z), Phi(z))
print(np.mean(np.abs(Z) > 3), 2 * (1 - Phi(3)))
print(np.abs(Z).max())Çıktı:
0.0008496661665800563 1.002624404897434
0.5 0.688608 0.691462461274013
1.0 0.838893 0.8413447460685428
2.0 0.977393 0.9772498680518209
0.002056 0.002699796063260207
4.557258791575404Örneklemin ortalaması ve varyansı \(0\) ve \(1\)’e yakındır. \(P(Z \le z)\) kestirimleri \(\Phi(z)\) ile iki basamak uyuşuyor: dağılımın orta kısmı için on iki terim yeterlidir. Kuyrukta durum farklıdır: \(P(|Z| > 3)\) için kestirim \(0{,}0021\), normal değer \(0{,}0027\)’dir. Bu fark simülasyon hatası değildir; toplamın tam dağılımından hesaplanan değer de \(0{,}0020\) civarındadır. \(Z\) hiçbir zaman \([-6, 6]\) aralığının dışına çıkamaz, oysa normal dağılım her değeri alabilir; bir milyon \(Z\) arasında en büyük mutlak değer yalnız \(4{,}56\)’dır. Merkezi limit teoremi dağılımın merkezini iyi açıklar, uç olayların olasılıklarını ise ancak büyük \(n\) için. \(\blacksquare\)
14.7 Monte Carlo ile İntegral
Bir integral de bir beklenen değer olarak yazılabilir; bu gözlem Monte Carlo yöntemini integral hesabının bir aracına çevirir. \(U \sim U(a, b)\) ise \(U\)’nun yoğunluğu \([a, b]\) üzerinde \(1/(b - a)\) olduğundan
\[ E[f(U)] = \int_a^b f(x)\,\frac{dx}{b - a}, \qquad \text{yani} \qquad \int_a^b f(x)\,dx = (b - a)\,E[f(U)] \]
olur. Aynı akıl yürütme birden çok boyutta da geçerlidir.
Tanım 14.4 (Monte Carlo İntegrali) \(D \subset \mathbb{R}^d\) hacmi \(|D|\) olan sınırlı bir bölge ve \(f\), \(D\) üzerinde integrallenebilir bir fonksiyon olsun. \(U_1, \ldots, U_N\), \(D\) içinde bağımsız ve düzgün dağılımlı noktalar ise
\[ \hat I_N = \frac{|D|}{N} \sum_{i=1}^{N} f(U_i) \]
sayısına \(I = \int_D f\) integralinin \(N\) noktalı Monte Carlo kestirimi denir.
Yani integral, bölgenin hacmiyle \(f\)’nin \(D\) üzerindeki ortalama değerinin (bkz. İntegral Calculus) çarpımıdır ve ortalama değeri rastgele noktalardaki değerlerin ortalamasıyla kestiriyoruz. Bu, Monte Carlo kestiricisinin (Tanım 14.3) \(g = |D|\,f\) için özel hâlidir. İlk örnek olarak \(e^{-x^2}\)’nin \([0, 1]\) üzerindeki integralini alalım. Bu fonksiyonun ilkeli elemanter fonksiyonlarla yazılamaz; integralin değeri hata fonksiyonuyla \(\frac{\sqrt{\pi}}{2}\operatorname{erf}(1) \approx 0{,}746824\)’tür.
import math
import numpy as np
def f(x):
return np.exp(-x**2)
rng = np.random.default_rng(14)
N = 10_000
u = rng.random(N) # [0, 1) aralığında düzgün noktalar
y = f(u)
I_hat = y.mean() # (b - a) = 1 olduğu için
print(u[:4].round(4))
print(y[:4].round(4))
exact = math.sqrt(math.pi) / 2 * math.erf(1)
print(I_hat, exact, abs(I_hat - exact))Çıktı:
[0.831 0.3609 0.7027 0.8601]
[0.5013 0.8778 0.6103 0.4772]
0.7476505933702069 0.746824132812427 0.0008264605577799067\(b - a = 1\) olduğundan kestirim, \(f\)’nin rastgele noktalardaki değerlerinin ortalamasıdır. On bin noktayla bulunan \(0{,}74765\), gerçek değerden yaklaşık \(8 \cdot 10^{-4}\) uzaktadır. Aşağıdaki şekil bu hesabı geometrik olarak gösteriyor.
default_rng(14) ile üretilen ilk 30 u sayısında e−u² değerleridir. Kesikli dikdörtgenin yüksekliği, 10 000 değerin ortalaması olan 0,7477 sayısıdır; dikdörtgenin alanı, eğrinin altında kalan taralı alanın (0,7468) kestirimidir. Eğrinin dikdörtgenin üstüne taşan parçasıyla dikdörtgenin içinde eğrinin üstünde kalan boşluğun alanları yaklaşık eşittir.Kestirimin hatasını kesin olarak bilemeyiz, çünkü bilseydik integrali de bilirdik. Ama hatanın tipik büyüklüğünü kestirebiliriz.
Önerme 14.2 (Monte Carlo Kestiriminin Hatası) \(\hat I_N\), \(I = \int_D f\) integralinin bir Monte Carlo kestirimi (Tanım 14.4) ve \(0 < \sigma^2 = \operatorname{Var}\big(f(U_1)\big) < \infty\) olsun. O zaman
\[ E(\hat I_N) = I, \qquad \operatorname{Var}(\hat I_N) = \frac{|D|^2 \sigma^2}{N} \]
olur. Ayrıca \(N \to \infty\) iken
\[ P\left( |\hat I_N - I| \le 1{,}96\, \frac{|D|\,\sigma}{\sqrt{N}} \right) \longrightarrow \Phi(1{,}96) - \Phi(-1{,}96) \approx 0{,}95 \]
olur.
İspat
\(Y_i = |D|\,f(U_i)\) diyelim. \(U_i\)’ler bağımsız ve aynı dağılımlı olduğundan \(Y_i\)’ler de öyledir ve \(\hat I_N = \bar Y_N\) olur. \(U_i\)’nin yoğunluğu \(D\) üzerinde \(1/|D|\) olduğundan
\[ E(Y_i) = |D| \int_D f(x)\,\frac{dx}{|D|} = I, \qquad \operatorname{Var}(Y_i) = |D|^2 \sigma^2 \]
olur. Bağımsız ve aynı ortalamalı, aynı varyanslı gözlemlerin ortalaması için \(E(\bar Y_N) = E(Y_1)\) ve \(\operatorname{Var}(\bar Y_N) = \operatorname{Var}(Y_1)/N\)’dir (bkz. Matematiksel İstatistik); bu ilk iki eşitliği verir.
Son iddia için \(Z_N = \sqrt{N}\,(\bar Y_N - I)/(|D|\,\sigma)\) diyelim; iddiadaki olay \(|Z_N| \le 1{,}96\) olayıdır. \(Y_i\)’lerin beklenen değeri \(I\), standart sapması \(|D|\,\sigma\) olduğundan merkezi limit teoremi (Teorem 14.2) her \(z\) için \(P(Z_N \le z) \to \Phi(z)\) verir. \(\Phi\) sürekli olduğu için \(P(Z_N < z) \to \Phi(z)\) de geçerlidir. Böylece
\[ \begin{aligned} P(|Z_N| \le 1{,}96) &= P(Z_N \le 1{,}96) - P(Z_N < -1{,}96) \\[1mm] &\quad \longrightarrow \Phi(1{,}96) - \Phi(-1{,}96) \end{aligned} \]
olur. Bu limitin değeri \(0{,}950004\ldots\), yani yaklaşık \(0{,}95\)’tir. \(\blacksquare\)
Yani kestirimin standart hatası (bkz. Matematiksel İstatistik) \(|D|\,\sigma/\sqrt{N}\)’dir: hata \(N\) ile değil \(\sqrt{N}\) ile küçülür. Hatayı 10 kat küçültmek için 100 kat fazla nokta gerekir. \(\sigma\)’yı bilmediğimiz için onu, zaten hesapladığımız \(f(U_i)\) değerlerinin örneklem standart sapması \(s\) ile kestiririz. Böylece
\[ \hat I_N \pm 1{,}96\, \frac{|D|\, s}{\sqrt{N}} \]
aralığı, \(I\) için yaklaşık %95 düzeyli bir güven aralığı (bkz. Matematiksel İstatistik) olur. Önceki hesaba standart hatayı ekleyelim:
import math
import numpy as np
def f(x):
return np.exp(-x**2)
rng = np.random.default_rng(14)
N = 10_000
y = f(rng.random(N))
I_hat = y.mean()
s = y.std(ddof=1) # sigma'nın kestirimi
se = s / math.sqrt(N) # standart hata
print(s, se)
print(f"I = {I_hat:.4f} ± {1.96 * se:.4f}")
print(I_hat - 1.96 * se, I_hat + 1.96 * se)Çıktı:
0.2013487113344417 0.002013487113344417
I = 0.7477 ± 0.0039
0.7437041586280518 0.751597028112362Standart hata \(0{,}0020\)’dir ve %95 güven aralığı \([0{,}7437;\ 0{,}7516]\) gerçek değer \(0{,}7468\)’i içerir. Gerçek hata \(0{,}0008\), standart hatadan küçüktür. Güven aralığının ne demek olduğunu bir deneyle görelim: aynı hesabı \(N = 1000\) noktayla 20 kez bağımsız olarak yapıp her seferinde aralığın gerçek değeri içerip içermediğine bakalım.
import math
import numpy as np
def f(x):
return np.exp(-x**2)
exact = math.sqrt(math.pi) / 2 * math.erf(1)
rng = np.random.default_rng(16)
N, runs = 1000, 20
hits = 0
for k in range(runs):
y = f(rng.random(N))
I_hat = y.mean()
h = 1.96 * y.std(ddof=1) / math.sqrt(N) # yarı genişlik
ok = I_hat - h <= exact <= I_hat + h
hits += ok
if k < 3:
print(f"[{I_hat - h:.4f}, {I_hat + h:.4f}]", ok)
print(hits, "/", runs)
# 10 000 tekrar: tek bir (10000, N) dizisiyle
Y = f(rng.random((10_000, N)))
I_hats = Y.mean(axis=1)
hs = 1.96 * Y.std(axis=1, ddof=1) / math.sqrt(N)
print(np.mean(np.abs(I_hats - exact) <= hs))Çıktı:
[0.7291, 0.7544] True
[0.7303, 0.7551] True
[0.7324, 0.7573] True
19 / 20
0.9492Yirmi aralıktan 19’u gerçek değeri içeriyor. İkinci kısımda aynı deneyi tek bir (10000, 1000) şekilli diziyle on bin kez yapıyoruz; aralıkların \(\%94{,}9\)’u gerçek değeri yakalıyor. “%95 güven” sözü tek bir aralık hakkında değil, aralığı üreten yöntem hakkındadır.
- Bölgeyi hacmi bilinen bir kutunun içine alın; \(f\)’yi kutunun bölge dışında kalan kısmında \(0\) kabul edin. Böylece integral kutu üzerinde bir integral olur.
- Kutuda
rng.uniformile \(N\) tane düzgün dağılımlı nokta üretin. - Fonksiyon değerlerini
valsdizisine yazıp kestirimivol * vals.mean()ile hesaplayın. - Standart hatayı
vol * vals.std(ddof=1) / np.sqrt(N)ile bulun ve sonucu kestirim \(\pm\, 1{,}96 \cdot\) standart hata aralığıyla bildirin.
Örnek 14.5 (Yarım Halka Üzerinde İki Katlı İntegral) \(R\), üst yarı düzlemde \(x^2 + y^2 = 1\) ve \(x^2 + y^2 = 4\) çemberleri arasında kalan bölge olsun. \(\iint_R (3x + 4y^2)\,dA\) integralini \(10^6\) noktalı Monte Carlo yöntemiyle, %95 güven aralığıyla kestirin ve kutupsal koordinatlarla bulunan tam değerle karşılaştırın.
Çözüm
Dört adımı uyguluyoruz.
- \(R\), \([-2, 2] \times [0, 2]\) kutusunun içindedir; kutunun alanı \(8\)’dir. Kutu üzerinde \(R\) içinde \(3x + 4y^2\), dışında \(0\) olan fonksiyonun integralini alacağız.
- Kutuda \(10^6\) nokta üretiyoruz; \(x\) ve \(y\) koordinatları ayrı ayrı düzgün dağılımlıdır.
np.wherenoktanın \(R\)’de olup olmadığına göre fonksiyon değerini ya da \(0\)’ı seçer.- Standart hatayı ve güven aralığını hesaplıyoruz.
import numpy as np
rng = np.random.default_rng(18)
N = 10**6
x = rng.uniform(-2, 2, N) # kutu: [-2, 2] x [0, 2]
y = rng.uniform(0, 2, N)
r2 = x**2 + y**2
in_R = (1 <= r2) & (r2 <= 4) # nokta yarım halkada mı?
vals = np.where(in_R, 3 * x + 4 * y**2, 0.0)
vol = 4 * 2 # kutunun alanı
I_hat = vol * vals.mean()
se = vol * vals.std(ddof=1) / np.sqrt(N)
print(in_R.mean(), 3 * np.pi / 16) # R'nin kutudaki payı
print(I_hat, se)
print(I_hat - 1.96 * se, I_hat + 1.96 * se)
print(15 * np.pi / 2)Çıktı:
0.589537 0.5890486225480862
23.623947891599368 0.038170669799821234
23.54913337879172 23.698762404407017
23.561944901923447İlk satır bir kontrol: noktaların \(R\)’ye düşme oranı \(0{,}5895\)’tir ve \(R\)’nin alanının kutunun alanına oranı olan \(\frac{(4\pi - \pi)/2}{8} = \frac{3\pi}{16} \approx 0{,}5890\)’a yakındır. Kestirim \(23{,}62\), standart hata \(0{,}038\)’dir. Kutupsal koordinatlarla integral
\[ \int_0^{\pi}\!\! \int_1^2 \left(3r\cos\theta + 4r^2\sin^2\theta\right) r\,dr\,d\theta = 0 + 4 \cdot \frac{\pi}{2} \cdot \frac{15}{4} = \frac{15\pi}{2} \approx 23{,}562 \]
olarak bulunur (ayrıntılı hesap için bkz. İntegral Calculus). Tam değer \([23{,}549;\ 23{,}699]\) güven aralığının içindedir. \(\blacksquare\)
14.8 Monte Carlo Ne Zaman Kazanır?
\(1/\sqrt{N}\) hızı yavaş bir hızdır. Bir boyutlu integrallerde Nümerik Analiz notlarındaki yamuk kuralı çok daha iyi sonuç verir. Aşağıdaki kod iki yöntemi aynı sayıda fonksiyon değeriyle karşılaştırıyor. Monte Carlo hatası rastgele olduğu için her \(N\) değerinde hesabı 400 kez tekrarlayıp hataların karelerinin ortalamasının karekökünü alıyoruz. Bunu \(\sigma/\sqrt{N}\) ile karşılaştırıyoruz; \(\sigma^2 = E[f(U)^2] - I^2\) ve \(E[f(U)^2] = \int_0^1 e^{-2x^2}\,dx\) yine hata fonksiyonuyla hesaplanır. Yamuk kuralı için Sayısal Türev ve İntegral bölümünde gördüğümüz np.trapezoid fonksiyonunu kullanıyoruz.
import math
import numpy as np
def f(x):
return np.exp(-x**2)
exact = math.sqrt(math.pi) / 2 * math.erf(1)
Ef2 = math.sqrt(math.pi / 8) * math.erf(math.sqrt(2)) # E[f(U)^2]
sigma = math.sqrt(Ef2 - exact**2)
print("sigma =", sigma)
rng = np.random.default_rng(17)
reps = 400 # her N için tekrar sayısı
print(f"{'N':>7}{'MC hatası':>11}{'sigma/kök N':>13}{'yamuk':>10}")
for N in [10, 100, 1000, 10_000, 100_000]:
est = np.array([f(rng.random(N)).mean() for _ in range(reps)])
rms = np.sqrt(np.mean((est - exact) ** 2))
x = np.linspace(0, 1, N) # yamuk kuralı da N nokta kullanır
trap = np.trapezoid(f(x), x)
print(f"{N:7d}{rms:11.2e}{sigma / math.sqrt(N):13.2e}"
f"{abs(trap - exact):10.2e}")Çıktı:
sigma = 0.2009918438899211
N MC hatası sigma/kök N yamuk
10 6.20e-02 6.36e-02 7.57e-04
100 1.96e-02 2.01e-02 6.26e-06
1000 6.37e-03 6.36e-03 6.14e-08
10000 2.01e-03 2.01e-03 6.13e-10
100000 6.33e-04 6.36e-04 6.13e-12Monte Carlo sütunu \(\sigma/\sqrt{N}\) sütunuyla neredeyse aynıdır: Önerme 14.2 tam olarak gözleniyor. \(N\) her 10 katına çıktığında Monte Carlo hatası yalnız \(\sqrt{10} \approx 3{,}16\) kat küçülür, yamuk kuralının hatası ise 100 kat küçülür; çünkü yamuk kuralının hatası \(h^2\) ile, yani \(N^{-2}\) ile orantılıdır. Bir boyutta Monte Carlo açık ara kaybeder.
Tablo yüksek boyutta tersine döner. \(d\) boyutlu bir küpte her eksende \(m\) nokta kullanan bir ızgara \(N = m^d\) fonksiyon değeri gerektirir ve yamuk kuralının hatası \(h^2 \sim m^{-2} = N^{-2/d}\) ile küçülür. \(d = 4\)’te bu Monte Carlo’nun \(N^{-1/2}\) hızına eşittir, daha yüksek boyutta ondan yavaştır. Monte Carlo hatasının hızı ise boyuttan bağımsızdır: standart hata formülünde (Önerme 14.2) \(d\) hiç geçmez. Izgara yöntemlerinin boyut arttıkça çaresiz kalmasına boyutun laneti (curse of dimensionality) denir.
Örnek 14.6 (On Boyutlu Birim Topun Hacmi) \(\mathbb{R}^{10}\)’daki birim topun hacmi \(\pi^5/120 \approx 2{,}5502\)’dir. Bu hacmi \([-1, 1]^{10}\) küpünde \(10^6\) noktalı Monte Carlo yöntemiyle kestirin ve sonucu, yaklaşık aynı sayıda (\(4^{10}\)) nokta kullanan orta nokta ızgarasının verdiği değerle karşılaştırın.
Çözüm
Topun hacmi, topun içinde \(1\), dışında \(0\) olan fonksiyonun küp üzerindeki integralidir. Küpün hacmi \(2^{10} = 1024\) olduğundan Monte Carlo kestirimi, topa düşen noktaların oranının \(1024\) katıdır. Izgarada ise her eksen \([-1, 1]\) aralığının dört eşit parçasının orta noktalarına, yani \(\pm 0{,}25\) ve \(\pm 0{,}75\)’e bölünür. Her ızgara noktası hacmi \((2/4)^{10} = 1/1024\) olan bir küçük küpü temsil eder.
import math
import numpy as np
d, N = 10, 10**6
rng = np.random.default_rng(19)
P = rng.uniform(-1, 1, size=(N, d)) # [-1, 1]^10 küpü, hacmi 2^10
inside = (P**2).sum(axis=1) <= 1
V_hat = 2**d * inside.mean()
se = 2**d * inside.std(ddof=1) / math.sqrt(N)
print(inside.sum(), V_hat, se)
print(math.pi**5 / 120) # tam değer
# Orta nokta kuralı: her eksende m = 4 nokta, toplam 4^10 nokta
m = 4
c = -1 + (np.arange(m) + 0.5) * 2 / m # -0.75, -0.25, 0.25, 0.75
G = np.stack(np.meshgrid(*[c] * d, indexing="ij"), axis=-1)
G = G.reshape(-1, d)
grid_in = (G**2).sum(axis=1) <= 1
print(len(G), grid_in.sum(), grid_in.sum() * (2 / m) ** d)Çıktı:
2480 2.53952 0.05093154142760173
2.550164039877345
1048576 1024 1.0Izgarayı kuran satırda [c] * d, c dizisinin 10 kopyasından oluşan bir listedir. Başındaki * bu listeyi açar ve 10 elemanını meshgrid’e ayrı ayrı argüman olarak verir; yani np.meshgrid(c, c, ..., c) yazmış oluruz. meshgrid her biri \((4, 4, \ldots, 4)\) şekilli 10 koordinat dizisi döndürür. np.stack(..., axis=-1) bunları yeni bir son eksen boyunca yığar ve \((4, \ldots, 4, 10)\) şekilli bir dizi kurar. reshape(-1, d) ise bu diziyi \((4^{10}, 10)\) şekline getirir: her satır bir ızgara noktasının 10 koordinatıdır.
Bir milyon noktanın yalnız \(2480\)’i topun içine düştü; topun küpteki payı \(2{,}55/1024 \approx 0{,}0025\) kadardır. Kestirim \(2{,}540\), standart hatası \(0{,}051\)’dir ve tam değer \(2{,}550\) güven aralığının içindedir. Izgara ise \(1{,}0\) veriyor, yani gerçek değerin yarısından az. Nedeni basittir: bir ızgara noktasının koordinatlarının karelerinin toplamı en az \(10 \cdot 0{,}25^2 = 0{,}625\)’tir ve tek bir koordinat \(\pm 0{,}75\) olursa toplam \(0{,}625 + 0{,}5 = 1{,}125 > 1\) olur. Topun içine yalnız bütün koordinatları \(\pm 0{,}25\) olan \(2^{10} = 1024\) nokta düşer. Her eksende yalnız dört nokta kullanabilen bir ızgara, topun biçimini göremez. Daha iyi bir ızgara için her eksende 10 nokta bile \(10^{10}\) fonksiyon değeri demektir. Topun hacmi için genel formül \(\pi^{d/2}/\Gamma(d/2 + 1)\)’dir; \(d = 10\) için \(\pi^5/5! = \pi^5/120\) olur. \(\blacksquare\)
14.9 Rastgele Yürüyüş
Şimdiye kadar bağımsız rastgele sayıları tek tek kullandık. Onları art arda toplayınca zamanla değişen rastgele bir büyüklük elde ederiz. Bunun en basit örneği rastgele yürüyüştür.
Tanım 14.5 (Basit Simetrik Rastgele Yürüyüş) \(X_1, X_2, \ldots\) bağımsız ve her biri \(1\) ve \(-1\) değerlerini \(\frac{1}{2}\)’şer olasılıkla alan rastgele değişkenler olsun. \(S_0 = 0\) ve \(n \ge 1\) için \(S_n = X_1 + \cdots + X_n\) ile tanımlanan \((S_n)\) dizisine basit simetrik rastgele yürüyüş (simple symmetric random walk) denir.
Yani sayı doğrusunda \(0\)’dan başlayan bir yürüyücü her adımda yazı tura atar; tura gelirse bir birim sağa, yazı gelirse bir birim sola gider. \(S_n\), \(n\) adım sonraki konumdur. Bir sonraki konum yalnız şimdiki konuma ve yeni adıma bağlı olduğundan \((S_n)\) bir Markov zinciridir (bkz. Raslantı Süreçleri). Beş bin yürüyüşü 400’er adımla tek bir diziyle simüle edip ilk sekizini çizelim:
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(20)
m, n = 5000, 400 # 5000 yürüyüş, her biri 400 adım
steps = rng.choice([-1, 1], size=(m, n))
S = np.zeros((m, n + 1), dtype=int)
S[:, 1:] = steps.cumsum(axis=1) # S[:, 0] = 0 başlangıç
k = np.arange(n + 1)
fig, ax = plt.subplots(figsize=(8, 4.5), layout="constrained")
for path in S[:8]: # ilk sekiz yürüyüş
ax.plot(k, path, linewidth=1)
for c in (1, 2):
ax.plot(k, c * np.sqrt(k), "k--", k, -c * np.sqrt(k), "k--")
ax.set_xlabel("n")
ax.set_ylabel("$S_n$")
fig.savefig("yuruyus.png", dpi=150)
print(S[0, :11])
end = S[:, n] # 400. adımdaki konumlar
print(end.mean(), (end**2).mean()) # teori: 0 ve 400
print(np.mean(np.abs(end) <= 2 * np.sqrt(n)))Çıktı:
[ 0 1 0 -1 -2 -1 -2 -3 -2 -3 -4]
0.0916 395.6552
0.9614rng.choice([-1, 1], size=(m, n)) her satırı bir yürüyüşün adımları olan bir dizi üretir; cumsum(axis=1) adımları satır boyunca toplayarak konumları verir. S dizisinin ilk sütunu başlangıç konumu \(S_0 = 0\)’dır. Son üç satır 400. adımdaki konumların ortalamasını, karelerinin ortalamasını ve \(|S_{400}| \le 40\) olan yürüyüşlerin oranını basıyor.
Önerme 14.3 (Yürüyüşün Ortalaması ve Varyansı) Basit simetrik rastgele yürüyüşte her \(n\) için
\[ E(S_n) = 0, \qquad \operatorname{Var}(S_n) = E(S_n^2) = n \]
olur ve \((S_n + n)/2\) rastgele değişkeni \(B(n, \frac{1}{2})\) dağılımlıdır.
Aşağıdaki şekil yukarıdaki kodun çizdiği sekiz yürüyüştür; kesikli eğriler \(\pm\sqrt{n}\) ve \(\pm 2\sqrt{n}\)’dir.
İspat
Her adım için \(E(X_i) = \frac{1}{2} \cdot 1 + \frac{1}{2} \cdot (-1) = 0\) ve \(E(X_i^2) = 1\), dolayısıyla \(\operatorname{Var}(X_i) = 1\)’dir. Beklenen değerin lineerliğinden \(E(S_n) = 0\) olur. Bağımsız değişkenlerin toplamının varyansı varyansların toplamı olduğundan (bkz. Olasılık Teorisi; iki terimli eşitlik tümevarımla \(n\) terime genişler) \(\operatorname{Var}(S_n) = n\)’dir. \(E(S_n) = 0\) olduğundan \(E(S_n^2) = \operatorname{Var}(S_n) = n\) olur.
\(B_i = (X_i + 1)/2\) değişkenleri bağımsızdır ve \(1\) ile \(0\) değerlerini \(\frac{1}{2}\)’şer olasılıkla alır, yani \(B(1, \frac{1}{2})\) dağılımlıdır. Toplamları \(\sum_{i=1}^{n} B_i = (S_n + n)/2\) olduğundan bağımsız Bernoulli değişkenlerinin toplamı hakkındaki önerme (bkz. Olasılık Teorisi) \((S_n + n)/2 \sim B(n, \frac{1}{2})\) verir. \(\blacksquare\)
Yani yürüyücü ortalamada başladığı yerdedir, ama başlangıçtan tipik uzaklığı \(\sqrt{E(S_n^2)} = \sqrt{n}\)’dir: uzaklık adım sayısıyla değil, onun kareköküyle büyür. Simülasyonda 400. adımdaki konumların ortalaması \(0{,}09\), karelerinin ortalaması \(395{,}7\)’dir; kuramsal değerler \(0\) ve \(400\)’dür. Merkezi limit teoremine göre \(S_n/\sqrt{n}\) yaklaşık \(N(0, 1)\) dağılımlıdır, bu yüzden yürüyüşlerin yaklaşık %95’i \(\pm 2\sqrt{n}\) bandında biter. Önermenin binom kısmından tam değer de hesaplanabilir: \(|S_{400}| \le 40\) olması \(180 \le (S_{400} + 400)/2 \le 220\) demektir ve bu olasılık \(0{,}9598\)’dir. Simülasyon \(0{,}9614\) buldu.
Örnek 14.7 (Kumarbazın İflası) Bir oyuncu 3 lirayla başlıyor ve her turda adil bir yazı tura oyununda 1 lira kazanıyor ya da kaybediyor. Parası \(0\)’a düşünce ya da \(10\) liraya ulaşınca oyun bitiyor. \(20\,000\) oyun simüle ederek oyuncunun 10 liraya ulaşma olasılığını kestirin.
Çözüm
Oyuncunun parası, \(3\)’ten başlayan ve \(0\) ya da \(10\)’a değince duran bir rastgele yürüyüştür. Oyunlar farklı sürelerde bittiği için sabit uzunlukta bir dizi kullanamayız. Bunun yerine bütün oyunları birlikte ilerletiyor ve hâlâ süren oyunları bir boolean maskeyle izliyoruz: her turda yalnız maskenin doğru olduğu oyunlara bir adım ekleniyor, maske de yeniden hesaplanıyor. Döngü bütün oyunlar bitince duruyor.
import numpy as np
rng = np.random.default_rng(21)
games, k, N = 20_000, 3, 10
pos = np.full(games, k) # her oyunda başlangıç sermayesi
steps = np.zeros(games, dtype=int) # oyunların süresi
active = np.ones(games, dtype=bool) # hâlâ süren oyunlar
while active.any():
pos[active] += rng.choice([-1, 1], size=active.sum())
steps[active] += 1
active = (0 < pos) & (pos < N)
print((pos == N).mean(), k / N)
print(steps.mean(), k * (N - k))
print(steps.max())Çıktı:
0.2997 0.3
21.1472 21
265Oyunların \(\%29{,}97\)’si 10 lirayla bitti. Tam değer \(k/N = 3/10\)’dur. Bunu görmek için \(p_k\), \(k\) lirayla başlayan oyuncunun 10’a ulaşma olasılığı olsun. İlk atışa göre koşullayınca \(0 < k < 10\) için
\[ p_k = \tfrac{1}{2}\, p_{k-1} + \tfrac{1}{2}\, p_{k+1}, \qquad p_0 = 0, \quad p_{10} = 1 \]
olur. Bu, \(p_{k+1} - p_k = p_k - p_{k-1}\) demektir: ardışık farklar sabittir, \(p_k\) doğrusal büyür ve sınır koşulları \(p_k = k/10\) verir. Aynı simülasyon oyunların süresini de veriyor: ortalama süre \(21{,}15\) tur, en uzun oyun \(265\) tur. Ortalama süre \(d_k\) için aynı koşullama \(d_k = 1 + \frac{1}{2}(d_{k-1} + d_{k+1})\), \(d_0 = d_{10} = 0\) denklemlerini verir; çözümü \(d_k = k(10 - k)\)’dir ve \(k = 3\) için \(21\) olur. Bu iki sonuç, durumu yutucu bir Markov zinciri olarak yazıp temel matrisi hesaplayarak da bulunabilir (bkz. Raslantı Süreçleri). \(\blacksquare\)
14.10 Alıştırmalar
Aşağıdaki alıştırmalarda her kod sabit bir tohumla çalışır; kendi denemelerinizde tohumu değiştirip sonuçların ne kadar oynadığına da bakın.
Alıştırma 14.1 (İki Zarın Toplamının Dağılımı) default_rng(31) üreteciyle iki zarı \(360\,000\) kez atın. Toplamın \(2, 3, \ldots, 12\) değerlerini alma sıklıklarını np.bincount ile bulun ve \(P(T = t) = \frac{6 - |t - 7|}{36}\) olasılıklarıyla karşılaştırın.
Çözüm
Önce formülü doğrulayalım: toplamı \(t\) olan \((a, b)\) çiftleri \(2 \le t \le 7\) için \((1, t - 1), \ldots, (t - 1, 1)\), yani \(t - 1\) tanedir; \(7 \le t \le 12\) için \(13 - t\) tanedir. İki durumda da sayı \(6 - |t - 7|\)’dir ve 36 eşit olasılıklı sonuca bölünür.
np.bincount(s, minlength=13) dizisinin \(j\). elemanı, s içinde \(j\) değerinin kaç kez geçtiğidir. \(0\) ve \(1\) toplamları imkânsız olduğundan [2:] ile atıyoruz.
import numpy as np
rng = np.random.default_rng(31)
N = 360_000
s = rng.integers(1, 7, N) + rng.integers(1, 7, N)
freq = np.bincount(s, minlength=13)[2:] / N # toplam 2, 3, ..., 12
exact = (6 - np.abs(np.arange(2, 13) - 7)) / 36
for total, fr, ex in zip(range(2, 13), freq, exact):
print(f"{total:2d} {fr:.4f} {ex:.4f}")
print(np.abs(freq - exact).max())Çıktı:
2 0.0275 0.0278
3 0.0553 0.0556
4 0.0837 0.0833
5 0.1115 0.1111
6 0.1391 0.1389
7 0.1656 0.1667
8 0.1390 0.1389
9 0.1107 0.1111
10 0.0837 0.0833
11 0.0560 0.0556
12 0.0280 0.0278
0.001063888888888892Sıklıklar olasılıklardan en çok \(0{,}0011\) sapıyor. Olasılığı \(p = 1/6\) olan bir olayın \(N = 360\,000\) denemedeki göreli sıklığının standart hatası \(\sqrt{p(1 - p)/N} \approx 0{,}0006\)’dır; sapmalar bu ölçektedir. \(\blacksquare\)
Alıştırma 14.2 (Hileli Zarı Kümülatif Toplamla Atmak) Bir hileli zarda \(1, 2, 3, 4, 5\) yüzlerinin her biri \(0{,}1\), \(6\) yüzü ise \(0{,}5\) olasılıkla geliyor. default_rng(32) ile ürettiğiniz \(10^5\) düzgün sayıyı, olasılıkların kümülatif toplamı ve np.searchsorted yardımıyla zar atışlarına çevirin; yüzlerin sıklıklarını ve ortalamayı kontrol edin.
Çözüm
Bu, kesikli bir dağılım için ters dönüşüm yöntemidir (Önerme 14.1). Kümülatif olasılıklar \(0{,}1;\ 0{,}2;\ 0{,}3;\ 0{,}4;\ 0{,}5;\ 1\)’dir. \(u\) sayısı \([0;\ 0{,}1)\) aralığındaysa \(1\), \([0{,}1;\ 0{,}2)\) aralığındaysa \(2\), …, \([0{,}5;\ 1)\) aralığındaysa \(6\) yüzünü seçeriz. Bu aralıkların uzunlukları tam olarak yüzlerin olasılıklarıdır. np.searchsorted(cdf, u, side="right") her \(u\) için \(u < \text{cdf}[j]\) olan ilk \(j\) indeksini döndürür; bu indeks de faces dizisinde yüzü seçer. Karşılaştırma için hazır yol olan rng.choice(..., p=p) sonucunu da yazdırıyoruz.
import numpy as np
faces = np.arange(1, 7)
p = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.5])
cdf = np.cumsum(p)
print(cdf)
rng = np.random.default_rng(32)
N = 100_000
u = rng.random(N)
idx = np.searchsorted(cdf, u, side="right") # u < cdf[j] olan ilk j
x = faces[idx]
print(np.bincount(x, minlength=7)[1:] / N)
print(x.mean(), (faces * p).sum())
y = rng.choice(faces, size=N, p=p) # hazır yol
print(np.bincount(y, minlength=7)[1:] / N)Çıktı:
[0.1 0.2 0.3 0.4 0.5 1. ]
[0.09994 0.09821 0.09993 0.10118 0.10079 0.49995]
4.50452 4.5
[0.09925 0.09895 0.10062 0.09792 0.09869 0.50457]Sıklıklar olasılıklara yakındır ve ortalama \(4{,}5045\), kuramsal değer \(0{,}1 \cdot 15 + 0{,}5 \cdot 6 = 4{,}5\)’tir. rng.choice aynı dağılımı verir. \(\blacksquare\)
Alıştırma 14.3 (Ters Dönüşümle Doğrusal Bir Yoğunluktan Örneklem) \([0, 1]\) aralığında yoğunluğu \(f(x) = 2x\) olan dağılımdan, default_rng(33) ile ters dönüşüm yöntemini kullanarak \(10^5\) elemanlı bir örneklem üretin. Örneklemin ortalamasını, varyansını ve \(P(X \le 1/2)\) kestirimini tam değerlerle karşılaştırın.
Çözüm
Üç adımı uyguluyoruz.
- \(0 \le x \le 1\) için \(F(x) = \int_0^x 2t\,dt = x^2\)’dir.
- \(u = x^2\) denkleminin \([0, 1]\)’deki çözümü \(x = \sqrt{u}\)’dur, yani \(F^{-1}(u) = \sqrt{u}\).
np.sqrt(rng.random(n))örneklemi verir.
Tam değerler: \(E(X) = \int_0^1 2x^2\,dx = \frac{2}{3}\), \(E(X^2) = \int_0^1 2x^3\,dx = \frac{1}{2}\), dolayısıyla \(\operatorname{Var}(X) = \frac{1}{2} - \frac{4}{9} = \frac{1}{18}\); ayrıca \(P(X \le \frac{1}{2}) = F(\frac{1}{2}) = \frac{1}{4}\).
import numpy as np
rng = np.random.default_rng(33)
n = 100_000
x = np.sqrt(rng.random(n)) # F^{-1}(u) = kök u
print(x.mean(), 2 / 3)
print(x.var(ddof=1), 1 / 18)
print(np.mean(x <= 0.5), 0.25)Çıktı:
0.666972682598439 0.6666666666666666
0.05534262966614402 0.05555555555555555
0.25007 0.25Üç kestirim de tam değerlerle üç basamak uyuşuyor. \(\blacksquare\)
Alıştırma 14.4 (İlk Altıya Kadar Atış Sayısı) Bir zar ilk kez 6 gelene kadar atılıyor ve \(T\) atış sayısı olsun. Bu deneyi default_rng(34) ile bir while döngüsü kullanarak \(10^5\) kez simüle edin; \(T\)’nin ortalamasını, varyansını ve \(P(T > 10)\) olasılığını kestirip tam değerlerle karşılaştırın.
Çözüm
Tek bir deneyi bir fonksiyon yapıyor: zar 6 gelmediği sürece sayacı artırıyor. \(T\), başarı olasılığı \(p = 1/6\) olan bağımsız denemelerde ilk başarıya kadar yapılan deneme sayısıdır, yani \(\operatorname{Geo}(1/6)\) dağılımlıdır. Bu dağılımda \(E(T) = 1/p = 6\), \(\operatorname{Var}(T) = q/p^2 = 30\) ve \(P(T > k) = q^k\) olur, burada \(q = 5/6\)’dır (bkz. Olasılık Teorisi). Son iki satırda aynı dağılımdan hazır rng.geometric metoduyla da örneklem alıyoruz.
import numpy as np
def rolls_until_six(rng):
"""İlk 6 gelene kadar atılan zar sayısı."""
count = 1
while rng.integers(1, 7) != 6:
count += 1
return count
rng = np.random.default_rng(34)
T = np.array([rolls_until_six(rng) for _ in range(100_000)])
print(T[:10])
print(T.mean(), T.var(ddof=1)) # teori: 6 ve 30
print(np.mean(T > 10), (5 / 6) ** 10) # P(T > 10) = (5/6)^10
G = rng.geometric(1 / 6, size=100_000) # hazır geometrik dağılım
print(G.mean(), G.var(ddof=1))Çıktı:
[ 4 10 2 4 2 7 1 2 1 1]
6.00385 30.22769745447455
0.16042 0.1615055828898458
6.00523 30.212884775947757Ortalama \(6{,}004\), varyans \(30{,}23\) ve \(P(T > 10)\) kestirimi \(0{,}1604\)’tür; tam değerler \(6\), \(30\) ve \((5/6)^{10} \approx 0{,}1615\)’tir. Döngülü simülasyon yavaş ama anlaşılırdır; rng.geometric aynı işi tek çağrıda yapar. \(\blacksquare\)
Alıştırma 14.5 (Monty Hall Problemi) Bir yarışmada üç kapıdan birinin arkasında araba, ikisinin arkasında keçi vardır. Yarışmacı bir kapı seçer. Arabanın yerini bilen sunucu, kalan iki kapıdan arkasında keçi olan birini açar; böyle iki kapı varsa birini rastgele seçer. Yarışmacı kapısını kapalı kalan diğer kapıyla değiştirirse arabayı kazanma olasılığı nedir? default_rng(35) ile \(10^5\) oyun simüle ederek kestirin.
Çözüm
Kapıları \(0, 1, 2\) diye numaralayalım. Her oyunda arabanın kapısını ve yarışmacının ilk seçimini düzgün dağılımla seçiyoruz. Sunucunun açabileceği kapılar, ne seçilen ne de arabalı olan kapılardır; bunlardan birini rng.choice ile seçiyoruz. Üç kapının numaralarının toplamı \(3\) olduğundan kapalı kalan diğer kapı 3 - pick - host’tur.
import numpy as np
rng = np.random.default_rng(35)
N = 100_000
car = rng.integers(0, 3, N) # arabanın arkasındaki kapı
pick = rng.integers(0, 3, N) # yarışmacının ilk seçimi
host = np.empty(N, dtype=int) # sunucunun açtığı kapı
for i in range(N):
options = [d for d in range(3) if d != pick[i] and d != car[i]]
host[i] = rng.choice(options)
other = 3 - pick - host # kalan kapı (0 + 1 + 2 = 3)
print(car[:6], pick[:6], host[:6], other[:6])
print(np.mean(pick == car), np.mean(other == car))Çıktı:
[0 1 2 1 2 2] [0 1 0 2 1 0] [2 2 1 0 0 1] [1 0 2 1 2 2]
0.33427 0.66573Kapıyı değiştirmeyen yarışmacı oyunların \(\%33{,}4\)’ünü, değiştiren \(\%66{,}6\)’sını kazanıyor. Tam değerler \(1/3\) ve \(2/3\)’tür. Nedeni şudur: değiştiren yarışmacı, ilk seçimi yanlışsa her zaman kazanır, çünkü sunucu öteki keçiyi açmıştır ve kalan kapı arabalıdır. İlk seçimin yanlış olma olasılığı \(2/3\)’tür. \(\blacksquare\)
Alıştırma 14.6 (Binoma Normal Yaklaşım) \(X \sim B(100, \frac{1}{2})\) için \(P(X \le 55)\) olasılığını binom toplamıyla tam olarak hesaplayın. Bu değeri default_rng(36) ile \(10^6\) simülasyonun kestirimiyle ve süreklilik düzeltmeli normal yaklaşımla karşılaştırın.
Çözüm
Tam değer \(P(X \le 55) = \sum_{j=0}^{55} \binom{100}{j} 2^{-100}\)’dür; bunu math.comb ile kesin tamsayılarla hesaplayabiliriz. \(X\), 100 bağımsız Bernoulli değişkeninin toplamı olduğundan merkezi limit teoremine göre yaklaşık \(N(\mu, \sigma^2)\) dağılımlıdır; burada \(\mu = np = 50\) ve \(\sigma = \sqrt{np(1 - p)} = 5\)’tir. \(X\) tamsayı değerli olduğundan \(P(X \le 55) = P(X < 55{,}5)\)’tir; normal yaklaşımda \(55\) yerine \(55{,}5\) kullanmaya süreklilik düzeltmesi denir (bkz. Matematiksel İstatistik).
import math
import numpy as np
def Phi(z):
"""Standart normal dağılım fonksiyonu."""
return 0.5 * (1 + math.erf(z / math.sqrt(2)))
n, p = 100, 0.5
exact = sum(math.comb(n, j) for j in range(56)) / 2**n # P(X <= 55)
rng = np.random.default_rng(36)
X = rng.binomial(n, p, size=10**6)
mu, sd = n * p, math.sqrt(n * p * (1 - p))
print(exact)
print(np.mean(X <= 55))
print(Phi((55 - mu) / sd), Phi((55.5 - mu) / sd))Çıktı:
0.8643734879630827
0.864361
0.8413447460685428 0.8643339390536173Tam değer \(0{,}86437\), simülasyon \(0{,}86436\)’dır. Düzeltmesiz normal yaklaşım \(\Phi(1) \approx 0{,}8413\) değeriyle yaklaşık \(0{,}023\) yanılıyor; düzeltmeli yaklaşım \(\Phi(1{,}1) \approx 0{,}86433\) ise dört basamak doğrudur. \(\blacksquare\)
Alıştırma 14.7 (Bir İntegralle Pi Sayısını Kestirmek) \(\int_0^1 \frac{4}{1 + x^2}\,dx = \pi\)’dir. Bu integrali Monte Carlo yöntemiyle \(10^{-4}\) standart hatayla hesaplamak için kaç nokta gerekir? default_rng(37) ile \(10^5\) noktalı bir ön deneme yapıp \(\sigma\)’yı kestirerek bulun.
Çözüm
İntegralin değeri \(4 \arctan 1 = \pi\)’dir. \(f(x) = 4/(1 + x^2)\) ve \([0, 1]\) aralığının uzunluğu \(1\) olduğundan Önerme 14.2 standart hatanın \(\sigma/\sqrt{N}\) olduğunu söyler. Ön denemenin örneklem standart sapması \(s\) ile \(\sigma\)’yı kestiririz; \(s/\sqrt{N} \le 10^{-4}\) koşulu \(N \ge (s/10^{-4})^2\) demektir.
import math
import numpy as np
rng = np.random.default_rng(37)
N = 10**5
y = 4 / (1 + rng.random(N) ** 2)
I_hat = y.mean()
s = y.std(ddof=1)
se = s / math.sqrt(N)
print(I_hat, se, abs(I_hat - math.pi))
print(s, math.ceil((s / 1e-4) ** 2)) # std. hata <= 1e-4 için gereken NÇıktı:
3.143196965868777 0.002028128191954597 0.0016043122789839437
0.6413504473375709 41133040Ön denemenin kestirimi \(3{,}14320\), standart hatası \(0{,}0020\)’dir; gerçek hata \(0{,}0016\) bu ölçektedir. \(s \approx 0{,}641\) olduğundan yaklaşık \(4{,}1 \cdot 10^7\) nokta gerekir. Bunu standart hatanın davranışından da görebiliriz: standart hatayı \(0{,}002\)’den \(0{,}0001\)’e, yani 20 kat küçültmek için nokta sayısını \(20^2 = 400\) katına çıkarmak gerekir. Kesin değer
\[ \sigma^2 = \int_0^1 \frac{16}{(1 + x^2)^2}\,dx - \pi^2 = 2\pi + 4 - \pi^2 \]
olduğundan \(\sigma \approx 0{,}643\)’tür. \(\blacksquare\)
Alıştırma 14.8 (Karşıt Değişkenlerle Varyansı Azaltmak) \(\int_0^1 e^x\,dx = e - 1\) integralini default_rng(38) ile iki yolla kestirin: \(N = 10^4\) noktalı sıradan Monte Carlo ile ve \(N/2\) tane düzgün \(V_j\) sayısı için \(\frac{1}{2}\big(e^{V_j} + e^{1 - V_j}\big)\) değerlerinin ortalamasıyla. İki yol da \(f\)’yi \(10^4\) kez hesaplar; hangisinin standart hatası daha küçüktür ve neden?
Çözüm
\(1 - V\) de \(U(0, 1)\) dağılımlı olduğundan \(E\big[\frac{1}{2}(e^V + e^{1-V})\big] = e - 1\)’dir; ikinci kestirim de doğru değeri hedefler. İkinci yöntemde birlikte kullanılan \(e^V\) ve \(e^{1-V}\) değerlerine karşıt değişkenler (antithetic variates) denir.
import math
import numpy as np
rng = np.random.default_rng(38)
N = 10_000 # f'nin hesaplanma sayısı
plain = np.exp(rng.random(N)) # sıradan Monte Carlo
v = rng.random(N // 2) # N/2 nokta, her biri iki kez
anti = (np.exp(v) + np.exp(1 - v)) / 2 # karşıt çiftlerin ortalaması
for name, y in [("sıradan", plain), ("karşıt", anti)]:
se = y.std(ddof=1) / math.sqrt(len(y))
print(f"{name:8}{y.mean():.5f} std. hata = {se:.2e}")
print(math.e - 1)Çıktı:
sıradan 1.71762 std. hata = 4.93e-03
karşıt 1.71803 std. hata = 8.84e-04
1.718281828459045Karşıt değişkenlerin standart hatası 5,6 kat küçüktür. Nedeni, \(e^V\) büyükken \(e^{1-V}\)’nin küçük olmasıdır: iki değerin kovaryansı negatiftir ve ortalamada hatalar birbirini kısmen götürür. Tam hesap:
\[ \begin{aligned} \operatorname{Var}(e^V) &= \frac{e^2 - 1}{2} - (e - 1)^2 \approx 0{,}2420, \\[1mm] \operatorname{Cov}(e^V, e^{1-V}) &= E[e^V e^{1-V}] - (e - 1)^2 = e - (e - 1)^2 \approx -0{,}2342. \end{aligned} \]
Bir çiftin ortalamasının varyansı \(\frac{1}{2}(\operatorname{Var} + \operatorname{Cov})\), yani \(\frac{1}{2}(0{,}2420 - 0{,}2342) \approx 0{,}0039\)’dur. Sıradan yöntemde \(10^4\) değerin ortalamasının varyansı \(0{,}2420/10^4\), karşıt yöntemde \(5000\) çiftin ortalamasının varyansı \(0{,}0039/5000\)’dir. Oran yaklaşık \(31\)’dir; standart hataların oranı da \(\sqrt{31} \approx 5{,}6\) olur. Aynı doğruluğa sıradan Monte Carlo ancak yaklaşık 31 kat fazla fonksiyon değeriyle ulaşır. \(\blacksquare\)
Alıştırma 14.9 (Diskte Orijine Ortalama Uzaklık) Birim diskte düzgün dağılımlı noktaları, default_rng(39) ile \([-1, 1]^2\) karesinde ürettiğiniz \(10^6\) noktadan diskin dışında kalanları atarak elde edin. Bu noktalarla bir noktanın orijine ortalama uzaklığını kestirin ve tam değerle karşılaştırın.
Çözüm
Karede düzgün dağılımlı bir noktanın, diske düştüğü bilindiğinde, diskte düzgün dağılımlı olduğunu kullanıyoruz: diskin içindeki her bölgenin koşullu olasılığı alanıyla orantılıdır. Bu yüzden diske düşmeyen noktaları atmak yeterlidir; bu yönteme reddetme yöntemi (rejection sampling) denir. Noktaların yaklaşık \(\pi/4\)’ü kalır.
import numpy as np
rng = np.random.default_rng(39)
N = 10**6
P = rng.uniform(-1, 1, size=(N, 2)) # [-1, 1]^2 karesinde noktalar
r = np.sqrt((P**2).sum(axis=1))
r = r[r <= 1] # diskin dışındakileri at
print(len(r), len(r) / N, np.pi / 4)
print(r.mean(), 2 / 3)
print(r.std(ddof=1) / np.sqrt(len(r)))Çıktı:
785471 0.785471 0.7853981633974483
0.6664194817138088 0.6666666666666666
0.0002660426055334039\(785\,471\) nokta kaldı. Ortalama uzaklık kestirimi \(0{,}66642\), standart hatası \(0{,}00027\)’dir. Tam değer, \(\sqrt{x^2 + y^2}\)’nin diskteki ortalama değeridir: kutupsal koordinatlarla \(\frac{1}{\pi}\int_0^{2\pi}\!\int_0^1 r \cdot r\,dr\,d\theta = \frac{2}{3}\) (bkz. İntegral Calculus). Fark \(0{,}00025\), bir standart hatadan küçüktür. \(\blacksquare\)
Alıştırma 14.10 (Kumarbazın İflası İçin Temel Matris) Kumarbazın İflası örneğindeki oyunu (Örnek 14.7) \(0, 1, \ldots, 10\) durumlu ve \(0\) ile \(10\) durumları yutucu olan bir Markov zinciri olarak yazın. \(Q\) ve \(R\) matrislerini kurup \((I - Q)\,t = \mathbf{1}\) ve \((I - Q)\,B = R\) sistemlerini np.linalg.solve ile çözerek, simülasyonun 3 liralık başlangıç için bulduğu \(0{,}2997\) ve \(21{,}15\) değerlerini tam değerlerle karşılaştırın.
Çözüm
Geçici durumlar \(1, \ldots, 9\), yutucu durumlar \(0\) ve \(10\)’dur. \(Q\), geçici durumlar arasındaki geçiş olasılıklarının \(9 \times 9\) matrisidir: \(i\) durumundan \(i - 1\) ve \(i + 1\) durumlarına \(\frac{1}{2}\) olasılıkla gidilir. \(R\), geçici durumlardan yutucu durumlara geçişlerin \(9 \times 2\) matrisidir; yalnız \(1 \to 0\) ve \(9 \to 10\) geçişleri sıfırdan farklıdır. Temel matris \(N = (I - Q)^{-1}\)’in \((i, j)\) elemanı, \(i\)’den başlayan zincirin yutulmadan önce \(j\) durumunda geçirdiği ortalama adım sayısıdır. Satır toplamları \(t = N\mathbf{1}\) ortalama süreleri, \(B = NR\) ise yutulma olasılıklarını verir (bkz. Raslantı Süreçleri). Ters matrisi açıkça hesaplamak yerine aynı sonuçları solve ile buluyoruz (bkz. NumPy ile Lineer Cebir).
import numpy as np
goal = 10 # durumlar 0, 1, ..., 10
T = goal - 1 # geçici durumlar 1, ..., 9
Q = np.zeros((T, T)) # i. satır: i + 1 durumu
for i in range(T - 1):
Q[i, i + 1] = Q[i + 1, i] = 0.5 # bir adım sağa ya da sola
R = np.zeros((T, 2)) # sütunlar: 0 ve 10 durumları
R[0, 0] = R[-1, 1] = 0.5
A = np.eye(T) - Q # I - Q
t = np.linalg.solve(A, np.ones(T)) # (I - Q)^{-1} 1: beklenen süre
B = np.linalg.solve(A, R) # (I - Q)^{-1} R: yutulma
print(t)
print(B[:, 1])
print(t[2], B[2, 1]) # başlangıç sermayesi 3Çıktı:
[ 9. 16. 21. 24. 25. 24. 21. 16. 9.]
[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]
21.0 0.2999999999999999Ortalama süreler \(t_k = k(10 - k)\), 10’a ulaşma olasılıkları \(k/10\) çıkıyor; Kumarbazın İflası örneğinde fark denklemleriyle bulduğumuz formüllerin aynısı. \(k = 3\) için tam değerler \(21\) ve \(0{,}3\)’tür (son basamaktaki \(0{,}2999999999999999\) yuvarlama hatasıdır). Simülasyonun \(21{,}15\) ve \(0{,}2997\) değerleri bunlara yakındır. \(\blacksquare\)
Alıştırma 14.11 (Düzlemde Rastgele Yürüyüş) Düzlemde orijinden başlayan bir yürüyücü her adımda sağa, sola, yukarı ya da aşağı, eşit olasılıkla bir birim gidiyor. default_rng(41) ile \(10^4\) yürüyüşü \(100\)’er adımla simüle edin ve 100 adım sonra orijine uzaklığın karesinin ortalamasını kuramsal değeriyle karşılaştırın.
Çözüm
Dört yönü dirs dizisinin satırlarına yazıyoruz. rng.integers(0, 4, size=(m, n)) her adım için bir yön indeksi seçer ve dirs[idx] gelişmiş indeksleme ile (m, n, 2) şekilli adım dizisini kurar (bkz. NumPy Dizileri). Adımları axis=1 boyunca toplayınca son konumları buluruz.
Kuramsal değer şöyle bulunur. Adım vektörleri \(X_1, \ldots, X_n\) ise \(|S_n|^2 = \sum_{i} \sum_{j} X_i \cdot X_j\)’dir. \(i \ne j\) için \(X_i\) ile \(X_j\) bağımsızdır ve her birinin beklenen değeri sıfır vektörüdür; bu yüzden \(E(X_i \cdot X_j) = E(X_i) \cdot E(X_j) = 0\) olur. \(i = j\) için \(X_i \cdot X_i = 1\)’dir. Dolayısıyla \(E|S_n|^2 = n = 100\)’dür.
import math
import numpy as np
rng = np.random.default_rng(41)
m, n = 10_000, 100
dirs = np.array([[1, 0], [-1, 0], [0, 1], [0, -1]])
idx = rng.integers(0, 4, size=(m, n)) # her adımda bir yön
steps = dirs[idx] # şekil (m, n, 2)
pos = steps.sum(axis=1) # n adım sonraki konum
d2 = (pos**2).sum(axis=1) # uzaklığın karesi
print(steps.shape, pos.shape)
print(d2.mean(), n)
print(np.mean(d2 == 0), (math.comb(n, n // 2) / 2**n) ** 2)Çıktı:
(10000, 100, 2) (10000, 2)
100.7924 100
0.0058 0.006334446707872695Kestirim \(100{,}79\)’dur. Son satır, yürüyüşlerin 100. adımda orijine dönme oranını \((\binom{100}{50} 2^{-100})^2 \approx 0{,}0063\) tam değeriyle karşılaştırıyor. Bu formül şöyle çıkar: \(u = x + y\) ve \(v = x - y\) koordinatlarında her adım \(u\)’yu ve \(v\)’yi \(\pm 1\) değiştirir ve dört yön dört işaret çiftine karşılık gelir. Dolayısıyla \(u\) ve \(v\) iki bağımsız basit simetrik yürüyüştür ve orijine dönmek ikisinin de \(0\)’da olması demektir. Tek boyutta \(P(S_{100} = 0) = \binom{100}{50} 2^{-100}\) olduğu ise yürüyüşün binom dağılımıyla ilişkisinden (Önerme 14.3) gelir. \(\blacksquare\)
Bu bölümde sözde rastgele sayı üreteçlerini, default_rng ile tekrarlanabilir simülasyon yazmayı ve temel dağılımlardan, gerekirse ters dönüşümle, örneklem üretmeyi gördük. Büyük sayılar yasası Monte Carlo kestirimlerinin doğru değere gittiğini, merkezi limit teoremi ise hatalarının \(\sigma/\sqrt{N}\) ölçekli ve yaklaşık normal olduğunu söyledi; bu da her kestirime bir standart hata ve bir güven aralığı eklememizi sağladı. Sıradaki Pandas ile Veri Analizi bölümünde, simülasyonların ve ölçümlerin ürettiği tablo biçimindeki verileri okumayı, düzenlemeyi, gruplamayı ve özetlemeyi öğreneceğiz.