12.3 Várható érték és Monte Carlo közelítés

12.3.1 Cél és konkrét helyzet

Egy ajánlómodell szimulációs önellenőrzésében a relevanciapontszámot $X \sim U(0,1)$-nek tekintjük, a modell jutalomfüggvénye pedig $g(x)=\sin(x)$ (radiánban mért, dimenzió nélküli érték). Ez mesterséges didaktikus feladat: a pontos várható érték $E[\sin X]=\int_0^1\sin x\,dx=1-\cos 1 \approx 0{,}4596977$. Ezt az ismert integrált Monte Carlo-mintavétellel becsüljük meg, és a becslési hibát több mintanagyságon, több különböző seeddel vizsgáljuk.

Bemenet: $g(x)=\sin(x)$, $x\in[0,1]$; $X\sim U(0,1)$; seed 0 és 199 közötti egész; $n\in\{100,400,1600,6400,25600\}$. Megfigyelési egység: egy $X$ realizáció. Kimenet: $\hat I_n=\frac1n\sum_{i=1}^n g(x_i)$; a hibát az ismert integrálhoz viszonyítjuk. Előzetes becslés: a becslés 0,46 körül mozog; $n=100$-nál a tipikus hiba nagyságrendje 0,025, ezért 0,435–0,485 közötti értékek gyakoriak.

12.3.2 Fogalmi magyarázat

A várható érték (expected value) a sűrűségfüggvény szerint súlyozott átlag. Folytonos $X$ esetén $$E[g(X)] = \int_{\mathbb R} g(x) f(x)\,dx,$$ ahol $f$ az $X$ sűrűségfüggvénye (density), $g(x)$ tetszőleges függvény, az integrál a teljes tartományon fut. Az $f(x)\,dx$ az $x$ körüli infinitezimális valószínűségi tömeg; az értelmezés tehát nem a $g$ egyszerű átlaga az $x$ tengely mentén, hanem az eloszlás szerint súlyozott átlag. Diszkrét esetben az integrált összeg helyettesíti: $E[g(X)]=\sum_k g(x_k)P(X=x_k)$. Az egyenletes $U(0,1)$ eloszlásnál $f(x)=1$ a $[0,1]$ intervallumon, máshol nulla, így $E[g(X)]=\int_0^1 g(x)\,dx$; itt a súlyozás láthatatlan, általános eloszlásnál azonban nem az.

A Monte Carlo-becslés ezt a területet véletlen pontokkal közelíti: ha $x_1,\dots,x_n$ független minta az $f$ eloszlásból, akkor $$\hat I_n = \frac1n\sum_{i=1}^n g(x_i).$$ A generált pontok sűrűsége hosszú távon az $f$ alakját követi, ezért az egyszerű mintaátlag a $\int f g$ integrált becsli. A Sűrűség és eloszlásfüggvény példája szerint a sűrűség magasság, a valószínűség pedig terület; itt a területet véletlen pontok átlagával mérjük. A formális bizonyítás a nagy számok törvénye, amelyet később, a Nagy számok és centrális határeloszlás alfejezetben vizsgálunk; most numerikus tapasztalatot szerzünk.

A seed a pszeudorandom generátor (pseudorandom number generator) kiinduló állapota. A NumPy modern API-ja a np.random.default_rng(seed) hívással hoz létre Generator objektumot; ettől a ponttól kezdve a generált sorozat determinisztikus. Ugyanaz a seed ugyanazt a sorozatot adja ugyanazon NumPy-verzión, de a bitazonosság nem garantált verziók vagy platformok között. A np.random.seed globális állapotot állít; tanulási projektben a lokális Generator biztonságosabb, mert nem szól bele más kód véletlen folyamába. A Polars lazy terv és streaming korlát lazy gondolatát felidézve: a default_rng(seed) még csak terv, a tényleges végrehajtás a random() vagy normal() hívás.

A Monte Carlo-hiba (Monte Carlo error) annak a szórása, amennyire $\hat I_n$ a pontos érték körül ingadozik. Egyetlen mintából a standard hiba (standard error, SE) becsülhető: $$s_n = \sqrt{\frac{1}{n-1}\sum_{i=1}^n\left(g(x_i)-\hat I_n\right)^2},\qquad \widehat{SE} = \frac{s_n}{\sqrt n}.$$ Itt $s_n$ a mintabeli szórás; a nevezőben $n-1$ szerepel, mert a mintaátlag eltávolít egy szabadsági fokot. A $\widehat{SE}$ azt mutatja, mekkora a becslés tipikus véletlen hibája. A képletben $g(x_i)$ shape-je $(n,)$, az összeg skalár; az $n$ mintaméret növelésével a hiba várhatóan $\sqrt n$-nel arányosan csökken, nem $n$-nel: négyszer akkora mintához körülbelül fele akkora hiba tartozik.

A független ismétlés (independent replication) több seeddel több külön becslést ad: $R$ darab, egyenként $n$ elemű becslést készítünk. A becslések $\hat I_n^{(1)},\dots,\hat I_n^{(R)}$ hibáiból a tapasztalati RMSE: $$RMSE_R = \sqrt{\frac1R\sum_{r=1}^R\left(\hat I_n^{(r)}-I\right)^2},$$ ahol $I=1-\cos 1$ a pontos érték. Ha a becslés torzítatlan, az RMSE a standard hiba empirikus párja; a sok seed éppen arra való, hogy az egyetlen realizáció mögötti bizonytalanságot láthatóvá tegye.

12.3.3 Teljes minimális munkapélda

A labor a lab/pyproject.toml alap numerikus csomagjait használja: NumPy és SciPy. A kódot mentsd mc_expected_value.py néven, majd futtasd a projektgyökérből:

python mc_expected_value.py

A parancs Windows PowerShell és Unix-szerű shell alatt is ugyanez. A kód:

import numpy as np
from scipy.integrate import quad

def mc_mean(g, rng, n):
    """Várható érték Monte Carlo-becslése X ~ U(0,1) esetén."""
    x = rng.random(size=n)          # shape: (n,), dtype: float64
    return float(np.mean(g(x)))     # g(x) shape: (n,), mean: skalár

exact = 1.0 - np.cos(1.0)           # E[sin(X)], X ~ U(0,1)
print("Numerikus kontroll quad-dal:",
      quad(lambda x: np.sin(x), 0, 1)[0])

sample_sizes = [100, 400, 1600, 6400, 25600]
repeats = 200

# Egy seedhez egy pszeudorandom folyam tartozik.
max_n = sample_sizes[-1]
x = np.vstack([
    np.random.default_rng(seed).random(max_n)
    for seed in range(repeats)
])                                   # shape: (repeats, max_n)
gx = np.sin(x)                       # shape: (repeats, max_n)

estimates = np.column_stack([
    gx[:, :n].mean(axis=1)
    for n in sample_sizes
])                                   # shape: (repeats, len(sample_sizes))

errors = estimates - exact
rmse_by_n = np.sqrt(np.mean(errors ** 2, axis=0))  # shape: (len(sample_sizes),)

demo = np.random.default_rng(0).random(1600)
demo_est = np.mean(np.sin(demo))
demo_sd = np.std(np.sin(demo), ddof=1)
demo_se = demo_sd / np.sqrt(len(demo))

print("sample_sizes:", sample_sizes)
print("rmse_by_n:", rmse_by_n)
print("single_seed n=1600 estimate, sd, SE:",
      demo_est, demo_sd, demo_se)

A x mátrix első dimenziója a seed, második dimenziója a mintaelem sorszáma. A gx[:, :n].mean(axis=1) minden seedhez az első $n$ elem átlagát adja, ezért minden mintaméretnél ugyanazt a folyamot használjuk. A lenti kimenet szemléltető, nem e szöveg ellenőrzött futtatása:

Numerikus kontroll quad-dal: 0.45969769413186023
sample_sizes: [100, 400, 1600, 6400, 25600]
rmse_by_n: [0.0248 0.0125 0.0061 0.0031 0.0016]
single_seed n=1600 estimate, sd, SE: 0.4601 0.248 0.0062

Kézi kis kontrollként vegyünk négy rögzített $u$-pontot: $0{,}10; 0{,}70; 0{,}35; 0{,}90$. A $\sin$ értékek összege körülbelül $0{,}09983+0{,}64422+0{,}34290+0{,}78333=1{,}87028$, az átlag $0{,}46757$. Ez egyetlen kis realizáció: a pontos $0{,}4597$ körül van, de $n=4$-nél a hiba $0{,}0079$; a kis minta nem bizonyíték.

12.3.4 Eredmény értelmezése

A quad értéke a pontos integrál, tehát a Monte Carlo-becslés célértéke. A rmse_by_n a 200 seedből számított empirikus hiba; azt várjuk, hogy minden négyszeres mintaméretnél megközelítőleg feleződjön. A fenti szemléltető sorozat ezt mutatja: $0{,}0248$, $0{,}0125$, $0{,}0061$, $0{,}0031$, $0{,}0016$. Az egyetlen seedből becsült $\widehat{SE}=s_n/\sqrt n$ ugyanezt a bizonytalanságot becsli egyetlen minta alapján; a példában $0{,}0062$ az $n=1600$-ra, ami közel van a $0{,}0061$ sok seedes RMSE-hez. Ez az egyezés maga is véletlen; más seeddel kissé eltér.

A három mennyiség – pontos integrál, mintaátlag, standard hiba – szerepe különbözik: a pontos integrál a paraméter, a mintaátlag a realizált becslés, a standard hiba pedig a becslési eljárás bizonytalanságát jellemzi. A $\widehat{SE}$ nem garantálja, hogy a pillanatnyi hiba pontosan ekkora; csak a tipikus nagyságrendet adja.

A szokásos Monte Carlo-standard hiba $s_g/\sqrt n$ független, azonos eloszlású és véges varianciájú mintákra alkalmazható. Korrelált szimulációknál, illetve nehéz farkú eloszlásoknál külön hibamodellre lehet szükség.

12.3.5 Feladatok

1. Értelmezés. Magyarázd el, miért közelíti a mc_mean függvényben szereplő np.mean(g(x)) az $E[g(X)]$ integrált, ha $X\sim U(0,1)$. Mi a szerepe a sűrűségfüggvénynek ebben az érvelésben?

2. Megvalósítás. Írj a mc_mean függvény felhasználásával szimulációt $g(x)=x^2$ függvényre, $X\sim U(0,1)$ mellett. A pontos várható érték $\int_0^1 x^2\,dx=1/3$. Használd az $n\in\{100,400,1600,6400,25600\}$ listát és $R=100$ seedet; minden $n$-re számold ki a hibák RMSE-ét, majd log-log skálán illesztett egyenes meredekségét. Milyen meredekséget vársz?

3. Hibadiagnózis. A következő kód abszolút hibákat rajzol az $n$ függvényében:

import numpy as np
import matplotlib.pyplot as plt

exact = 1.0 - np.cos(1.0)
sizes = np.logspace(2, 5, 50, dtype=int)
errs = []
for n in sizes:
    x = np.random.default_rng(7).random(n)
    errs.append(abs(np.mean(np.sin(x)) - exact))

plt.plot(sizes, errs)

Miért félrevezető ez az ábra a módszer bizonytalanságára nézve? Hogyan javítanád úgy, hogy a Monte Carlo-hiba eloszlását mutassa?

4. Megváltozott feltétel. Változtasd az eloszlást standard normálisra: $X\sim N(0,1)$, $g(x)=x^2$. A pontos várható érték $E[X^2]=1$. Generálj mintákat np.random.default_rng(seed).normal(size=n) hívással; $n=1000$, $R=500$ seeddel. Számold ki a becslések RMSE-ét, egyetlen minta $\widehat{SE}$-jét, és hasonlítsd az elméleti $\sqrt{2/n}$ értékhez. Milyen korlátot kell szem előtt tartani a normális eloszlás farokterülete miatt?

12.3.6 Megoldások és indoklás

1. Az egyenletes $U(0,1)$ sűrűsége $f(x)=1$ a $[0,1]$-en, máshol nulla. Ha $x_1,\dots,x_n$ egyenletesen szórja be a $[0,1]$ intervallumot, akkor a függvényértékek átlaga Riemann-szerűen közelíti a $\int_0^1 g(x)\,dx$ területet. A sűrűség az általános esetben azt mondja meg, hogy mely tartományok kapnak több véletlen pontot; az $E[g(X)]=\int g(x)f(x)\,dx$ éppen ilyen súlyozott átlag. A np.mean(g(x)) akkor becsli ezt jól, ha az $x$ minta valóban $f$-ből származik és a pontok függetlenek.

2. Lehetséges megoldás:

import numpy as np

def mc_mean(g, rng, n):
    x = rng.random(size=n)
    return float(np.mean(g(x)))

g = lambda x: x * x
I_true = 1.0 / 3.0
sizes = [100, 400, 1600, 6400, 25600]
R = 100
rmse = []
for n in sizes:
    est = [mc_mean(g, np.random.default_rng(seed), n)
           for seed in range(R)]
    est = np.array(est)
    rmse.append(np.sqrt(np.mean((est - I_true) ** 2)))
print("rmse:", rmse)
slope = np.polyfit(np.log(sizes), np.log(rmse), 1)[0]
print("log-log meredekség:", slope)

Elvárt meredekség hozzávetőleg $-0{,}5$, mert a standard hiba $1/\sqrt n$ rendű. A tényleges érték $-0{,}6$ és $-0{,}4$ között alakulhat; a meredekség becslésének is van Monte Carlo-hibája. Ellenőrző invariant: a számított RMSE-értékeknek csökkenő sorozatot kell alkotniuk, és a legnagyobb $n$-re az RMSE nagyjából az első érték tizenhatoda közelében kell lennie.

3. A hiba az, hogy minden $n$-hez ugyanazt a seed 7-et használjuk, ezért az ábra egyetlen $n$ szerinti realizációt mutat, nem az eljárás hibaeloszlását. A seed rögzítése reprodukálhatóságot ad, de a $|hiba|$ görbéje ettől még véletlen; más seeddel más képet kaphatunk. A javítás: minden $n$-hez több különböző seeddel ismételjünk, és a hibák RMSE-ét vagy szórását rajzoljuk. Például:

rmse = []
for n in sizes:
    est = []
    for seed in range(200):
        x = np.random.default_rng(seed).random(n)
        est.append(np.mean(np.sin(x)))
    est = np.array(est)
    rmse.append(np.sqrt(np.mean((est - exact) ** 2)))
plt.plot(sizes, rmse)

4. Lehetséges megoldás:

import numpy as np
from scipy import stats

I_true = 1.0
n = 1000
R = 500
x = np.vstack([
    np.random.default_rng(seed).normal(size=n)
    for seed in range(R)
])                            # shape: (R, n)
estimates = np.mean(x ** 2, axis=1)   # shape: (R,)
rmse = np.sqrt(np.mean((estimates - I_true) ** 2))
single_se = np.std(x[0] ** 2, ddof=1) / np.sqrt(n)
print("RMSE:", rmse)
print("egy minta SE:", single_se)
print("elméleti SE:", np.sqrt(2.0 / n))
print("stats kontroll:", stats.norm.expect(lambda z: z ** 2))

Az elméleti szórás: $Var(X^2)=E[X^4]-(E[X^2])^2=3-1=2$, tehát $SE=\sqrt{2/n}\approx0{,}0447$. A futás RMSE-je körülbelül $0{,}04$–$0{,}05$; az egy minta $\widehat{SE}$-je is hasonló. A korlát: a normális farka miatt az $X^2$ eloszlása ferde, így a minta szórása lassabban stabilizálódik, mint szimmetrikus kis szórású esetben; de mivel a variancia véges, a Monte Carlo-becslés továbbra is konzisztens.

12.3.7 Gyakori tévedés és önellenőrzés

A leggyakoribb tévedés: egyetlen seed abszolút hibagörbéjét a módszer hibájának tekintjük. A seed 7-es görbe simának és csökkenőnek tűnhet, de egy másik seed más, akár nem monoton görbét ad. Az egyetlen realizáció nem mutatja a bizonytalanságot. Hasonlóan, a Láncszabály és számítási gráf numerikus gradiens-ellenőrzéséhez hasonlóan a becsült mennyiség bizonytalanságát a módszer bizonytalanságától el kell választani.

Önellenőrző kérdések: - Ha a mintaméretet 25-szörösére növelem, hányad részére csökken a standard hiba? Válasz: körülbelül ötödére, mert $1/\sqrt{25}=1/5$. - A seed megadása garantálja-e, hogy a becslés közelebb kerül a pontos értékhez? Nem; csak reprodukálható folyamot ad. - Mi a különbség $E[g(X)]$ és $g(E[X])$ között? Az első a függvényértékek súlyozott átlaga, a második a várható értéken kiértékelt függvény; általában nem egyenlők, például a 4. feladatnál $E[X^2]=1$, de $(E[X])^2=0$.

12.3.8 Rövid összefoglaló

A várható érték a sűrűségfüggvény szerint súlyozott integrál; a Monte Carlo-becslés ezt véletlen mintából vett átlaggal közelíti. A default_rng(seed) reprodukálható pszeudorandom folyamot ad. A becslés standard hibája a minta szórása osztva $\sqrt n$-nel, így a tipikus hiba $1/\sqrt n$ arányban csökken. Egyetlen seed hibagörbéje csak egy realizáció; a módszer bizonytalanságát több független seedből számolt RMSE vagy becsült SE mutatja. A következő ellenőrző pont ezt a hibakezelést egy veszteségfüggvény teljes értelmezéséhez kapcsolja.