15.3 Bootstrap a megfelelő mintavételi egységgel

15.3.1 Helyzet és előzetes becslés

A Northstar orders.csv 2400 rendelési sora 300 ügyfélhez tartozik. Az üzemeltetési csapat azt kérdezi: mennyire bizonytalan a delivery_days átlaga és mediánja? A rendelési sor kézenfekvő megfigyelési egység, de ugyanannak az ügyfélnek több rendelése van, tehát a sorok nem függetlenek — ezt már a Függetlenség, kovariancia és függő adat alfejezetben láttuk. A kérdés tehát nemcsak az, hogy mekkora a pontbecslés, hanem az is, hogy melyik szinten érdemes újramintázni.

Előzetes becslés: ha a szállítási napok szórása kb. 3 nap, akkor a függetlennek vélt sorokra számolt standard error $s/\sqrt{2400} \approx 0{,}06$ nap lenne. Ha viszont az ügyfelen belüli szállítási idők hasonlóak, a hatékony mintaméret kisebb, és a valódi standard error akár 0,1–0,2 nap is lehet. A feladatunk ezt bootstrap segítségével ellenőrizhetővé tenni.

15.3.2 Fogalmi magyarázat: percentile és BCa

A bootstrap (újramintázásos becslés) a megfigyelt mintát tekinti a sokaság legjobb helyettesítőjének. Ha az eredeti adatvektor $x=(x_1,\dots,x_n)$, akkor $B$-szer húzunk $n$ elemet visszatevéssel, és minden $b$ húzásra kiszámítjuk a $\hat\theta^*_b = g(x^*_b)$ statisztikát. A $\hat\theta^*$ replikációk vektora $(B,)$ alakú, és a cél nem új adat előállítása, hanem a $\hat\theta$ becslő mintavételi eloszlásának közelítése. Ez a megközelítés akkor működik jól, ha a minta reprezentatív a kérdéses sokaságra nézve, és az $n$ elem cserélhető, tehát a sorrend nem hordoz releváns csoportszerkezetet.

A legegyszerűbb, percentilis (percentile) intervallum a replikációk eloszlásából számított kvantilis:

$$\mathrm{CI}_{\mathrm{perc}} = \left[q_{\alpha/2}(\hat\theta^*),\, q_{1-\alpha/2}(\hat\theta^*)\right],$$

ahol $\alpha=0{,}05$ esetén $q_{0{,}025}$ és $q_{0{,}975}$ a 2,5%-os és 97,5%-os kvantilis. Ez az intervallum ott működik jól, ahol $\hat\theta^*$ eloszlása megközelítőleg szimmetrikus a $\hat\theta$ pontbecslés körül, és a becslő varianciája nagyjából állandó a lehetséges $\theta$-értékek mentén.

A BCa (bias-corrected and accelerated, torzítás-korrigált és gyorsított) intervallum két korrekciót vezet be. Az első a torzítás-korrekció:

$$\hat z_0 = \Phi^{-1}\!\left(\frac{1}{B}\sum_{b=1}^B I[\hat\theta^*_b < \hat\theta]\right),$$

ahol $\Phi^{-1}$ a standard normális eloszlás kvantilisfüggvénye, $I[\cdot]$ pedig indikátor. Ha a bootstrap-replikációk fele kisebb $\hat\theta$-nál, akkor $\hat z_0=0$. Ha a replikációk mediánja eltolódott, $\hat z_0$ nem nulla.

A második, $\hat a$ gyorsító (acceleration) tag a becslő érzékenységét méri az egyes megfigyelésekre. A jackknife becslés $\hat\theta_{(i)}$ az $i$-edik megfigyelés elhagyásával számolt statisztika, $\hat\theta_{(\cdot)}$ ezek átlaga:

$$\hat a = \frac{\sum_{i=1}^{n}\left(\hat\theta_{(\cdot)}-\hat\theta_{(i)}\right)^3} {6\left[\sum_{i=1}^{n}\left(\hat\theta_{(\cdot)}-\hat\theta_{(i)}\right)^2\right]^{3/2}}.$$

A BCa a replikációk eloszlásán a $\Phi$-transzformált kvantiliseket igazítja: az $1-\alpha$ szintű intervallum végpontjait a replikációk $\alpha_1$ és $\alpha_2$ kvantilisei adják, ahol

$$\alpha_1 = \Phi\!\left(\hat z_0+\frac{\hat z_0+z_{\alpha/2}}{1-\hat a(\hat z_0+z_{\alpha/2})}\right), \quad \alpha_2 = \Phi\!\left(\hat z_0+\frac{\hat z_0+z_{1-\alpha/2}}{1-\hat a(\hat z_0+z_{1-\alpha/2})}\right).$$

Ha $\hat z_0=\hat a=0$, a BCa visszaesik a percentilis intervallumra. A $\hat a$ számítása ugyanazt a gondolatot használja, mint a derivált: kis, egyetlen megfigyelésre szorítkozó perturbációra adott érzékenységet vizsgál — visszautalva a Láncszabály és számítási gráf alfejezet lokális változási rátájára. A BCa jobban kezeli a ferde eloszlásokat, de érzékenyebb a jackknife-becslések simaságára is.

15.3.3 Fogalmi magyarázat: páros és klaszteres újramintázás

Ha a megfigyelési egység pár, a páros újramintázás (paired resampling) a párokat, nem az egyedi értékeket keveri. Ha minden $i$ ügyfélhez tartozik egy előtte/utána pár: $d_i = y_{i,1} - y_{i,0}$, akkor a bootstrap az $i$ párindexekből húz visszatevéssel, és a $\hat\theta^*$ statisztikát a kiválasztott $d_i$ értékeken számolja. A $d_i$ vektor alakja $(n_p,)$, ahol $n_p$ az értékelhető párok száma. A párok felbontása — tehát az előtte-érték és az utána-érték külön-külön történő visszatevéses húzása — tönkreteszi a párhatást, mert olyan kombinációkat is megenged, amelyek sosem fordultak elő egy ügyfélnél.

A klaszer-bootstrap (cluster bootstrap) értelemszerűen akkor szükséges, ha a mintában klaszterek, például ügyfelek vannak. Itt a rendelési sorok nem cserélhetők szabadon, mert a customer_id közös hatást visz be. A húzott egység ekkor a klaszter, nem a sor. Ha $m$ ügyfél van, és a $j$-edik ügyfélnek $n_j$ rendelése, akkor minden bootstrap-replikációban $m$ klasztert húzunk visszatevéssel, majd az összes kiválasztott ügyfél összes sorát bevesszük a mintába:

$$\hat\theta^*_b = \frac{\sum_{j \in \mathcal J^*_b}\sum_{i=1}^{n_j} y_{ji}}{\sum_{j \in \mathcal J^*_b} n_j}.$$

Itt $\mathcal J^*_b$ a $b$-edik húzásban szereplő klaszterindexek multihalmaza, elemszáma $m$. Figyeljünk arra, hogy a húzott sorok száma reprodukciónként változhat, mert az ügyfelek rendelésszáma különböző. A pontbecslés továbbra is a rendelésszintű átlag; a standard error viszont már tükrözi a klasztereken belüli függést. Ha a delivery_days ügyfélen belül hasonló, a sor-szintű bootstrap alábecsüli a bizonytalanságot, a klaszter-bootstrap pedig szélesebb intervallumot ad.

15.3.4 RNG és degeneráltság

A reprodukálhatósághoz a véletlenszám-generátort (rng, random number generator) explicitül meg kell adni. A np.random.default_rng(seed) modern Generator példányt hoz létre; a scipy.stats.bootstrap random_state argumentuma elfogad ilyen példányt is. A régi np.random.RandomState használata kerülendő, mert elavult a NumPy új projektjeiben. A seed rögzítése mellett ugyanazon a gépen és csomagverzión belül a húzások reprodukálhatók; más NumPy-verzió vagy platform esetén bitazonos eredményt nem ígérhetünk. Statisztikailag a lényeg az, hogy a Monte Carlo ingadozás $O(1/\sqrt B)$ nagyságrendű legyen, ezért a $B=1000$–$2000$ replikáció ésszerű intervallumstabilitást ad.

Degeneráltság (degeneracy) akkor lép fel, ha a bootstrap-replikációk eloszlása túl kevés egyedi értékre koncentrálódik. Klasszikus példa a két elemből álló minta mediánja: az $x=(1,100)$ elemből húzott $n=2$ bootstrap-minták mediánjai közt csak az $1$, $50{,}5$ és $100$ értékek fordulnak elő. A bootstrap eloszlás ilyenkor lépcsős, a percentilis intervallum nem folytonos, és a lefedettség értelmezése sérül. Ugyanezt látjuk extrém kvantiliseknél, például a minimum bootstrap becslésénél. A degeneráltság figyelmeztetés: a bootstrap nem varázsol folytonos valószínűségi eloszlást egy túl kicsi mintából.

15.3.5 Teljes minimális munkapélda

A következő szkript a lab/ munkakönyvtárból futtatható: python bootstrap_delivery.py. Az input a tiszta orders.csv; a cél a delivery_days átlagának és mediánjának intervalluma. A seed 20261020, a sor-szintű replikációk száma $B=2000$, a klaszter-bootstrapé $B=1000$. A SciPy bootstrap hívás magas szintű API, de a klaszter-változathoz kézi hurkot írunk, mert a SciPy nem tud customer_id szerinti klasztert húzni.

import numpy as np
import pandas as pd
from scipy import stats

orders = pd.read_csv("orders.csv", parse_dates=["order_ts"])
y = orders["delivery_days"].to_numpy(dtype=float)
customer = orders["customer_id"].to_numpy()
n = y.size

print(f"rows={n}")
print(f"mean={y.mean():.3f}, median={np.median(y):.1f}, "
      f"sd={y.std(ddof=1):.3f}")

# Sor-szintű bootstrap: itt minden rendelési sort függetlennek tételezünk fel.
rng = np.random.default_rng(20261020)
B = 2000

boot_mean_perc = stats.bootstrap(
    (y,), np.mean, n_resamples=B, method="percentile",
    random_state=rng, confidence_level=0.95,
)
boot_mean_bca = stats.bootstrap(
    (y,), np.mean, n_resamples=B, method="BCa",
    random_state=rng, confidence_level=0.95,
)
boot_med_perc = stats.bootstrap(
    (y,), np.median, n_resamples=B, method="percentile",
    random_state=rng, confidence_level=0.95,
)

for name, res in [
    ("row mean percentile", boot_mean_perc),
    ("row mean BCa", boot_mean_bca),
    ("row median percentile", boot_med_perc),
]:
    lo, hi = res.confidence_interval
    print(f"{name:22s} [{lo:.3f}, {hi:.3f}] "
          f"se={res.standard_error:.3f}")

# Klaszter-bootstrap: a húzott egység az ügyfél, az összes rendelésével együtt.
customers, inverse = np.unique(customer, return_inverse=True)
m = customers.size
cluster_rows = [np.flatnonzero(inverse == j) for j in range(m)]

B_cluster = 1000
cluster_mean = np.empty(B_cluster)
cluster_median = np.empty(B_cluster)
for b in range(B_cluster):
    chosen = rng.integers(0, m, size=m)
    rows = np.concatenate([cluster_rows[j] for j in chosen])
    sample = y[rows]
    cluster_mean[b] = sample.mean()
    cluster_median[b] = np.median(sample)

for name, arr in [
    ("cluster mean percentile", cluster_mean),
    ("cluster median percentile", cluster_median),
]:
    lo, hi = np.quantile(arr, [0.025, 0.975])
    print(f"{name:22s} [{lo:.3f}, {hi:.3f}] "
          f"se={arr.std(ddof=1):.3f}")

A pontos számok a helyi futtatástól függenek, ezért az alábbi megjegyzés szemléltető, nem mért kimenet: ha a sor-szintű átlag-intervallum kerekítve 4,6–4,8 körüli, a klaszter-változat ezt gyakran 0,1–0,2 nappal szélesebbre nyitja. A medián intervalluma a nagy minta miatt már elfogadhatóan sima, de enyhe fokon továbbra is lépcsős lehet.

15.3.6 Eredmény értelmezése

A számok mögötti döntés a húzott egység. A sor-szintű bootstrap minden replikációban $n$ rendelést húz, és úgy kezeli a sorokat, mintha függetlenek lennének. A klaszter-bootstrap ezzel szemben minden replikációban $m$ ügyfelet húz, és az ügyfél teljes rendeléstörténetét viszi be. A két eljárás ugyanazt az eredeti rendelési átlagot használja kiinduló pontbecslésként. A klaszter-bootstrap replikációinak rendelésenként súlyozott átlaga azonban véletlen klaszterméretek mellett nem szükségképpen torzítatlan a megfigyelt átlagra; az értelmezés a célpopulációtól és a klaszterméretek kezelésétől is függ. A standard error viszont akkor különbözik, ha az ügyfelen belüli szállítási idők korrelálnak.

Ellenőrző invariánsok:

• A klaszter-bootstrap minden replikációban pontosan $m$ klasztert húz, de a kiválasztott sorok száma nem feltétlenül $n$.
• A replikációk átlaga közelíti a megfigyelt pontbecslést; nagy eltérés hiba jele.
• A cluster_rows lista minden eleme nem üres, mert minden customer_id jelen van az orders táblában.
• A BCa akkor tér el érdemben a percentilistől, ha $\hat z_0$ vagy $\hat a$ nullától távol van. A delivery_days enyhén jobbra ferde lehet, ezért a két átlag-intervallum nem feltétlenül egyezik meg.

Klaszteres újramintázásnál külön rögzítsük, hogy ügyfélátlagot vagy rendelésátlagot kívánunk becsülni. Egyenlőtlen klaszterméretek mellett ez más súlyozást és eltérő célparamétert jelenthet.

15.3.7 Feladatok

1. Értelmezés. Érveljen amellett, hogy a delivery_days esetén miért az ügyfél a helyes bootstrap-egység! Mondjon egy olyan feltételt, amely mellett a sor-szintű és a klaszter-bootstrap azonos standard errorhoz vezet.

2. Implementálás. A Northstar orders.csv-ből válassza ki azokat az ügyfeleket, akiknek legalább két rendelésük van. Párosítsa minden ilyen ügyfél első két rendelését időrendben, számítsa ki a $d_i$ = második delivery_days − első delivery_days különbséget, majd adjon páros resamplinggel 95%-os percentilis intervallumot az átlagos $d$-re. Használjon $B=2000$ replikációt és rögzített seedet.

3. Hibakeresés. Valaki a páros feladatot így próbálta megoldani:

first_boot = first[rng.integers(0, first.size, size=first.size)]
second_boot = second[rng.integers(0, second.size, size=second.size)]
repl_b = second_boot.mean() - first_boot.mean()

Mi a hiba? Hogyan változtatja ez meg az átlagos különbség bizonytalanságát? Javítsa a kódot.

4. Megváltozott feltétel. A külön kísérleti assignments.csv és outcomes.csv 300 ügyfélen a conversion_7d kimenetet tartalmazza. A beavatkozási cél a treatment és control értékek átlagának különbsége. Ebben az adatban minden ügyfél pontosan egyszer randomizált, tehát a húzott egység az ügyfél. Készítsen csoporton belüli bootstrapet $B=2000$ replikációval a hatáskülönbség 95%-os percentilis intervallumára. Ellenőrizze, hogy a nevező minden replikációban tartalmazhatja ugyanazt az ügyfelet többször is, de nem használhatja a rendelési sorokat.

15.3.8 Megoldások és indoklás

1. megoldás. Az ügyfél a helyes egység, mert a rendelési sorok ugyanazon ügyfélhez tartozó kimenetei közös ügyfélhatást hordozhatnak; ezt a függetlenség és csoport alfejezetben az ICC fogalmával láttuk. A sor-szintű bootstrap akkor adna hasonló standard errort, ha a delivery_days ügyfélen belül független volna, és az ügyfélhatás nem befolyásolná a kimenetet. Ha a klaszteren belüli korreláció nulla és a húzott sorok száma nagyjából azonos marad, a két replikációeloszlás közelíti egymást. Az ellenőrző gondolat: tegyük fel, hogy minden ügyfélhez konstans szállítási idő tartozna; ilyenkor a sor-szintű szórás kicsi, a klaszter-bootstrap viszont az ügyfelek közötti szórást tükrözi.

2. megoldás. A groupby művelet először azokat az ügyfeleket tartja meg, akiknek legalább két rendelése van, majd a rendezett order_ts alapján az első két rendelést párosítja.

import numpy as np
import pandas as pd

orders = pd.read_csv("orders.csv", parse_dates=["order_ts"])
customer_id = "customer_id"
eligible = orders.groupby(customer_id).size()
eligible = eligible[eligible >= 2].index
sub = orders[orders[customer_id].isin(eligible)].sort_values(
    [customer_id, "order_ts"]
)
first = sub.groupby(customer_id)["delivery_days"].first().to_numpy()
second = sub.groupby(customer_id)["delivery_days"].nth(1).to_numpy()
d = second - first

rng = np.random.default_rng(20261020)
B = 2000
repl = np.empty(B)
for b in range(B):
    idx = rng.integers(0, d.size, size=d.size)
    repl[b] = d[idx].mean()

lo, hi = np.quantile(repl, [0.025, 0.975])
print(f"pairs={d.size}, mean diff={d.mean():.3f} "
      f"CI=[{lo:.3f}, {hi:.3f}]")

Az ellenőrző invariant: d.size pontosan azoknak az ügyfeleknek a száma, akiknek legalább két rendelésük van; a replikációk átlaga közel van d.mean()-hez. A szelektált részsokaság nem reprezentálja a csak egy rendeléssel rendelkező ügyfeleket — ezt explicit feltételként kell kimondani a jelentésben.

3. megoldás. A hibás kód a first és a second tömböket külön-külön húzza újra. Ezzel a párok összerendelése elvész: egy ügyfél első rendeléséhez egy másik ügyfél második rendelése is kerülhet. Ha a páron belüli értékek pozitívan korrelálnak, a párosítás csökkenti a különbség szórását, tehát a hibás verzió szélesebb, konzervatívabb intervallumot ad, de rossz célt mér: olyan kombinációk is szerepelnek, amelyek a valós ügyfélúton nem fordultak elő. A javítás az, hogy a $d_i$ különbségvektort húzzuk újra, mint a 2. megoldásban. Párhuzamos ellenőrzés: a javított kód replikációinak szórása kisebb, mint a hibás kódból számolt szórás, ha a párkorreláció pozitív.

4. megoldás. A külön kísérletben minden ügyfél egyszer szerepel, így a megfelelő húzás a customer_id szintű, csoporton belüli bootstrap:

import numpy as np
import pandas as pd

assign = pd.read_csv("assignments.csv")
out = pd.read_csv("outcomes.csv")
df = assign.merge(out, on="customer_id", validate="one_to_one")
control = df.loc[df["variant"] == "control", "conversion_7d"].to_numpy()
treatment = df.loc[df["variant"] == "treatment", "conversion_7d"].to_numpy()

rng = np.random.default_rng(20261020)
B = 2000
repl = np.empty(B)
for b in range(B):
    c = control[rng.integers(0, control.size, size=control.size)]
    t = treatment[rng.integers(0, treatment.size, size=treatment.size)]
    repl[b] = t.mean() - c.mean()

lo, hi = np.quantile(repl, [0.025, 0.975])
print(f"ITT estimate={treatment.mean()-control.mean():.4f} "
      f"CI=[{lo:.4f}, {hi:.4f}]")

Az ITT-nevező az összes randomizált ügyfél; a replikációk csak ügyfeleket húznak. A replikációk átlaga közelíti a becsült hatás nagyságát; ha ettől a repl átlaga nagyon eltér, a húzási logikában hiba van. A merge validate="one_to_one" ellenőrzi, hogy a táblák valóban ügyfélszintűek.

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

A leggyakoribb félreértés az, hogy a bootstrap korrigálja a rossz adatfelvételi folyamatot. Ez nem igaz. Ha egy elégedettségi mérést csak panaszos ügyfelek töltenek ki, akkor a 95%-os bootstrap-intervallum továbbra is a panaszos rétegre vonatkozik, nem az összes ügyfélre. A $B$ növelése csak a Monte Carlo hibát csökkenti, a szelekciós torzítást nem.

A bootstrap a nagyon kicsi mintát sem varázsolja megbízhatóvá. Az $x=(1,100)$ kételemes minta mediánjának bootstrap-eloszlása mindössze három lehetséges értékre támaszkodik; bármilyen nagy $B$ mellett is degenerált marad. A delivery_days laborban is ellenőrizni kell, hogy a replikációk eloszlása elég sok egyedi értéket tartalmaz-e; a kis minta vagy extrém kvantilis degeneráltsága megbízhatatlanná teszi a percentilis intervallumot.

15.3.10 Önellenőrzés

Kérdezze meg magától: mi a húzott egység a saját feladatomban, sor vagy ügyfél? A scipy.stats.bootstrap magas szintű hívása a sor- vagy elemszintű függetlenséget feltételezi; klaszter esetén kézi hurkot kell írni. A BCa-t használva ellenőrizze, hogy a jackknife-becslések simán változnak-e; a medián ilyen szempontból nem ideális, mert a jackknife érzékenysége ugrálhat. A seed rögzítése a statisztikai reprodukálhatóság minimális feltétele; a Bináris csomag, BLAS és conda-döntés alfejezet tanulsága szerint a natív SciPy-környezet verzióját is rögzíteni kell, ha valaki bitazonosan szeretné újrafuttatni a kódot. A konfidenciaszint továbbra is az eljárás ismételt mintavételi lefedettségére vonatkozik, nem a konkrét intervallum valószínűségére, ahogyan ezt a Standard error és konfidenciaintervallum alfejezetben tárgyaltuk.

15.3.11 Rövid összefoglaló

A bootstrap a megfigyelt mintából visszatevéses húzással közelíti a becslő mintavételi eloszlását. A percentilis intervallum a replikációk kvantiliseit veszi, a BCa emellett torzítás- és ferdeségkorrekciót alkalmaz. Páros adatnál a párokat, ismételt ügyféladatnál az ügyfelet kell újramintázni; a sor-szintű húzás a klaszteren belüli függés miatt alábecsüli a standard errort. A seed rögzítése szükséges, a degenerált replikációeloszlás pedig figyelmeztetés a kis mintára vagy extrém statisztikára. A következő alfejezetben a predikciós, paraméter- és Monte Carlo intervallumok közötti fogalmi határokat tisztázzuk.