7.2 Bootstrap for standardfeil og simulert dekning

Estimér standardfeilen du ikke kan regne — parametrisk og ikke-parametrisk bootstrap — og simuler dekningsgraden til et konfidensintervall.

60 min

8 oppgaver

Bootstrap for standardfeilsimulert dekning

Din fremgang i kapitlet

0 / 8 oppgaver

Kapitlets plass i kurset

Bygger på

7.1Inversjonsmetoden og numpy-verktøykassen 5.1Estimatorer og egenskaper: forventningsretthet, varians og skjevhet

Brukes videre i

7.PPrøver til del 7: Simulering, bootstrap og Python

Bootstrap er ikke et nytt fenomen på UiO slik det er ved NTNU — det har vært fast inventar i et helt tiår. Kan du algoritmen og koden, plukker du opp poeng de fleste går glipp av.

- Bootstrap for standardfeil (sjanger I): ikke-parametrisk allerede i V2015 og V2016; parametrisk i V2019 og V2021; full kjørbar Python-kode fra V2020.
- Simulert dekningsgrad (eksakt pivotal-KI mot Wald-KI) i V2024 — den bygger direkte på dekningstallene fra kap. 6.2.
- «Sjanger I» = simulering, bootstrap og Python (samme kode som i kap. 7.1).
- Sensorkrav 10: bootstrap-algoritmen skal beskrives presist og koden være kjørbar med ddof=1. All kode her er faktisk kjørt.

Prioritet: kunne (nivå 2). Sammen med kap. 7.1 er dette det som skiller en A fra en B når den analytiske standardfeilen er vanskelig eller umulig å regne ut for hånd.

Forkunnskaper: Bygger på kap. 7.1 (simulering, np.random.uniform, vektorisering, ddof=1) og kap. 5.1 (standardfeil). Dekningsdelen henter inn det eksakte pivotal-KI-et fra kap. 6.2 og den analytiske standardfeilen fra kjeden i kap. 6.3.

Sist du var her — det du må ha friskt:

1. Standardfeilen til en estimator er standardavviket til estimatoren: $\text{SE}(\hat\theta) = \sqrt{\text{Var}(\hat\theta)}$ (fra kap. 5.1). Den forteller hvor mye estimatet ville variert fra utvalg til utvalg.
2. ddof=1 i np.std/np.var gir den forventningsrette empiriske variansen (divisor $n-1$ ) — vi bruker den til å regne ut bootstrap-SE-en.
3. Eksakt pivotal-KI for eksponensial (forventning $\mu$ , fra kap. 6.2): siden $2n\bar X/\mu \sim \chi^2_{2n}$ , er et $95\%$ -KI $\left[\dfrac{2n\bar X}{\chi^2_{2n,\,0{,}025}},\ \dfrac{2n\bar X}{\chi^2_{2n,\,0{,}975}}\right]$ , der $\chi^2_{2n,\,0{,}025}$ og $\chi^2_{2n,\,0{,}975}$ er øvre- og nedre-halekvantilene. Dette bruker vi i dekningsstudien til slutt.

Iblant har du en estimator $\hat\theta$ der du enkelt kan regne ut estimatet, men der den analytiske standardfeilen er stygg eller umulig å utlede — for eksempel medianen, eller en komplisert funksjon av flere parametre. Bootstrap løser dette ved en enkel idé: lat som utvalget ditt er hele befolkningen, trekk mange nye «utvalg», og se hvor mye estimatoren spriker. Spredningen i disse gjentatte estimatene er standardfeilen.

Vi møter to varianter (ikke-parametrisk og parametrisk), pluss det nært beslektede Monte Carlo-prinsippet: en sannsynlighet er en forventet indikator, så den estimeres med et gjennomsnitt. Til slutt setter vi alt sammen i en dekningsstudie som viser hvorfor et eksakt KI treffer 95 % mens et Wald-KI bommer for liten $n$ — mekanismen bak tallene i kap. 6.2.

Fire læringsløkker à ~13 min — ta en pause etter løkke 2.

Løkke 1 — Ikke-parametrisk bootstrap (~14 min)

Bootstrap-algoritmen (presist): For

b = 1, 2, \ldots, B

1. Trekk et nytt «utvalg» $x^*_1, \ldots, x^*_n$ av samme størrelse $n$ som dataene.
2. Regn ut estimatoren på nytt: $\hat\theta^*_b = \text{estimator}(x^*_1, \ldots, x^*_n)$ .

Til slutt: estimer standardfeilen som det empiriske standardavviket til de $B$ verdiene:

$\widehat{\text{SE}}(\hat\theta) = \sqrt{\frac{1}{B-1}\sum_{b=1}^{B}\big(\hat\theta^*_b - \bar{\hat\theta^*}\big)^2} = \texttt{np.std(theta\_b, ddof=1)}.$

Ikke-parametrisk bootstrap trekker de nye utvalgene med tilbakelegging fra selve dataene (np.random.choice(data, size=n, replace=True)). Vi antar altså ingen bestemt fordeling — dataene får «tale for seg selv». Nøkkelen er replace=True: uten tilbakelegging ville hvert trekk bare vært en omstokking av de samme tallene, og all variasjon forsvant.

Merk skillet $B$ mot $n$ : $n$ er utvalgsstørrelsen (fast, gitt av dataene); $B$ er hvor mange bootstrap-trekk vi gjør (vi velger den selv, typisk 1 000–10 000). Flere $B$ gir bare et mer presist SE-estimat — det endrer ikke $n$ .

✏️Eksempel 1: Ikke-parametrisk bootstrap-SE for medianen

Ti målte konsentrasjoner (mg/l) er $2{,}1,\ 5{,}4,\ 1{,}2,\ 8{,}9,\ 3{,}3,\ 6{,}7,\ 0{,}8,\ 4{,}5,\ 12{,}1,\ 2{,}9$ . Du vil estimere standardfeilen til medianen — men det finnes ingen enkel formel. Bruk ikke-parametrisk bootstrap med $B = 10000$ .

Medianen har ingen pen analytisk standardfeil, så bootstrap er det naturlige verktøyet. Vi trekker med tilbakelegging fra de ti tallene, regner medianen på nytt hver gang, og tar standardavviket:

import numpy as np

data = np.array([2.1, 5.4, 1.2, 8.9, 3.3, 6.7, 0.8, 4.5, 12.1, 2.9])
n = len(data)

np.random.seed(0)
B = 10000
theta_b = np.empty(B)
for b in range(B):
    xb = np.random.choice(data, size=n, replace=True)   # med tilbakelegging!
    theta_b[b] = np.median(xb)

print(np.median(data), np.std(theta_b, ddof=1))
# 3.9  1.3634...

Medianen i de opprinnelige dataene er $3{,}9$ , og bootstrap-standardfeilen er $\approx 1{,}36$ . Tolkning: hadde vi samlet inn nye datasett av samme størrelse, ville medianestimatet typisk variert med rundt $1{,}4$ mg/l.

Kritisk detalj: replace=True og ddof=1. Fjerner du replace=True, blir hvert xb bare en omstokking av de samme ti tallene, medianen er (nesten) konstant, og SE-en kollapser feilaktig mot 0.

📝Oppgave 1

(Innsteg.) Beskriv i ord de tre stegene i ikke-parametrisk bootstrap for å estimere standardfeilen til en estimator $\hat\theta$ , og si eksplisitt hvordan de nye utvalgene trekkes.

📝Oppgave 2

En student skriver denne bootstrap-koden for SE-en til gjennomsnittet, men får en mistenkelig liten verdi:

theta_b = np.empty(B)
for b in range(B):
    xb = np.random.choice(data, size=n, replace=False)
    theta_b[b] = np.mean(xb)
se = np.std(theta_b, ddof=1)

a) Hva er feilen, og hvorfor blir SE-en (nesten) 0?

b) Rett koden.

Løkke 2 — Parametrisk bootstrap (~14 min)

— naturlig pausepunkt etter denne løkka —

Parametrisk bootstrap bruker samme algoritme, men trekker de nye utvalgene fra den tilpassede fordelingen — med parameterestimatet

\hat\theta

satt inn — i stedet for fra dataene direkte. Har du grunn til å tro at dataene er (for eksempel) eksponensialfordelte, tilpasser du fordelingen (

\hat\mu = \bar x

), og trekker nye utvalg fra eksponensial

(\hat\mu)

med inversjonsmetoden fra kap. 7.1.

Trekker nye utvalg fra …
Ikke-parametrisk	dataene, med tilbakelegging
Parametrisk	den tilpassede fordelingen (med $\hat\theta$ satt inn)

Når velge hva? Parametrisk bootstrap er mer presis hvis fordelingsantakelsen er riktig (den utnytter modellen); ikke-parametrisk er tryggere når du er usikker på fordelingen. Begge gir SE-en som np.std(theta_b, ddof=1).

✏️Eksempel 2: Parametrisk bootstrap for eksponensialmodellen

De ti målingene fra eksempel 1 antas nå eksponensialfordelte med forventning $\mu$ (forventnings-parametrisering, $\displaystyle f(x) = \frac{1}{\mu}e^{-x/\mu}$ ). Estimatoren er $\hat\mu = \bar X$ . Bruk parametrisk bootstrap ( $B = 10000$ ) til å estimere standardfeilen til $\hat\mu$ , og sammenlign med den analytiske $\text{SE}(\bar X) = \mu/\sqrt{n}$ .

Vi tilpasser modellen ( $\hat\mu = \bar x$ ), trekker nye eksponensiale utvalg fra eksponensial $(\hat\mu)$ med inversjonsformelen $-\hat\mu\ln(1-u)$ , og regner $\bar X$ på nytt:

import numpy as np

data = np.array([2.1, 5.4, 1.2, 8.9, 3.3, 6.7, 0.8, 4.5, 12.1, 2.9])
n = len(data)
muhat = np.mean(data)

np.random.seed(0)
B = 10000
theta_b = np.empty(B)
for b in range(B):
    xb = -muhat * np.log(1 - np.random.uniform(size=n))   # fra tilpasset eksp(muhat)
    theta_b[b] = np.mean(xb)

print(muhat, np.std(theta_b, ddof=1), muhat / np.sqrt(n))
# 4.79  1.5072...  1.5147...

Estimatet er $\hat\mu = 4{,}79$ . Parametrisk bootstrap-SE er $\approx 1{,}51$ , og den analytiske $\hat\mu/\sqrt{n} = 4{,}79/\sqrt{10} \approx 1{,}515$ . De to stemmer nesten perfekt — som de skal, siden vi bootstrapper fra nøyaktig den modellen den analytiske formelen forutsetter. Her kan vi regne SE-en analytisk; poenget er at bootstrap gir samme svar og fungerer også når den analytiske formelen mangler.

📝Oppgave 3

Forklar forskjellen på ikke-parametrisk og parametrisk bootstrap i én setning hver, og gi ett argument for når du ville foretrukket den ikke-parametriske varianten.

📝Oppgave 4

Levetider antas eksponensialfordelte med forventning $\mu$ , og du vil estimere standardfeilen til estimatoren $\hat\theta = \hat\mu \ln 2$ (et estimat av medianen, siden medianen er $\mu\ln 2$ ). Du har $n = 8$ observasjoner med gjennomsnitt $\bar x = 6{,}0$ .

a) Skriv parametrisk bootstrap-kode ( $B = 10000$ ) som estimerer $\text{SE}(\hat\theta)$ .

b) Den analytiske SE-en er $\ln 2 \cdot \mu/\sqrt{n}$ . Hva bør bootstrap-anslaget ligge nær?

Løkke 3 — Monte Carlo som inferensverktøy (~14 min)

Bak all simulering ligger ett prinsipp: en sannsynlighet er en forventet indikator. Hvis

I

er 1 når en hendelse inntreffer og 0 ellers, er

P(\text{hendelse}) = E(I)

. Og en forventning estimeres med et gjennomsnitt. Derfor:

$\widehat{P(\text{hendelse})} = \frac{1}{B}\sum_{b=1}^{B} I_b = \texttt{np.mean(betingelse)}.$

Dette Monte Carlo-estimatet har fine egenskaper:

- Det er forventningsrett: $E\big[\widehat P\big] = p$ (gjennomsnittet av forventningsrette indikatorer).
- Variansen er $\dfrac{p(1-p)}{B}$ (indikatorene er uavhengige med varians $p(1-p)$ ).
- Derfor er standardfeilen $\sqrt{p(1-p)/B}$ , som avtar som $1/\sqrt B$ . Vil du halvere usikkerheten, må du firedoble $B$ .

Merk igjen: presisjonen styres av $B$ (antall trekk), ikke av $n$ .

✏️Eksempel 3: Monte Carlo-estimat av en halesannsynlighet

En komponent har eksponensialfordelt levetid med forventning $\mu = 5$ år. Estimer $P(X > 8)$ med Monte Carlo ( $B = 10^6$ ), oppgi Monte Carlo-standardfeilen, og sammenlign med den eksakte verdien.

Vi trekker $10^6$ eksponensiale levetider og tar gjennomsnittet av indikatoren x > 8:

import numpy as np

np.random.seed(0)
B = 10**6
x = -5 * np.log(1 - np.random.uniform(size=B))
p = np.mean(x > 8)
print(p, np.sqrt(p * (1 - p) / B))
# 0.20242...  0.000402...

Monte Carlo-estimatet er $\widehat P \approx 0{,}2024$ med standardfeil $\approx 0{,}0004$ . Eksakt: $P(X > 8) = e^{-8/5} = e^{-1{,}6} \approx 0{,}2019$ . Estimatet ligger godt innenfor et par standardfeil av sannheten.

La merke til $1/\sqrt B$ -loven: med $B = 10^6$ er SE-en på fjerde desimal. Ville vi ned til $0{,}0002$ , måtte $B$ opp til $4\cdot 10^6$ .

📝Oppgave 5

Et Monte Carlo-estimat av en sannsynlighet $p$ bygger på $B = 10000$ uavhengige trekk. Estimatet blir $\widehat p = 0{,}30$ .

a) Hva er Monte Carlo-standardfeilen?

b) Omtrent hvor stor må $B$ være for å halvere denne standardfeilen?

📝Oppgave 6

(Kodelesing.) Forklar hva denne koden estimerer, og hvorfor np.mean av en betingelse gir en sannsynlighet:

import numpy as np
from scipy import stats
np.random.seed(0)
y = stats.norm.rvs(10, 2, size=10**6)
print(np.mean((y > 8) & (y < 12)))

Løkke 4 — Simulert dekningsgrad: eksakt mot Wald (~14 min)

Et $95\%$ -konfidensintervall skal, over mange gjentatte utvalg, dekke den sanne parameteren $95\%$ av gangene. Simulert dekningsgrad sjekker om det stemmer:

1. Generer mange datasett fra en fordeling med kjent sann parameter.
2. Bygg KI-et i hvert datasett.
3. Tell andelen KI-er som faktisk dekker den sanne parameteren.

Dette er mekanismen bak dekningstallene fra kap. 6.2: det eksakte pivotal-KI-et treffer $95\%$ (fordi $2n\bar X/\mu \sim \chi^2_{2n}$ er en eksakt pivotal), mens Wald-KI-et $\bar X \pm 1{,}96\,\text{SE}$ underdekker for liten $n$ — fordi $\bar X$ for eksponensiale data er høyreskjevt, ikke symmetrisk normalfordelt slik Wald antar.

«Dekker»-indikatoren er bare lo <= mu <= hi, og dekningsgraden er np.mean av den — igjen: en sannsynlighet er en forventet indikator.

✏️Eksempel 4: Dekningsstudie for eksponensialmodellen, n = 10

For eksponensialfordelte data med forventning $\mu = 5$ og $n = 10$ : simuler dekningsgraden til (i) det eksakte pivotal-KI-et og (ii) Wald-KI-et $\bar X \pm 1{,}96\,\bar X/\sqrt{n}$ , over $M = 20000$ datasett. Forklar forskjellen.

Vi lager $20000$ datasett fra eksponensial $(\mu=5)$ , bygger begge KI-ene i hvert, og teller andelen som dekker $\mu = 5$ :

import numpy as np
from scipy import stats

np.random.seed(0)
n, mu, M = 10, 5.0, 20000
dekk_eksakt = 0
dekk_wald = 0
for _ in range(M):
    x = -mu * np.log(1 - np.random.uniform(size=n))
    xbar = np.mean(x)
    # eksakt pivotal:  2*n*xbar/mu ~ chi^2_{2n}
    lo = 2*n*xbar / stats.chi2.ppf(0.975, 2*n)
    hi = 2*n*xbar / stats.chi2.ppf(0.025, 2*n)
    if lo <= mu <= hi:
        dekk_eksakt += 1
    # Wald:  xbar +/- 1.96 * SE, med SE = xbar/sqrt(n)
    se = xbar / np.sqrt(n)
    if xbar - 1.96*se <= mu <= xbar + 1.96*se:
        dekk_wald += 1

print(dekk_eksakt/M, dekk_wald/M)
# 0.9508  0.9015

Resultat: det eksakte KI-et dekker $\approx 95{,}1\%$ — praktisk talt nominelt $95\%$ , som det skal siden pivotalen $2n\bar X/\mu \sim \chi^2_{2n}$ er eksakt uansett $\mu$ og $n$ . Wald-KI-et dekker bare $\approx 90{,}2\%$ — det underdekker klart.

Hvorfor? Wald antar at $\bar X$ er tilnærmet normalfordelt og symmetrisk rundt $\mu$ . For $n = 10$ eksponensiale observasjoner er $\bar X$ fortsatt merkbart høyreskjevt, så det symmetriske $\pm 1{,}96\,\text{SE}$ -intervallet bommer oftere enn de lovede $5\%$ . Dette er nøyaktig asymmetrien kap. 6.2 tallfestet (der ca. $95{,}9\%$ mot $91{,}1\%$ for den oppgavens oppsett) — samme mekanisme, små forskjeller skyldes ulikt antall gjentak og tilfeldig frø. Konklusjonen er den samme: for liten $n$ er det eksakte pivotal-KI-et å foretrekke.

📝Oppgave 7

I dekningskoden over står linjen

lo = 2*n*xbar / stats.chi2.ppf(0.975, 2*n)
hi = 2*n*xbar / stats.chi2.ppf(0.025, 2*n)

a) Forklar hvorfor den nedre grensen lo bruker den øvre kvantilen chi2.ppf(0.975, 2n).

b) Hvor mange frihetsgrader har kjikvadratfordelingen her, og hvorfor?

📝Oppgave 8

Du skal undersøke hvordan Wald-KI-ets underdekning avhenger av utvalgsstørrelsen for eksponensiale data.

a) Beskriv hvordan du ville endre dekningskoden for å studere $n = 10, 30, 100$ i stedet for bare $n = 10$ .

b) Hva forventer du skjer med Wald-dekningen når $n$ vokser, og hvorfor? Hva skjer med den eksakte dekningen?

Typiske feil

- Glemme ddof=1 i bootstrap-SE-en. np.std(theta_b, ddof=1) deler på $B-1$ og gir standardfeilen; uten ddof=1 deles det på $B$ .
- Beskrive bootstrap upresist. Si eksplisitt at ikke-parametrisk bootstrap trekker med tilbakelegging fra dataene, og at parametrisk bootstrap trekker fra den tilpassede fordelingen (med $\hat\theta$ satt inn). En beskrivelse uten dette taper poeng (sensorkrav 10).
- Forveksle $B$ og $n$ . $B$ = antall bootstrap-/Monte Carlo-trekk (styrer presisjonen, avtar som $1/\sqrt B$ ); $n$ = utvalgsstørrelsen (fast, gitt av dataene). Flere $B$ gjør ikke datasettet ditt større.
- replace=False i ikke-parametrisk bootstrap. Uten tilbakelegging blir hvert «utvalg» en ren omstokking, og SE-en kollapser mot 0.
- Simulere dekning under feil fordeling. Skal du studere dekningen til et KI, må datasettene genereres fra den fordelingen KI-et forutsetter — ellers måler du noe annet.

Begrepsbank

Flashcard-/repetisjonsstoff — hopp trygt over ved førstegangslesing; tidsanslaget gjelder kjernestoffet. Kortene samler algoritmene og idiomene til rask repetisjon.

Bootstrap (idéen)

Å estimere usikkerheten til en estimator ved å trekke mange nye «utvalg», regne estimatoren på nytt hver gang, og måle spredningen. Man later som utvalget er hele befolkningen. Brukes særlig når den analytiske standardfeilen er vanskelig eller umulig.

Ikke-parametrisk bootstrap

Bootstrap der de nye utvalgene trekkes med tilbakelegging fra selve dataene, uten noen fordelingsantakelse. Dataene «taler for seg selv». I kode: np.random.choice(data, size=n, replace=True) .

Parametrisk bootstrap

\hat\theta

`np.random.choice(..., replace=True)`

Trekker med tilbakelegging fra et array — motoren i ikke-parametrisk bootstrap. replace=True er avgjørende: uten den blir trekket bare en omstokking og all variasjon forsvinner.

Bootstrap-standardfeil

B

Bootstrap-algoritmen

b = 1, \ldots, B

B

mot

n

B

Estimatoren regnes på nytt

\hat\mu\ln 2

Monte Carlo-prinsippet

P(\text{hendelse}) = E(I)

Monte Carlo-estimat

B

Forventningsrett MC-estimat

E[\widehat P] = p

Monte Carlo-varians

p(1-p)/B

1/\sqrt B

-presisjon

1/\sqrt B

Standardfeil (repetisjon)

\text{SE}(\hat\theta) = \sqrt{\text{Var}(\hat\theta)}

Simulert dekningsgrad

Metode for å sjekke et KI: generer mange datasett med kjent sann parameter, bygg KI-et i hvert, og tell andelen som dekker parameteren. Andelen bør ligge nær den nominelle graden (f.eks. 95 %).

Dekker-indikatoren

For hvert simulerte datasett: lo <= param <= hi — 1 om KI-et dekker den sanne parameteren, 0 ellers. Dekningsgraden er np.mean av disse indikatorene over alle datasettene.

Eksakt pivotal-KI dekker nominelt

2n\bar X/\mu \sim \chi^2_{2n}

Wald-KI underdekker

\bar X \pm 1{,}96\,\text{SE}

Parametrisk vs. ikke-parametrisk — valget

Velg parametrisk når du stoler på fordelingsmodellen (mer presis, utnytter modellen); velg ikke-parametrisk når du er usikker på fordelingen (mer robust mot feilspesifikasjon).

`ddof=1` i bootstrap-SE

B-1

Repetisjonsoppgaver

Symbol- og formelliste

Dette kapitlet er skrevet av Anthropics toppmodeller (Claude Opus og Claude Fable) og er foreløpig ikke manuelt gjennomgått — kvalitetskontrollen gjøres av uavhengige KI-agenter, og innmeldte feil rettes fortløpende. Funnet en feil? Meld fra, så retter vi den. Les mer om hvordan innholdet lages.

Skolesaga er en uavhengig læringsressurs og er ikke tilknyttet eller godkjent av Universitetet i Oslo. Dette er ikke offisielt studiemateriell. Les mer.