19.2 Credible interval és posterior predictive check

19.2.1 Cél és konkrét helyzet

A Bayesi frissítés rácson alfejezetben rácsos posteriorral becsültük egy Bernoulli-folyamat $p$ siker-valószínűségét. Most nem elégszünk meg a pontbecsléssel: intervallumot és modell-ellenőrzést készítünk.

A konkrét bemenet a 15 elemű megfigyelt sorozat:

$$y_{\mathrm{obs}} = [0,0,1,0,0,0,1,0,0,1,0,1,0,1,1]$$

A megfigyelési egység egyetlen bináris kimenet (0 vagy 1). Az $n=15$ elemű mintában $s=6$ siker és 9 kudarc van. A cél:

1. 95%-os credible interval a $p\in[0,1]$ paraméterre;
2. posterior predictive eloszlás új, azonos szerkezetű adatra;
3. posterior predictive check (PPC), amely az előrejelzett adatot két diagnosztikai statisztikával veti össze a megfigyelttel.

Mielőtt a kódot megnéznéd, adj magadnak gyors becslést. A flat priorral a Bayesi frissítés rácson szerint a posterior Beta$(7,10)$, várható értéke $7/17\approx0.4118$. A 15 új megfigyelés várható sikere tehát $15\cdot 7/17\approx6.18$, a megfigyelt $s=6$ ennek közelében van. A 95%-os intervallum nagyságrendje a posterior Beta-szórása alapján körülbelül $[0.20,0.64]$ lehet; ezt a kódban ellenőrizzük.

19.2.2 Credible interval: a posterior kvantilisei

A bayesi keretben a paraméter valószínűségi változó. A posterior sűrűség $\pi(p\mid y)$ azt fejezi ki, hogy a prior és a likelihood együttes elfogadása után mely $p$-értékek a legplauzibilisebbek. Egy $1-\alpha$ szintű credible interval, magyarul hiteles intervallum, olyan $[L,U]$ tartomány, amelyre

$$\int_L^U \pi(p\mid y)\,dp = 1-\alpha .$$

A leggyakoribb, kvantilisen alapuló változat az equal-tailed interval: a posterior eloszlás $\alpha/2$ és $1-\alpha/2$ kvantilise határolja. Ha a posterior sűrűséget a $g$ rácson diszkréten ismerjük, akkor a kumulált tömegek alapján keressük meg a kvantiliseket:

$$L = q_{\alpha/2}(p\mid y), \quad U = q_{1-\alpha/2}(p\mid y).$$

A Bayesi frissítés rácson rácsos feladatában flat priorral a posterior Beta$(7,10)$. A pontos equal-tailed 95%-os intervallum a SciPy beta.ppf([0.025, 0.975], 7, 10) hívásával kapható meg; az ott megadott ponthoz közeli rácspont pedig a normalizált posteriorból kereshető ki. A két értéknek $1/2001$ nagyságrendű rácshibán belül egyeznie kell.

Az equal-tailed intervallum nem az egyetlen választás. A highest posterior density interval (HPDI) azt a legrövidebb tartományt keresi, amelyre a posterior tömeg $1-\alpha$, ezért ferde eloszlásnál eltérhet az equal-tailed változattól. Itt a Beta$(7,10)$ csak enyhén ferde, ezért a különbség kicsi; a tanulási célhoz az equal-tailed kvantilisintervallum egyszerűbb és közvetlenül kapcsolódik a Standard error és konfidenciaintervallum alfejezetben használt kvantilis fogalmához.

A bayesi intervallum állítása feltételes a priorra és a modellre: adott $y$ mellett a posterior valószínűség $0.95$, hogy $p$ az intervallumban van. Ez nem ugyanaz, mint a konfidenciaintervallum ismételt mintavételes lefedettségi garanciája. A különbséget érdemes azonnal kimondani, még a kód előtt:

A kettő tehát nem automatikusan felcserélhető, még akkor sem, ha egy adott példában a számértékük hasonló.

19.2.3 Predictive distribution és prior predictive

A posteriorból nemcsak a paraméterre, hanem új megfigyelésre is tehetünk valószínűségi állítást. Egy új $\tilde y$ megfigyelés posterior predictive eloszlása a likelihood posterior szerinti keveréke:

$$f(\tilde y \mid y) = \int_0^1 f(\tilde y \mid p)\,\pi(p\mid y)\,dp .$$

A keverés miatt nem elég p-t a posterior átlaggal vagy módusszal helyettesíteni, mert akkor elvész a paraméterbizonytalanság. Egyetlen új Bernoulli-megfigyelésre a posterior átlag jelenik meg:

$$\Pr(\tilde Y=1\mid y) = \int_0^1 p\,\pi(p\mid y)\,dp = \mathbb E\,[p\mid y] = \frac{7}{17}\approx0.4118.$$

Ha viszont $m=15$ új elemű összeget tekintünk, akkor a béta-posteriorral a sikerösszeg eloszlása béta-binomiális:

$$\Pr(\tilde K=k\mid y)=\int_0^1\binom{15}{k}p^k(1-p)^{15-k}\,\mathrm{Beta}(p;7,10)\,dp .$$

Ennek várható értéke $15\cdot7/17=6.18$, de a varianciája nagyobb, mint egy rögzített $p=7/17$ melletti $\mathrm{Binomial}(15,7/17)$ eloszlás varianciája, mert a $p$ bizonytalansága is hozzáadódik.

A prior predictive eloszlás ugyanezt a keveréket az adat megfigyelése előtt, a priorral számítja:

$$f(\tilde y) = \int_0^1 f(\tilde y \mid p)\,\pi(p)\,dp .$$

Ez a prior ellenőrzésére szolgál: mielőtt a likelihood frissítené a modellt, megmutatja, milyen adatokat tart a prior eleve valószínűnek. Flat Beta$(1,1)$ prior esetén a 15 kísérletből álló sikerösszeg prior predictive eloszlása egyenletes:

$$\Pr(\tilde K=k) = \int_0^1 \binom{15}{k}p^k(1-p)^{15-k}\,dp = \frac{1}{16}, \quad k=0,\ldots,15.$$

Ez jól ellenőrizhető kézi invariant: a flat prior nem részesít előnyben semmilyen sikergyakoriságot, ezért a lehetséges 16 összeg mindegyike ugyanolyan plauzibilis.

19.2.4 Posterior predictive check: replikált adat és két diagnosztika

A posterior predictive check (PPC) a modell illeszkedését a következő módon vizsgálja. Minden $r=1,\ldots,S$ replikációban:

1. sorsolunk egy $p^{(r)}$ paramétert a posteriorból;
2. abból generálunk egy $y^{\mathrm{rep},(r)}$ adatsort, amelynek mérete és megfigyelési egysége megegyezik az eredeti $y_{\mathrm{obs}}$ adatéval;
3. kiszámítunk egy-egy $T(y)$ diagnosztikai statisztikát mind a replikált, mind a megfigyelt adatra.

Itt két statisztikát használunk. Az első a sikerösszeg:

$$T_1(y)=\sum_{i=1}^{n} y_i,$$

amely az átlagos szintet ellenőrzi. A második a leghosszabb egyesekből álló futam hossza, $T_2(y)$; ez érzékeny arra, ha a generált sorozatban a sikerek csomósodnak vagy éppen túl szabályosan váltakoznak. A PPC lokális $p$-értéke, a $p_B$, azt méri, hogy a replikált statisztikák mekkora része éri el a megfigyelt értéket:

$$p_B = \frac{1}{S}\sum_{r=1}^{S}\mathbf 1\!\left[T(y^{\mathrm{rep},(r)}) \ge T(y_{\mathrm{obs}})\right].$$

A $p_B \approx 0$ vagy $p_B \approx 1$ azt jelzi, hogy a megfigyelt adat a posterior predictive eloszlás valamelyik szélére esik, tehát a modell az adott statisztikát rosszul írja le. A $p_B$ azonban nem klasszikus nullhipotézis-p-érték: ugyanazt az adatot használja a posterior illesztéséhez és az ellenőrzéshez is, ezért a gyakorisági kalibrációja nem automatikus. Kétoldali diagnosztikaként a $2\min(p_B, 1-p_B)$ érték szolgálhat, de küszöb helyett a modellhibák közötti rangsorolásra érdemes használni.

A calibration szó itt két dolgot jelenthet. Az egyik a posterior predictive intervallumok kalibrációja: ha a paramétert valóban a priornak megfelelően sorsoljuk, az adatot a likelihood szerint generáljuk, majd a posteriorból 95%-os predikciós intervallumot készítünk, a jövőbeli megfigyelés a sok ismétlés során közel 95%-ban esik bele. Ez a joint prior-data folyamatra vonatkozó állítás. Egy rögzített, ismeretlen $p$ mellett viszont egy credible interval nem garantál frekventista lefedettséget; ezt az alábbi gyakori tévedésnél számpéldával is megmutatjuk. A másik jelentés a prediktív valószínűségek pontossága: akkor nevezzük jól kalibráltnak a modellt, ha az 50%-os valószínűséggel jelzett események körülbelül 50%-ban következnek be.

19.2.5 Minimális futtatható munkapélda

A következő, NumPy/SciPy alapú minimal script a rácsos posteriorból kiindulva számítja ki a credible intervalt, generál posterior predictive replikációkat, és két PPC-statisztikát hasonlít össze. A labor környezetében, a projekt gyökeréből futtatható:

# Linux/macOS terminál, a tankönyvi lab projekt gyökerében
uv run python ppc_grid.py

# Windows PowerShell, a projekt gyökerében
uv run python ppc_grid.py

A kód kimenetét itt nem ellenőrzött futtatásként közöljük; a seed megismételhetőséget segít egy gépen belül, de nem garantál bitazonos eredményt minden NumPy-verzión.

import numpy as np
from scipy.stats import beta as beta_dist

rng = np.random.default_rng(20261020)

y_obs = np.array([0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1])
n = len(y_obs)              # 15
s = int(y_obs.sum())        # 6


def longest_run(x: np.ndarray) -> int:
    best = cur = 0
    for v in x:
        if v == 1:
            cur += 1
            best = max(best, cur)
        else:
            cur = 0
    return best


# Rácsos posterior: flat Beta(1,1) prior
grid = np.linspace(0.0, 1.0, 2001)
prior = np.ones_like(grid)
lik = grid**s * (1.0 - grid) ** (n - s)
post = prior * lik
post /= post.sum()                       # normalizálás: sum(post) == 1

# 95%-os equal-tailed credible interval a rácsból
cdf = np.cumsum(post)
lo = grid[int(np.searchsorted(cdf, 0.025))]
hi = grid[int(np.searchsorted(cdf, 0.975))]

# Posterior predictive replikációk
S = 5000
p_draws = rng.choice(grid, size=S, p=post)   # shape: (S,)
y_rep = rng.binomial(1, p_draws[:, None], size=(S, n))  # shape: (S, 15)

sum_rep = y_rep.sum(axis=1)                              # shape: (S,)
run_rep = np.array([longest_run(row) for row in y_rep])  # shape: (S,)

t1_obs = s
t2_obs = longest_run(y_obs)
p_b_sum_right = float(np.mean(sum_rep >= t1_obs))
p_b_run_right = float(np.mean(run_rep >= t2_obs))
p_b_sum_two = float(2.0 * min(p_b_sum_right, 1.0 - p_b_sum_right))
p_b_run_two = float(2.0 * min(p_b_run_right, 1.0 - p_b_run_right))

print(f"posterior mean: {np.sum(grid * post):.4f} (elvárás: {7/17:.4f})")
print(f"grid CI:        [{lo:.3f}, {hi:.3f}]")
print(f"Beta(7,10) CI:  {beta_dist.ppf([0.025, 0.975], 7, 10).round(3)}")
print(f"pred. sum mean: {sum_rep.mean():.2f} (elvárás: {n*7/17:.2f})")
print(f"observed: sum={t1_obs}, longest_run={t2_obs}")
print(f"PPC sum: right={p_b_sum_right:.3f}, two-sided={p_b_sum_two:.3f}")
print(f"PPC run: right={p_b_run_right:.3f}, two-sided={p_b_run_two:.3f}")

A kód minden replikációban a posteriorból sorsol p-t, majd abból generálja a 15 Bernoulli-megfigyelést. A replikált sorok száma $S$, az oszlopok száma $n$, tehát a replikáció nem változtatja meg az eredeti megfigyelési egységet. A statisztikák közül sum_rep a sikerösszeg, run_rep a leghosszabb sikeres futam.

19.2.6 Eredmény értelmezése és kézi kontroll

A kód futtatása nélkül is rögzíthetők az elméletileg várt ellenőrző pontok:

• A normalizált rácstömeg összege pontosan $1$ a post /= post.sum() lépés után; a float kerekítés miatt $10^{-16}$ nagyságrendű eltérés elfogadható.
• A posterior átlag np.sum(grid * post) a Beta$(7,10)$ várható értékéhez, $7/17=0.4118$-hoz tart. A rács és a folytonos béta közötti numerikus eltérés tipikusan $10^{-4}$ alatt marad.
• A posterior predictive sikerösszeg átlaga a béta-binomiális várható érték $15\cdot7/17\approx6.18$ körül várható. A Monte Carlo standard error a predikciós szórással ($\approx2.5$) és $S=5000$-rel számolva körülbelül $0.04$, ezért a $\pm0.05$ tűrés elfogadható a futás ellenőrzéséhez.
• A gridből számolt credible intervalnak a beta.ppf intervallumhoz közelinek kell lennie; a rácshiba a 2001 pontból adódóan kisebb $0.001$-nél. A megfigyelt $T_1=6$ a prediktív átlag alatt van, de a prediktív szórás tartományában.

A PPC-számok értelmezésénél ne egyetlen gépi küszöböt nézzünk. Ha például a p_b_sum_two 0.2 és 0.8 között van, a sikerösszeg nem mutat durva ütközést a modellel. Ha viszont a leghosszabb futamra a kétoldali érték 0.05 alá esik, akkor a modell nem jól ragadja meg a megfigyelt sorozat szomszédsági mintázatát; ilyenkor érdemes az időbeli függőséget vagy az egyenlő valószínűség feltételét gyanúba fogni. A p_B tehát diagnosztika, nem a parameter igazságáról szóló döntés.

Intervallumértelmezés: a 95%-os credible interval a modell és a prior feltételei mellett 95% posterior tömeget tartalmaz. Ettől különbözik a gyakorisági konfidenciaintervallum lefedettségi állítása. A posterior predictive check modellspecifikációs diagnosztika, nem az általános modelligazság igazolása.

19.2.7 Feladatok

1. Értelmezés. A futás után a 95%-os equal-tailed credible interval körülbelül $[0.20,0.64]$. Kollégád ezt mondja: „Ebből következik, hogy a valódi $p$ 95% valószínűséggel 0.20 és 0.64 között van, tehát ha 100 hasonló mintát gyűjtenénk, akkor 95 ilyen intervallum tartalmazná a valódi $p$-t.” Melyik mondatrész igaz, melyik hamis? Pontosan nevezd meg, milyen feltétel mellett igaz az első rész, és miért nem következik belőle a második.

2. Megvalósítás. Írj egy quantile_interval(post, grid, level) függvényt, amely adott level mellett equal-tailed credible intervalt ad vissza. Futtasd 80%-os szinten, és növeld a posterior predictive replikációk számát $S=10{,}000$-re. Ellenőrizd a p_draws, y_rep, sum_rep alakját, és add meg mindkét PPC-statisztika jobb oldali és kétoldali értékét. Az input a fenti y_obs és a flat prior.

3. Hibakeresés. Az alábbi kódrészlet a posterior predictive eloszlást a posterior móduszával próbálja generálni:

p_map = grid[np.argmax(post)]
y_wrong = rng.binomial(1, p=p_map, size=(S, n))
wrong_mean = y_wrong.sum(axis=1).mean()
print(wrong_mean)

Mi a statisztikai hiba? Melyik varianciatag hiányzik a generált eloszlásból? Milyen következménye lehet ennek a PPC-re? Javítsd úgy a sort, hogy a posterior bizonytalanságát is megőrizze.

4. Transzfer megváltozott priorra. Ismételd meg a teljes számítást a szkeptikus Beta$(1,4)$ priorral, amelynek sűrűsége arányos $(1-p)^3$-nel. Mi lesz a posterior család és a posterior predictive sikerösszeg várható értéke? Hasonlítsd össze a flat és a szkeptikus prior predictive eloszlást is a 15 kísérletből álló, még nem látott mintára. Az eredmény hogyan változtatja meg a megfigyelt $T_1=6$ értelmezését?

19.2.8 Megoldások és indoklás

1. A mondat első fele a bayesi Markov-feltételekkel pontosítható: ha elfogadjuk a Beta$(1,1)$ priort és a Bernoulli-likelihoodot, akkor adott $y$ mellett a posterior valószínűség $0.95$, hogy $p\in[L,U]$. Ez posterior valószínűségi állítás a paraméterről. A második fele hamis: abból, hogy a posterior tömeg 95%-ot ad az intervallumnak, nem következik, hogy az eljárás rögzített, ismeretlen $p$ mellett ismételt mintavételnél 95%-os lefedettséget ad. A konfidenciaintervallum lefedettsége az eljárás ismételt tulajdonsága, és nem függ a priortól; a credible interval ezt nem garantálja. A két fogalom tehát más kérdésre válaszol.

2. Lehetséges megoldás:

def quantile_interval(post, grid, level):
    lo_q = (1.0 - level) / 2.0
    hi_q = 1.0 - lo_q
    cdf = np.cumsum(post)
    lo = grid[int(np.searchsorted(cdf, lo_q))]
    hi = grid[int(np.searchsorted(cdf, hi_q))]
    return lo, hi

lo80, hi80 = quantile_interval(post, grid, 0.80)
print(f"80% grid CI: [{lo80:.3f}, {hi80:.3f}]")
print("80% Beta CI:", beta_dist.ppf([0.10, 0.90], 7, 10).round(3))

S = 10000
p_draws = rng.choice(grid, size=S, p=post)
y_rep = rng.binomial(1, p_draws[:, None], size=(S, n))
sum_rep = y_rep.sum(axis=1)

assert p_draws.shape == (S,)
assert y_rep.shape == (S, n)
assert sum_rep.shape == (S,)

p_b_right = float(np.mean(sum_rep >= s))
p_b_two = 2.0 * min(p_b_right, 1.0 - p_b_right)
print(f"PPC sum: right={p_b_right:.3f}, two-sided={p_b_two:.3f}")

A 80%-os intervallum szűkebb, mint a 95%-os, mert a posterior eloszlás középső tömegét fogja közre; a megfigyelt adat nem változott. A replikációk alakja a befoglaló mátrix dimenzióival ellenőrizhető.

3. A p_map a posterior módusza, azaz egyetlen legvalószínűbb pont. Ha minden replikációt ezzel az egy értékkel generálunk, akkor valójában $\mathrm{Binomial}(15,p_{\mathrm{map}})$ eloszlást kapunk, nem a posterior predictive-eloszlást. A béta-binomiális varianciája két tagra bomlik:

$$\operatorname{Var}(K) = \mathbb E[nP(1-P)] + \operatorname{Var}(nP),$$

ahol az első tag a mintavételi zaj, a második a paraméter posterior bizonytalansága. A rögzített $p_{\mathrm{map}}$ mellett a második tag nulla, ezért a szórás alulbecsült lesz. A PPC-ben ez mesterségesen szűk prediktív eloszlást ad, és a megfigyelt statisztikát indokolatlanul szélsőségesnek mutathatja. A javítás az, hogy p-t a teljes posteriorból sorsoljuk:

p_draws = rng.choice(grid, size=S, p=post)
y_rep = rng.binomial(1, p_draws[:, None], size=(S, n))

Az így kapott eloszlás már a posterior predictive-et közelíti.

4. A szkeptikus prior sűrűsége $\pi(p)\propto(1-p)^3$, ami Beta$(1,4)$. A likelihood továbbra is $p^6(1-p)^9$. A posterior tehát arányos $p^6(1-p)^{12}$-nel, ami Beta$(7,13)$. Várható értéke:

$$\mathbb E\,[p\mid y] = \frac{7}{7+13} = 0.35.$$

A 15 jövőbeli elemű posterior predictive sikerösszeg várható értéke $15\cdot0.35=5.25$. A prior predictive eloszlás flat priorral az egyenletes $k=0,\ldots,15$ eloszlás, várható értéke $7.5$. A szkeptikus prior predictive Beta-binomiális$(15,1,4)$, várható értéke $15\cdot1/(1+4)=3$. A prior tehát alapból kevés sikert tart plauzibilisnek; a megfigyelt $T_1=6$ ebben a keretben már nem a prediktív eloszlás alsó oldala, hanem a szkeptikus prior utáni posteriorhoz képest kissé magasabb, kevésbé szélsőséges érték. A kódban ezt például a scipy.stats.betabinom.pmf segítségével ellenőrizhetjük:

from scipy.stats import betabinom

k_vals = np.arange(n + 1)
prior_pred_flat = betabinom.pmf(k_vals, n, 1, 1)      # közel egyenletes
prior_pred_skept = betabinom.pmf(k_vals, n, 1, 4)
post_pred_skept = betabinom.pmf(k_vals, n, 7, 13)

print("prior pred mean flat:  ", float(np.sum(k_vals * prior_pred_flat)))
print("prior pred mean skept: ", float(np.sum(k_vals * prior_pred_skept)))
print("post  pred mean skept: ", float(np.sum(k_vals * post_pred_skept)))

A kézi kontroll szerint 7.5, 3.0 és 5.25 körüli értékeket várunk.

19.2.9 Gyakori tévedés ellenpéldával

A leggyakoribb hiba az, hogy a credible intervalt automatikusan konfidenciaintervallumként értelmezzük. Ellenpélda: tegyük fel, hogy a prior degenerált, azaz teljes tömegét egyetlen $p=0.5$ pontra teszi. Ekkor bármilyen adat után a 95%-os credible interval is a $[0.5,0.5]$ pont. Ha a valódi $p=0.3$, akkor ez az intervallum ráadásul soha nem tartalmazza a paramétert, tehát az ismételt mintavételes lefedettsége 0. A posterior állítás a modellen belül továbbra is „igaz”, de a prior hibás volt a külvilághoz képest. Ezért a credible interval nem állít frekventista lefedettséget, és a konfidenciaintervallum sem állít posterior valószínűséget.

A második hiba a PPC-p-érték mechanikus 0.05-ös küszöbölése. A $p_B$ nem nullmodell szerinti p-érték, és ugyanaz az adat kétszeresen szerepel: először a posterior illesztésében, másodszor az ellenőrzésben. Ezért a $p_B$ eloszlása akkor sem feltétlenül egyenletes, ha a modell igaz. A gyakorlatban a nagyon kicsi vagy nagyon nagy érték irányíthatja a figyelmet a modellhibára, de döntési küszöbként nem ugyanolyan státuszú, mint a Standard error és konfidenciaintervallum kontrollált ismételt mintavételes eljárása.

19.2.10 Önellenőrzés

• Le tudom-e írni, hogy a 95%-os credible interval melyik feltétel mellett állít posterior valószínűséget 0.95-öt?
• A replikált adatomnak ugyanaz a megfigyelési egysége és sor-/oszlopszáma, mint az eredeti adatnak, tehát nem keverem a szemcsézettséget?
• A predictive eloszlásom a posterior teljes keveréke, nem pedig egyetlen pontbecslés körüli binomiális eloszlás?
• A kalibrációs állításomnál meg tudom nevezni, hogy a véletlent a prior-data közös folyamatból vettem, nem rögzített paraméterből?

19.2.11 Rövid összefoglaló

A credible interval a posterior kvantiliseiből készített tartomány: flat priorral a Beta$(7,10)$ posterior 95%-os intervalluma körülbelül $[0.20,0.64]$. A posterior predictive eloszlás a likelihood posterior szerinti keveréke; 15 új Bernoulli-megfigyelés sikerösszegének várható értéke $15\cdot7/17\approx6.18$. A prior predictive a prior ellenőrzésére való, flat priorral egyenletes sikerösszeget ad. A PPC replikált adatokat generál a posteriorból, és két statisztikát, itt a sikerösszeget és a leghosszabb futást veti össze. A $p_B$ diagnosztikai érték, nem klasszikus p-érték. A credible és confidence interval értelmezése különbözik. Következő alkalmazásként a PyMC modell és MCMC diagnosztika alfejezetben ugyanezt PyMC és MCMC segítségével ismételjük meg.