19.3 PyMC modell és MCMC diagnosztika

19.3.1 Helyzet és cél

Folytatjuk a Credible interval és posterior predictive check rácsos bayesi példáját. Egy mobilalkalmazás kísérleti csoportjában 16 véletlenül kiválasztott ügyfelet figyelünk meg; minden ügyfél egyetlen megfigyelési egység. A bináris kimenet $y_i \in \{0,1\}$ azt jelöli, hogy az $i$-edik ügyfél 7 napon belül konvertált-e. A megfigyelt adatban 6 siker és 10 kudarc van.

Célunk a konverziós valószínűség $p$ paraméterének posterior eloszlását MCMC-mintavétellel becsülni, majd csak akkor értelmezni, ha a konvergencia-diagnosztika rendben van. A kézi kontrollt már ismerjük: $p \sim \mathrm{Beta}(1,1)$ flat prior és 6/10 siker/kudarc esetén a konjugált posterior $p \mid y \sim \mathrm{Beta}(7,10)$. Ennek várható értéke $\frac{7}{17} \approx 0.412$, a Credible interval és posterior predictive check alfejezetben számított 95%-os credible interval körülbelül $[0.20, 0.64]$. A PyMC-futtatásnak ehhez kell közelítenie a Monte Carlo hiba mértékén belül.

Az eszköz itt bővítő: NumPy, SciPy mellett PyMC és ArviZ. Az alábbi kód CPU-n fut; GPU-ra vagy távoli gépre nincs szükség.

19.3.2 Miért lánc, és mit jelent a mintavétel?

A posterior eloszlásból közvetlenül gyakran nem tudunk független mintákat húzni. Az MCMC (Markov chain Monte Carlo) ehelyett egy Markov-láncot futtat: a következő $p_{t+1}$ paraméterérték csak a jelenlegi $p_t$ értéktől függ. Elég hosszú futtatás után a lánc állapotainak eloszlása a posteriorhoz közelít; a bemelegedés (tuning/warmup) után megtartott lépések pedig közelítő poszteriori mintaként használhatók.

A lépések nem függetlenek. A Bootstrap a megfelelő mintavételi egységgel bootstrapjében a független újramintázási egységért aggódtunk; itt éppen azt mérjük, hogy a lánc előző lépése milyen erősen befolyásolja a következőt. Emiatt a nyers draws szám nem azonos az információtartalommal.

A PyMC folytonos paraméterű modelleknél a NUTS (No-U-Turn Sampler) nevű Hamilton-féle Monte Carlo módszert használja. A lánc egy lépése:

1. javasol egy új paraméterértéket a posterior geometriája alapján,
2. elfogadja vagy elutasítja azt,
3. az új állapotot rögzíti,
4. a következő iteráció ebből az állapotból indul.

Több láncot (chain) különböző kezdőértékből indítunk. Így a láncok közötti egyezés vagy eltérés maga is diagnosztikai jel. A kontrollpélda még egyszerű: a konjugált posterior miatt a láncok futtatása „önellenőrző” futás, mert az elméleti eloszlást külön is ismerjük.

19.3.3 A diagnosztikai fogalmak kis kézi kontrollal

Az ArviZ InferenceData objektumba gyűjti a láncok értékeit. A fő jellemzők:

• draws: a bemelegedő fázis után, lánconként megtartott minták száma. Ha chains=4 és draws=1000, akkor összesen 4000 poszteriori mintánk van.
• tune: a bemelegedő lépések száma. Az itt keletkező mintát nem használjuk a posterior értelmezéséhez.
• R-hat (Gelman–Rubin statisztika): a láncok közötti és a láncon belüli szóródást hasonlítja össze. Az 1,01 fölötti $\hat R$ diagnosztikai figyelmeztetés; a küszöb nem bizonyítja önmagában a konvergenciát vagy annak hiányát. A diagnosztikát ESS-szel, divergenciákkal és trace-plotokkal együtt kell értelmezni.
• ESS (effective sample size, effektív mintaméret): a lánc mintáinak autokorrelációját korrigáló szám. A ess_bulk a poszterior középső részére, az ess_tail a farokra jellemző. A cél az, hogy legalább néhány száz legyen; a példánkban több ezer várható.
• Mcse: Monte Carlo standard error; azt mondja meg, a finite lánc miatt mekkora a becsült posterior-átlag numerikus bizonytalansága.
• divergences: a NUTS futása közben a lépés geometriája elfajulhat. A divergens lépések darabszáma legyen nulla, mert hullámzó lefolyás esetén a lánc könnyen félrevezető területre sodródhat.

Kis kézi ellenőrzés a fogalomra. Képzeljünk el két láncot, mindegyik 2000 poszteriori mintával. Az A lánc minden értéke 0.42 körül, a B lánc minden értéke 0.66 körül mozog. A láncokon belül kicsi a szórás, a láncok között nagy. A közös posterior-átlag 0.54 körül lenne, de ez nem a kívánt $p \approx 0.41$ poszterior összegzése, csak két megszorult lánc átlaga. Az R-hat ilyenkor jóval 1.01 fölött lenne: a sok draw ellenére rossz a keveredés. A konvergencia tehát nem a mintaszámból, hanem a láncok közötti kiegyenlítettségből és az autokorrelációból olvasható le.

19.3.4 Teljes minimális munkapélda

Először opcionális PyMC/ArviZ környezetet hozunk létre. A pontos, elérhető verziót a saját futtatási napló rögzíti; a célterjedelem a PyMC 5.x család. Az alábbi parancscsomag szemléltető, nem ellenőrzött lock.

# macOS/Linux; munkakönyvtár: lab/
python -m venv .venv-s075
source .venv-s075/bin/activate
python -m pip install "pymc>=5,<6" "arviz>=0.19"

Windows PowerShellben ugyanez a környezet:

# PowerShell; munkakönyvtár: lab\
python -m venv .venv-s075
.\.venv-s075\Scripts\Activate.ps1
python -m pip install "pymc>=5,<6" "arviz>=0.19"

A munkapélda a $p$ paramétert és a megfigyelt Bernoulli-adatokat definiálja. A with pm.Model() blokk a PyMC modellkörnyezete: benne a véletlen változók függőségi gráfja épül fel.

# minimal_mcmc.py
import numpy as np
import pymc as pm
import arviz as az

# 16 ügyfél, 6 siker, 10 kudarc
y = np.array([1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

with pm.Model() as model:
    p = pm.Beta("p", alpha=1.0, beta=1.0)      # flat prior a [0,1] intervallumon
    pm.Bernoulli("y_obs", p=p, observed=y)     # likelihood: yi ~ Bernoulli(p)

    idata = pm.sample(
        draws=1000,       # lánconkénti poszteriori mintaszám
        tune=1000,        # bemelegedő fázis lánconként
        chains=4,         # független Markov-láncok száma
        cores=1,          # CPU-n futtatva egyértelműbb naplózni
        random_seed=20261020,
        return_inferencedata=True,
    )

# Diagnosztika és poszteriori összegzés
print(az.summary(idata, var_names=["p"], hdi_prob=0.95))
print("R-hat:", az.rhat(idata, var_names=["p"]).to_array().item())
print("ESS bulk:", az.ess(idata, var_names=["p"], method="bulk").to_array().item())
div = idata.sample_stats["diverging"].sum().item()
print("Divergences:", int(div))

A p változó shape-je a posteriorban (4, 1000):

• 0. tengely: lánc azonosítója, azaz 4 független futtatás
• 1. tengely: láncon belüli sorszám, azaz 1000 megtartott lépés

A pm.Bernoulli("y_obs", ...) likelihood az egyes megfigyeléseket egy Bernoulli- véletlen változón keresztül köti a p paraméterhez. A PyMC a modell log-posteriorját építi fel ebből a gráfból, és a NUTS ezen a felületen javasol lépéseket.

Az elméletileg várt, szemléltető kimenet:

# posterior summary (szemléltető, nem helyi mért output)
        mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat
p      0.412  0.116   0.199    0.635      0.002    0.001      2200      2100    1.001
R-hat: 1.001
ESS bulk: 2200
Divergences: 0

A számok az elméleti Beta(7,10) poszterior környékén mozognak. A pontos eredmények futásonként némileg eltérhetnek, mert a MLMC nem bitazonos az eszközök között; az ellenőrzési küszöbök viszont stabilak.

19.3.5 Az eredmények kézi értelmezése

A posterior-összegzés a poszterior közelítő átlagát és szórását mutatja. A 0.412 közel van a konjugált várható értékhez: $\frac{7}{17} \approx 0.412$. A láncok keveredése megfelelő, mert

• R-hat $\approx 1.001 \ll 1.01$,
• a láncok külön-külön is ugyanazt a posterior módusz körüli tartományt fedték le,
• divergens lépés nincs,
• az ESS bulk 2200 körüli, tehát a négyzetes MC hiba az átlagban hasonló ahhoz, mintha, óvatosan fogalmazva, körülbelül 2200 kevésbé autokorrelált mintánk lenne. A pontos független ekvivalencia nem állítható; az ESS mérése maga is statisztikai becslés. Ez a tartomány elég a napi tantételhez.

A 95%-os HDI a [0.199, 0.635] körüli intervallum, ami a Credible interval és posterior predictive check kézi kontrolljával összhangban van. Az eltérés az MC-véletlen miatt elfogadható. Ha ehelyett a posterior-átlag például 0.46 vagy 0.36 lenne, és az R-hat 1.05 fölé menne, akkor a konvergenciát nem vehetnénk készpénznek; a modell építését kellene felülvizsgálni.

A „divergences” értelmezése: a divergens lépés nem maga a HDI-intervallum divergenciája, hanem a HMC-lépés numerikus hibájának jelzése. A nulla jó jel. Ha néhány százalék divergens, a posterior farokbecslései, például a felső 97.5%-os kvantilis, torzultak lehetnek, még ha az átlag stabilnak is tűnik.

A modellhez most nincs szükség a Megfigyeléses becslés és érzékenység megfigyeléses elemzésében megismert confounding-kérdésre, mert ez a 16 ügyfél egyetlen randomizált kísérleti csoportjának tekinthető. Az MCMC önmagában nem igazítja meg a nem mért zavaró változóból eredő torzítást: a Megfigyeléses becslés és érzékenység érzékenységi logikája a bayesi modellnél is érvényes. A lánc lépései pedig nem hasonlítanak a Bootstrap a megfelelő mintavételi egységgel ügyfélszintű bootstrapjére: ott ügyfélrekordokat újramintázunk, itt paraméterértékeket mozgatunk a láncon. Az időben rendezett adathoz pedig a CTE, window és előzményből számolt feature ablak-gondolkodása hasonlít: a lánc állapota az előző állapotból indul, ezért az előzményfüggőséget mérni kell, nem egyszerűen a sorok számát számolni.

19.3.6 Feladatok

19.3.6.1 1. feladat — Értelmezés

Egy munkatárs a következő PyMC kimenetet adja át egy másik konverziós adathoz:

mean=0.414, sd=0.116, hdi=(0.199,0.635), ess_bulk=1200, r_hat=1.002, divergences=0

Írj mellé két–három mondatos értelmezést, amelyben

• megadod, mit mond a posterior a konverzióról,
• indoklod, hogy a lánc miért alkalmas a következtetésre,
• megnevezed, milyen korlátot mégis figyelembe kell venni.

19.3.6.2 2. feladat — Implementálás

Írj egy sample_conversion(successes, failures, prior_a=1, prior_b=1) függvényt, amely a fenti mintához hasonlóan Bernoulli-modellt illeszt. A függvény a successes és failures egészekből állítsa elő az y tömböt, futtasson 2 láncot draws=500, tune=500 beállítással, és adja vissza az idata objektumot.

Ellenőrzés: 3 siker és 4 kudarc esetén, flat prior mellett a konjugált posterior $p \sim \mathrm{Beta}(4,5)$, várható értéke $\frac{4}{9} \approx 0.444$. Futtasd a függvényt, és hasonlítsd össze az az.summary(idata) átlagát az elméleti 0.444-gyel. Adj meg egy toleranciát is, és indokold, miért nem várhatunk hajszálpontos egyezést.

19.3.6.3 3. feladat — Hibadiagnózis

Egy másik munkatárs így futtatja a modellt:

pm.sample(draws=100, tune=0, chains=2, cores=1)

Az az.summary(idata) kimenetében r_hat=1.6, ess_bulk=8. A posterior-átlag 0.31, pedig az adat szerint a nyers sikerarány 0.42. Sorold fel a hibalehetőségeket, és adj két konkrét javítást a munkatársnak. Melyik diagnosztikai értékből látja azonnal, hogy a posterior-átlagot nem szabad a 0.31-gyel összefoglalni?

19.3.6.4 4. feladat — Megváltozott feltétel

Az üzleti csapat információja szerint a korábbi kísérletben a konverzió már 0.3 körül volt, ezért most a flat prior helyett gyengén informatív $\mathrm{Beta}(2,5)$ priort használunk, ahol a prior várható értéke $\frac{2}{7} \approx 0.286$. Új adatként 14 siker és 3 kudarc érkezik. Mi a pontos konjugált posterior alakja, várható értéke és a prior–posterior elmozdulás iránya? Ha a PyMC futásban r_hat=1.000, ess_bulk=2500, divergences=0, elegendő-e a lánc a következtetéshez? Indokold a küszöbértékekkel.

19.3.7 Megoldások és indoklás

19.3.7.1 1. feladat megoldása

A posterior-átlag szerint a konverziós valószínűség körülbelül 0.41, a 95%-os HDI pedig 0.20 és 0.64 közé esik, tehát a becslés jelentős bizonytalanságot hordoz akkor is, ha a lánc jól keveredik. A lánc diagnosztikája jó, mert R-hat $\approx 1.002$ messze az 1.01-es határ alatt van, divergens lépés nincs, és az ESS bulk 1200 számottevő. Korlát: a modell egy egyszerű Bernoulli-likelihoodot használ, tehát az eredmény nem feltétlenül ügyfélszegmensenként is igaz, és a posterior csak a választott prior és likelihood mellett érvényes.

19.3.7.2 2. feladat megoldása

import numpy as np
import pymc as pm

def sample_conversion(successes, failures, prior_a=1, prior_b=1):
    y = np.array([1] * successes + [0] * failures, dtype=int)
    with pm.Model() as model:
        p = pm.Beta("p", alpha=prior_a, beta=prior_b)
        pm.Bernoulli("y_obs", p=p, observed=y)
        idata = pm.sample(
            draws=500,
            tune=500,
            chains=2,
            cores=1,
            random_seed=42,
            return_inferencedata=True,
        )
    return idata

idata = sample_conversion(3, 4)
print(az.summary(idata, var_names=["p"]))

A konjugált várható érték 0.444. A megtartott minták 2 lánc $\times$ 500 minta, tehát az MC hiba nagyságrendje a poszterior szórásnak körülbelül az ESS négyzetgyökével osztva adódik. Ha ESS bulk tipikusan néhány száz, a 0.01–0.02 körüli átlagbeli eltérés elfogadható. Ezért bemondhatunk egy $\pm 0.02$ abszolút toleranciát a posteriior-átlag ellenőrzésére, és egy 0.01-nél kisebb R-hat küszöböt. Ha a becsült átlag 0.41 vagy 0.46, az MCMC-véletlen; de ha 0.30 és R-hat magas, akkor a lánc nem konvergált.

Pontosan: ellenőrzéskor a lánconkénti átlagoknak is közel kell esniük egymáshoz; az ArviZ idata.posterior.mean(dim="draw") ábrája mutathatja a kiegyensúlyozottságot.

19.3.7.3 3. feladat megoldása

A hiba négy forrásra vezethető vissza:

1. tune=0: a lánc kezdőállapotai nem estek a posterior tömegébe, így az elfogadási fázis torzult maradt.
2. draws=100, chains=2: nagyon kevés mintavétel, különösen autokorrelált láncban.
3. A r_hat=1.6 egyenesen jelzi, hogy a két lánc különböző tartományokat talált meg; a közös átlag értelmetlen, mert a láncok közötti variancia nem tűnt el a futam alatt.
4. Az ess_bulk=8 pedig azt mutatja, hogy 200 nyers mintából csupán körülbelül 8 effektív lépés maradt; ez a 0.31 magyarázata is lehet, ha az egyik lánc megszorult egy dudorban.

Konkrét javítások:

• állítsa be tune=1000, draws=2000, chains=4 értékeket;
• adjon meg target_accept=0.9 értéket, hogy a NUTS óvatosabb lépésekkel dolgozzon;
• ezután ismét ellenőrizze az R-hat, ESS és divergenc metrikákat.

Ha a kis mintás paraméterezés valóban nehezen keveredik, a modell újraparaméterezése vagy a prior szűkítése is segíthet.

19.3.7.4 4. feladat megoldása

A konjugált prior–likelihood frissítés:

• prior: $\mathrm{Beta}(2,5)$
• megfigyelés: 14 siker, 3 kudarc
• posterior: $p \mid y \sim \mathrm{Beta}(2 + 14,\; 5 + 3) = \mathrm{Beta}(16, 8)$

A posterior várható értéke:

$$E[p \mid y] = \frac{16}{24} = \frac{2}{3} \approx 0.667.$$

A posterior módusza $\frac{16-1}{24-2} = \frac{15}{22} \approx 0.68$, ami eligazodási pont az egycsúcsú posteriorhoz. A prior 0.286-ról a posterior 0.667-re mozdul el, mert a megfigyelt 14/17 sikerarány jelentősen meghaladja a korábbi előzetes várakozást.

A diagnosztika elegendő: r_hat=1.000 és divergences=0 azt jelzi, hogy a láncok keveredtek, ess_bulk=2500 pedig a poszterior középső tartományában komoly effektív mintát jelent. A gyengén informatív prior nem akadályozta a lánc mozgását; a posterior mégis a 0.67 köré koncentrál, mert a likelihood messze tolta a frissítést.

19.3.8 Gyakori tévedés és ellenpélda

Sokan azt hiszik, hogy a nagy draws érték automatikusan elegendő. Ellenpélda: futtassunk két láncot, egyenként 50 000 poszteriori mintával. Az A lánc minden lépésben 0.42 körül marad, a B lánc 0.66 körül. A raktárban a posterior.mean() például 0.54 lenne, és az ess_bulk sem feltétlenül lenne nulla. Az R-hat azonban magas: hiába 100 000 draw, a láncok közötti variancia nem csökkent, mert a Markov-lánc nem vándorolt a teljes poszterioron. A jó diagnosztika ezért nem a minta számát, hanem a láncok közötti egyezést és az autokorrelációt ellenőrzi.

19.3.9 Önellenőrzés

1. Mi a különbség a posterior valódi eloszlása és az MCMC-közelítés között?
2. Milyen küszöb alatt fogadunk el egy R-hat értéket a napi munkában?
3. Miért nem egyenlő az ESS a nyers minta méretével?
4. Hány független lánccsal dolgozik a minimális példa, és ez miért jobb, mint egyetlen lánc 4000 mintával?
5. Ha 20 divergens lépés jelenik meg 4000-ből, mi a helyes reakció, mielőtt a posterior felső kvantilisét jelentenénk?

A válaszok a szövegben megtalálhatók: a posterior a cél, a láncok közelítik; R-hat jellemzően 1.01 alatt; az ESS az autokorreláció miatt csökken; 4 lánccal diagnosztizálható a nem keveredés; divergens lépés esetén a farok mennyiségeit nem szabad lejelenteni, hanem a mintavételezőt kell javítani.

MCMC-diagnosztika: pusztán a mintaszám és a láncok lefutása nem bizonyít konvergenciát. Ellenőrizendő a rank-normalized R-hat, bulk/tail effektív mintanagyság, divergenciák, energiamutatók és a trace plot; hibás diagnosztikával posterior intervallum nem tekinthető megbízhatónak.

19.3.10 Rövid összefoglaló

A PyMC-modellben a p paramétert priorral és Bernoulli-likelihooddal kötjük a megfigyelt adathoz, majd a NUTS-láncok a poszterior közelítő mintáját adják. A 4 lánc és a bemelegedő fázis után az ArviZ r_hat, ess_bulk és divergenciaszámlálója ellenőrzi a keveredést. Az R-hat 1.01 alatt, a divergenciaszám nulla, az effektív minta legalább több száz legyen. Sok nyers draw nem pótolja a jó keveredést. A PyMC modell és MCMC diagnosztika következő alkalmazása a Hierarchia és részleges pooling alfejezetben hierarchikus modellek és részleges pooling, ahol ugyanezen diagnosztikát csoportonkénti paraméterekre kell kiterjeszteni.