Tilbake
9.2
Alle fag for Høyskole

9.2 Monte Carlo som inferensverktøy: simulert sannsynlighet, p-verdi og styrke

Estimér det du ikke kan regne: sannsynligheter, forventninger, p-verdier, styrke og feilforplantning ved gjentatte trekk.

65 min
8 oppgaver
Monte Carlo som inferensverktøysimulert sannsynlighetp-verdistyrke
Din fremgang i kapitlet
0 / 8 oppgaver
Kapitlets plass i kurset
Forkunnskaper: Dette kapitlet bygger på kap. 9.1 (inversjonsmetoden, np.random.uniform/normal, np.mean som estimat, np.var(..., ddof=1)), på styrke- og testritualet i kap. 7.3 (kritisk grense, styrke, re-standardisering under μ1\mu_1) og på deltametoden i kap. 5.3 (feilforplantning gjennom en funksjon). Maks/min av flere levetider er fra kap. 4.1.

Sist du var her — de tre resultatene MC-inferensen hviler på (ferdig oppfrisket):

1. np.mean(betingelse) estimerer P(betingelse)P(\text{betingelse}) — andelen True (fra kap. 9.1).
2. Styrke (7.3): for den ensidige z-testen H0:μ=μ0H_0:\mu=\mu_0 mot H1:μ>μ0H_1:\mu>\mu_0 forkaster vi når Z=Xˉ−μ0σ/n>zα\displaystyle Z=\frac{\bar{X}-\mu_0}{\sigma/\sqrt{n}}>z_\alpha; styrken regnes med data trukket under den sanne μ1\mu_1.
3. Deltametoden (5.3): Var(g(θ^))≈[g′(μ)]2Var(θ^)\text{Var}(g(\hat{\theta}))\approx[g'(\mu)]^2\text{Var}(\hat{\theta}) — en tilnærming som kan kontrolleres ved simulering.

Du trenger ikke ny teori — bare å kombinere disse med numpy.

Hverdagsanker: vil du vite hvor ofte to terninger gir sum over 9, kan du enten telle utfall — eller kaste tusen ganger og telle andelen. Den andre veien er Monte Carlo: når regnestykket er tungt eller umulig, simulerer vi mange ganger og leser av andelen. Det virker fordi en sannsynlighet er en forventet indikator, og en forventning kan alltid estimeres med et gjennomsnitt.

Kapitlet lærer deg fire oppskrifter som SKAL sitte: (i) simulert sannsynlighet/forventning for en sammensatt størrelse, (ii) simulert p-verdi, (iii) simulert styrke, og (iv) feilforplantning ved simulering — med deltametoden som fasit å kontrollere mot. Vi starter med grunnprinsippet (løkke 1) og tar så én oppskrift per løkke.

Samlet kjernetid ≈ 65 min. Naturlige pausepunkter er markert mellom løkkene. All kode er komplett og kjørt; Monte Carlo-svar varierer litt med seed, så bruk stor BB (antall simuleringer) og se på om tallet stemmer med den analytiske fasiten.

Løkke 1 — Grunnprinsippet: sannsynlighet som forventet indikator (~12 min)

Grunnprinsippet: sannsynlighet = forventet indikator
La 1(A)\mathbf{1}(A) være indikatoren til en hendelse AA: den er 1 hvis AA inntreffer, ellers 0. Da er forventningen til indikatoren nettopp sannsynligheten:

E[1(A)]=1⋅P(A)+0⋅P(Ac)=P(A).E[\mathbf{1}(A)]=1\cdot P(A)+0\cdot P(A^c)=P(A).

En sannsynlighet er altså en forventning. Og en forventning estimerer vi med et gjennomsnitt. Simulerer vi hendelsen BB ganger, er Monte Carlo-estimatet

p^=1B∑i=1B1(A i trekk i),\hat{p}=\frac{1}{B}\sum_{i=1}^{B}\mathbf{1}(A \text{ i trekk } i),

som i numpy er nøyaktig np.mean(betingelse) — andelen trekk der betingelsen er sann. Samme idé gir en forventning E[g(X)]E[g(X)] som np.mean(g(x)).

MC-estimatet: forventningsrett, varians p(1−p)/Bp(1-p)/B
Antallet «treff» blant BB uavhengige trekk er binomisk, så p^\hat{p} arver egenskapene fra kap. 5.1:

E(p^)=p(forventningsrett),Var(p^)=p(1−p)B.E(\hat{p})=p\quad(\textbf{forventningsrett}),\qquad \text{Var}(\hat{p})=\frac{p(1-p)}{B}.

Standardfeilen er p(1−p)/B\sqrt{p(1-p)/B}, som avtar som 1/B1/\sqrt{B}: firedobler du antall simuleringer, halveres usikkerheten. Det er derfor vi bruker stor BB (ofte 10510^5–10610^6). Merk skillet: BB er antall simuleringer (styrer presisjonen i estimatet), mens en eventuell utvalgsstørrelse nn inne i hver simulering er en helt annen størrelse (den hører til problemet vi simulerer).

✏️Eksempel 1: Simulert sannsynlighet med kontroll

En måling er X∼N(100, 225)X\sim N(100,\,225) (varians 225, altså σ=15\sigma=15). Estimer P(X>120)P(X>120) med Monte Carlo, og sammenlign med fasit.

import numpy as np

def simuler_sanns(B):
    x = np.random.normal(loc=100, scale=15, size=B)   # scale = SD = sqrt(225) = 15
    return np.mean(x > 120)                            # andel over 120 = estimat paa P(X>120)

np.random.seed(0)
print(simuler_sanns(10**6))   # ~ 0.0912

Kjørt gir dette ≈0,0912\approx 0{,}0912. Fasit: P(X>120)=P ⁣(Z>120−10015)=P(Z>1,333)=1−Φ(1,333)≈0,0912\displaystyle P(X>120)=P\!\left(Z>\frac{120-100}{15}\right)=P(Z>1{,}333)=1-\Phi(1{,}333)\approx 0{,}0912. ✓

Usikkerheten: med B=106B=10^6 og p≈0,091p\approx 0{,}091 er standardfeilen 0,091⋅0,909/106≈0,00029\sqrt{0{,}091\cdot 0{,}909/10^6}\approx 0{,}00029 — tre desimaler er trygge. Hadde vi brukt B=1000B=1000, ville SE vært ≈0,009\approx 0{,}009, altfor grovt.

> Sensornotat: np.mean(x > 120) — andel = sannsynlighet. Og scale=15, ikke 225.

📝Oppgave 1

(Innstegsoppgave.) La AA være hendelsen «X>50X>50».

a) Hva er forventningen til indikatoren 1(A)\mathbf{1}(A)?

b) Hvilket numpy-uttrykk estimerer P(A)P(A) fra en vektor x av trekk?

c) Hvordan endres standardfeilen i estimatet om du firedobler antall simuleringer BB?

— naturlig pausepunkt —

Løkke 2 — Oppskrift (i): sannsynlighet/forventning for en sammensatt størrelse (~12 min)

Oppskrift (i): simulert sannsynlighet/forventning

Når størrelsen du spør om er en transformasjon eller kombinasjon av flere variable — maks av levetider, sum, forhold — er det ofte enklere å simulere enn å utlede fordelingen. Oppskriften:

1. Trekk inngangsvariablene vektorisert (én vektor per variabel, lengde BB).
2. Regn den sammensatte størrelsen elementvis (np.maximum, +, /, …).
3. Estimer: np.mean(størrelse > c) for en sannsynlighet, np.mean(størrelse) for en forventning.

For et parallellsystem er levetiden M=max⁡(X1,X2,… )M=\max(X_1,X_2,\dots) (systemet lever til den siste komponenten dør — kobling til kap. 4.1); for et seriesystem er den min⁡(… )\min(\dots). Slik unngår du å regne ut Fmax⁡=[F]nF_{\max}=[F]^n for hånd når du bare trenger et tall.

✏️Eksempel 2: Hvilken komponent varer lengst (parallellsystem)

To uavhengige komponenter har eksponensialfordelte levetider med forventning β1=4\beta_1=4 og β2=6\beta_2=6 år. De er parallellkoblet, så systemet lever til den siste dør: M=max⁡(T1,T2)M=\max(T_1,T_2). Estimer P(M>8)P(M>8) og E(M)E(M), og kontroller mot fasit.

import numpy as np

def simuler_maks(B):
    t1 = np.random.exponential(scale=4, size=B)   # eksponensial, forventning beta=4
    t2 = np.random.exponential(scale=6, size=B)   # forventning beta=6
    m = np.maximum(t1, t2)                         # parallell = maks
    return m

np.random.seed(0)
m = simuler_maks(10**6)
print(np.mean(m > 8))   # ~ 0.363
print(np.mean(m))       # ~ 7.60

(Merk: np.random.exponential(scale=beta) trekker eksponensial med forventning beta — samme fordeling som inversjonsmetoden -beta*np.log(1-u) fra kap. 9.1.)

Kontroll av P(M>8)P(M>8): M≤8M\le 8 betyr at begge er ≤8\le 8, så P(M>8)=1−F1(8)F2(8)=1−(1−e−8/4)(1−e−8/6)=1−0,8647⋅0,7364≈0,3633P(M>8)=1-F_1(8)F_2(8)=1-(1-e^{-8/4})(1-e^{-8/6})=1-0{,}8647\cdot 0{,}7364\approx 0{,}3633. Kjørt gir ≈0,363\approx 0{,}363. ✓

Kontroll av E(M)E(M): for maks av to uavhengige eksponensiale er E(M)=β1+β2−β1β2β1+β2=4+6−2410=7,6\displaystyle E(M)=\beta_1+\beta_2-\frac{\beta_1\beta_2}{\beta_1+\beta_2}=4+6-\frac{24}{10}=7{,}6. Kjørt gir ≈7,60\approx 7{,}60. ✓

> Sensornotat: parallell →max⁡\to\max, serie →min⁡\to\min. Å simulere det tunge E(M)E(M)-integralet på tre kodelinjer er hele poenget.

📝Oppgave 2

To uavhengige komponenter har eksponensielle levetider med forventning β1=3\beta_1=3 og β2=5\beta_2=5. De er seriekoblet, så systemet dør når den første dør: M=min⁡(T1,T2)M=\min(T_1,T_2).

a) Skriv en komplett funksjon simuler_min(B).

b) Hvilket kall estimerer E(M)E(M), og hvilket estimerer P(M>4)P(M>4)?

c) Det er kjent at min⁡\min av to uavhengige eksponensiale igjen er eksponensiell med rate lik summen av ratene. Regn E(M)E(M) og P(M>4)P(M>4) analytisk. (Rate =1/β=1/\beta.)

— naturlig pausepunkt —

Løkke 3 — Oppskrift (ii): simulert p-verdi (~13 min)

Oppskrift (ii): simulert p-verdi
P-verdien er sannsynligheten, beregnet under H0H_0, for å få en testobservator minst så ekstrem som den observerte. Når nullfordelingen er tung å regne ut, simulerer vi den:

1. Regn den observerte testobservatoren tobst_{\text{obs}} fra dataene.
2. Generer mange datasett under H0H_0 (parametrene satt til nullverdiene) og regn testobservatoren i hvert.
3. P-verdien er andelen simulerte observatorer minst så ekstreme som tobst_{\text{obs}}:
p^=np.mean(t_sim >= t_obs)(ensidig oppover).\hat{p}=\texttt{np.mean(t\_sim >= t\_obs)}\quad(\text{ensidig oppover}).

Retningen på «mer ekstrem» må stemme med H1H_1: ensidig oppover bruker >=, ensidig nedover <=, og tosidig bruker np.mean(np.abs(t_sim) >= abs(t_obs)). Nøkkelkravet: generer under H0H_0 — bruker du feil fordeling, får du feil p-verdi.

✏️Eksempel 3: Simulert p-verdi for en z-test

Vi tester H0:μ=100H_0:\mu=100 mot H1:μ>100H_1:\mu>100 med kjent σ=15\sigma=15 og n=20n=20 observasjoner. Vi observerte xˉ=106\bar{x}=106, altså zobs=106−10015/20=1,789\displaystyle z_{\text{obs}}=\frac{106-100}{15/\sqrt{20}}=1{,}789. Estimer p-verdien ved simulering og sammenlign med fasit.

import numpy as np

def simuler_pverdi(B):
    n, mu0, sigma, z_obs = 20, 100, 15, 1.789
    x = np.random.normal(loc=mu0, scale=sigma, size=(B, n))  # generer UNDER H0
    xbar = x.mean(axis=1)
    z = (xbar - mu0) / (sigma / np.sqrt(n))                  # testobservator i hvert datasett
    return np.mean(z >= z_obs)                               # andel minst saa ekstrem

np.random.seed(0)
print(simuler_pverdi(200000))   # ~ 0.037

Kjørt gir ≈0,037\approx 0{,}037. Fasit: siden zz er standardnormal under H0H_0, er p=P(Z≥1,789)=1−Φ(1,789)≈0,0368p=P(Z\ge 1{,}789)=1-\Phi(1{,}789)\approx 0{,}0368. ✓

Her kunne vi lest p-verdien rett av tabellen, og det er nettopp derfor eksemplet er lærerikt: simuleringen treffer fasiten. I en ekte eksamensoppgave er testobservatoren så sammensatt at tabellen ikke rekker til — da er dette den eneste veien.

> Sensornotat: matrisen size=(B, n) gir BB datasett à nn observasjoner; x.mean(axis=1) gir de BB gjennomsnittene. Alt under μ0\mu_0, fordi p-verdien måles under H0H_0.

📝Oppgave 3

Vi tester tosidig H0:μ=100H_0:\mu=100 mot H1:μ≠100H_1:\mu\ne 100, kjent σ=15\sigma=15, n=20n=20, og observerer zobs=2,10z_{\text{obs}}=2{,}10.

a) Hvilken fordeling skal datasettene genereres under?

b) Endre siste linje i funksjonen fra Eksempel 3 slik at den gir den tosidige p-verdien.

c) Regn den tosidige fasiten analytisk. (Bruk Φ(2,10)=0,9821\Phi(2{,}10)=0{,}9821.)

— naturlig pausepunkt —

Løkke 4 — Oppskrift (iii): simulert styrke (~14 min)

Oppskrift (iii): simulert styrke

Styrken er sannsynligheten for å forkaste H0H_0 når alternativet er sant (kap. 7.3). Simuleringen speiler definisjonen:

1. Generer mange datasett under alternativet — sett den sanne parameteren til μ1\mu_1 (ikke μ0\mu_0!).
2. Utfør testen i hvert datasett: regn testobservatoren og avgjør om den havner i forkastningsområdet (f.eks. z>zαz>z_\alpha).
3. Styrken er andelen datasett der vi forkaster: np.mean(forkast).

Den kritiske forskjellen fra p-verdien: p-verdien genererer under H0H_0; styrken genererer under alternativet μ1\mu_1. Bytter du om de to, får du enten α\alpha (om du feilaktig trekker under H0H_0) i stedet for styrken, eller motsatt. Når antall simuleringer BB vokser, konvergerer den simulerte styrken mot den analytiske fasiten fra kap. 7.3.

✏️Eksempel 4: Simulert styrke (samme case som kap. 7.3)

En z-test har μ0=200\mu_0=200, kjent σ=20\sigma=20, n=25n=25, α=0,05\alpha=0{,}05 ensidig oppover, så vi forkaster når z>z0,05=1,645z>z_{0{,}05}=1{,}645. Estimer styrken ved den sanne verdien μ1=210\mu_1=210, og sammenlign med den analytiske fasiten.

import numpy as np

def simuler_styrke(B):
    n, mu0, mu1, sigma, z_alpha = 25, 200, 210, 20, 1.645
    x = np.random.normal(loc=mu1, scale=sigma, size=(B, n))   # generer UNDER alternativet mu1
    z = (x.mean(axis=1) - mu0) / (sigma / np.sqrt(n))         # testobservator (mot mu0)
    return np.mean(z > z_alpha)                               # andel forkastninger = styrke

np.random.seed(0)
print(simuler_styrke(200000))   # ~ 0.804

Kjørt gir ≈0,804\approx 0{,}804. Fasit fra kap. 7.3: kritisk grense k=200+1,645⋅20/25=206,58k=200+1{,}645\cdot 20/\sqrt{25}=206{,}58, og styrke =1−Φ ⁣(1,645−210−20020/25)=1−Φ(−0,855)=Φ(0,855)≈0,804\displaystyle =1-\Phi\!\left(1{,}645-\frac{210-200}{20/\sqrt{25}}\right)=1-\Phi(-0{,}855)=\Phi(0{,}855)\approx 0{,}804. ✓

Legg merke til de to ulike rollene til μ\mu: dataene trekkes med loc=mu1 (sannheten), men testobservatoren standardiseres mot mu0 (påstanden vi tester). Det er nettopp re-standardiseringen fra kap. 7.3, nå gjort automatisk av simuleringen.

> Sensornotat: trekk under μ1\mu_1, test mot μ0\mu_0. Trekker du under μ0\mu_0, får du bare α=0,05\alpha=0{,}05 tilbake.

📝Oppgave 4

Fyll ut de to manglende linjene i funksjonen under, som skal estimere styrken til en ensidig z-test på α=0,05\alpha=0{,}05 (z0,05=1,645z_{0{,}05}=1{,}645) når μ0=500\mu_0=500, kjent σ=40\sigma=40, n=16n=16, og den sanne forventningen er μ1=520\mu_1=520.

def simuler_styrke(B):
    n, mu0, mu1, sigma, z_alpha = 16, 500, 520, 40, 1.645
    x = ______                                    # (1) generer datasettene
    z = (x.mean(axis=1) - mu0) / (sigma/np.sqrt(n))
    return ______                                 # (2) andel forkastninger

a) Fyll ut linje (1) og (2).

b) Forklar med én setning hvorfor svaret nærmer seg den analytiske fasiten når BB øker.

c) Hva ville gått galt om linje (1) hadde loc=mu0?

— naturlig pausepunkt —

Løkke 5 — Oppskrift (iv): feilforplantning ved simulering (~14 min)

Oppskrift (iv): feilforplantning ved simulering

Deltametoden (kap. 5.3) gir en tilnærmet varians for en avledet størrelse g(⋅)g(\cdot) via deriverte. Simulering gir det samme — og fungerer som kontroll:

1. Trekk inngangsvariablene fra fordelingene sine (bruk scale=varians\sqrt{\text{varians}} for normale innganger).
2. Regn g(⋅)g(\cdot) per trekk, elementvis.
3. Estimer variansen med np.var(g_verdier, ddof=1) og standardfeilen som kvadratroten.

Stemmer den simulerte variansen med delta-formelen [g′(μ)]2Var(θ^)[g'(\mu)]^2\text{Var}(\hat\theta) (eller summen ∑gi′2Var\sum g_i'^2\text{Var} for flere innganger), har du kontrollert begge. Avviker de mye, er enten delta-tilnærmingen dårlig (sterkt krum gg eller stor usikkerhet) eller så er det en regnefeil. ddof=1 er obligatorisk her — det er en empirisk varians.

✏️Eksempel 5: Areal fra to målinger — simulering mot delta

Et panel har uavhengige, målte sider XX (forventning 88, varians 0,250{,}25) og YY (forventning 55, varians 0,160{,}16), begge normalfordelte. Arealet er A=XYA=XY. Deltametoden ga Var(A)≈16,49\text{Var}(A)\approx 16{,}49 (kap. 5.3). Kontroller med simulering.

import numpy as np

def simuler_areal(B):
    x = np.random.normal(8, np.sqrt(0.25), size=B)   # scale = SD = sqrt(0.25) = 0.5
    y = np.random.normal(5, np.sqrt(0.16), size=B)   # scale = sqrt(0.16) = 0.4
    a = x * y                                         # g(x,y) = xy per trekk
    return np.mean(a), np.var(a, ddof=1)              # forventning og empirisk varians

np.random.seed(0)
middel, varians = simuler_areal(10**6)
print(middel)    # ~ 40.0
print(varians)   # ~ 16.5

Kjørt gir aˉ≈40,0\bar{a}\approx 40{,}0 og empirisk varians ≈16,5\approx 16{,}5. Delta-fasit: Var(A)≈μY2Var(X)+μX2Var(Y)=52⋅0,25+82⋅0,16=6,25+10,24=16,49\text{Var}(A)\approx \mu_Y^2\text{Var}(X)+\mu_X^2\text{Var}(Y)=5^2\cdot 0{,}25+8^2\cdot 0{,}16=6{,}25+10{,}24=16{,}49. Simuleringen (≈16,5\approx 16{,}5) bekrefter delta-svaret. ✓

> Sensornotat: pass på scale=np.sqrt(varians) (her 0,5 og 0,4), og ddof=1 i variansen. Uten ddof=1 blir avviket ubetydelig ved stor BB, men vanen skal sitte — det er en empirisk varians.

📝Oppgave 5

En strømningsrate er forholdet R=VTR=\dfrac{V}{T} mellom uavhengige, normalfordelte målinger: volum VV (forventning 100100, varians 44) og tid TT (forventning 2020, varians 11). Deltametoden ga Var(R)≈0,0725\text{Var}(R)\approx 0{,}0725 (kap. 5.3).

a) Skriv en funksjon simuler_R(B) som returnerer forventning og empirisk varians til R^\hat{R}.

b) Hvilke scale-verdier skal inn i de to np.random.normal-kallene?

c) Omtrent hvilket tall forventer du at np.var(..., ddof=1) gir, og hva bekrefter det?

Begrepsbank

Kjernekortene for Monte Carlo-inferens, samlet.

Flashcard-stoff — hopp trygt over ved førstegangslesing; tidsanslaget gjelder kjernestoffet.

Kort: np.mean(betingelse) som p^\hat{p}

En sannsynlighet er en forventet indikator: E[1(A)]=P(A)E[\mathbf{1}(A)]=P(A). Derfor er p^=1B∑1(A)=\displaystyle \hat{p}=\frac{1}{B}\sum\mathbf{1}(A)= np.mean(betingelse) et estimat på P(A)P(A). En forventning E[g(X)]E[g(X)] estimeres tilsvarende med np.mean(g(x)).

Kort: presisjon 1/B1/\sqrt{B} og valg av BB
p^\hat{p} er forventningsrett med Var(p^)=p(1−p)/B\text{Var}(\hat{p})=p(1-p)/B, så standardfeilen p(1−p)/B\sqrt{p(1-p)/B} avtar som 1/B1/\sqrt{B}. Firedobling av BB halverer usikkerheten. Velg BB stort (10510^5–10610^6) når du trenger flere sikre desimaler; MC-svar varierer med seed.
Feilkort: H0H_0 (p-verdi) vs. alternativ (styrke)
P-verdi: generer under H0H_0, andel minst så ekstrem som observert. Styrke: generer under alternativet μ1\mu_1, andel forkastninger. Å bytte de to er den vanligste feilen: trekker du styrken under H0H_0, får du bare α\alpha igjen.
Feilkort: BB (simuleringer) vs. nn (utvalg)
BB = antall Monte Carlo-simuleringer, styrer presisjonen i estimatet (1/B1/\sqrt{B}). nn = utvalgsstørrelsen inne i hver simulering, hører til problemet (f.eks. testens nn). Å øke BB gjør estimatet skarpere; å endre nn endrer hva du estimerer (f.eks. styrken selv).
Repetisjonsoppgaver
Din fremgang
0 / 3 oppgaver
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 Norges teknisk-naturvitenskapelige universitet. Dette er ikke offisielt studiemateriell. Les mer.