7.1 Inversjonsmetoden og numpy-verktøykassen

Sist du var her — de tre tingene du må ha friskt:

1. Fordelingsfunksjonen $F(x) = P(X \le x)$ — voksende, fra 0 til 1.
2. Den inverse $F^{-1}$: løser du $F(x) = u$ for $x$, får du $x = F^{-1}(u)$. For en voksende, kontinuerlig $F$ finnes den entydig for $0 < u < 1$.
3. Eksponensialfordelingen (forventnings-parametrisering, forventning $\mu$): tetthet $\displaystyle f(x) = \frac{1}{\mu}e^{-x/\mu}$, fordelingsfunksjon $F(x) = 1 - e^{-x/\mu}$ for $x > 0$. (Rate-parametriseringen $f(x) = \lambda e^{-\lambda x}$, $F(x) = 1 - e^{-\lambda x}$ bruker vi der oppgaven gjør det — les alltid oppgavens egen definisjon.)

Tenk deg at du vil vite hvor ofte en komponent med tilfeldig levetid varer i mer enn 8 år, eller hva forventet verdi av en komplisert funksjon av flere variable er — men integralet er stygt eller finnes ikke i lukket form. I stedet for å regne kan du la datamaskinen «kaste terning» én million ganger og telle. Det er hele idéen bak simulering: en sannsynlighet er en andel, og en forventning er et gjennomsnitt, så trekker vi mange nok realisasjoner, leser vi svaret rett ut av tallene.

Men datamaskinen kan i utgangspunktet bare trekke uniforme tall mellom 0 og 1 (np.random.uniform). Spørsmålet blir da: hvordan lager vi tall fra vår fordeling — eksponensial, Pareto, Rayleigh — ut fra uniforme tall? Svaret er inversjonsmetoden, og den er kapitlets første halvdel. Andre halvdel er numpy-verktøykassen: de få funksjonene og idiomene sensor forventer å se, og hvordan du leser en simuleringssnutt.

Kapitlet er bygget som fire korte læringsløkker — teori → gjennomregnet eksempel → øvingsoppgave — så du kan ta én om gangen. Med fire løkker à ~13 min passer det med en kort pause etter løkke 2.

Løkke 1 — Inversjonsmetoden (~15 min)

$$P(X \le x) = P\big(F^{-1}(U) \le x\big) = P\big(U \le F(x)\big) = F(x).$$

Intuisjon steg 1: Vi bruker at $F$ er voksende, så $F^{-1}(U) \le x$ er det samme som $U \le F(x)$ (vi «kjører $F$ på begge sider»).

Intuisjon steg 2: Siden $U$ er uniform på $(0,1)$, er $P(U \le c) = c$ for enhver $c$ mellom 0 og 1 — og $F(x)$ ligger nettopp mellom 0 og 1. Derfor blir $P(U \le F(x)) = F(x)$.

Intuisjon steg 3: Resultatet $P(X \le x) = F(x)$ er per definisjon fordelingsfunksjonen til $X$. Ferdig — $X$ har fordelingen vi ville ha.

Arbeidsflyten blir dermed alltid den samme:

1. Skaff fordelingsfunksjonen $F(x)$ (utledet i kap. 2.3).
2. Løs likningen $F(x) = u$ for $x$ — det gir formelen $x = F^{-1}(u)$.
3. Implementer: trekk uniforme $u$-er og sett dem inn i formelen.

Løkke 2 — numpy-verktøykassen (~15 min)

— naturlig pausepunkt etter denne løkka —

Inversjonsmetoden lager tallene; numpy og scipy.stats regner på dem. Her er hele verktøykassen sensor forventer — det er ikke mange funksjoner:

- np.random.uniform(size=n) — trekker $n$ uniforme tall på $(0,1)$. Dette er råstoffet til inversjonsmetoden.
- stats.norm.rvs(mu, sigma, size=n) — trekker $n$ normalfordelte tall. Advarsel: annet argument er STANDARDAVVIKET $\sigma$, ikke variansen. Skriver du variansen der, får du feil spredning (en klassisk sensorfelle). For $X \sim N(\mu, \sigma^2)$ med $\sigma^2 = 9$ skriver du stats.norm.rvs(mu, 3, size=n) — altså 3, ikke 9.
- np.mean(x) — gjennomsnittet, estimatet av $E(X)$.
- np.var(x, ddof=1) og np.std(x, ddof=1) — empirisk varians og standardavvik. ddof=1 betyr «del på $n-1$», som gir den forventningsrette variansen (kobling kap. 5.1). Uten ddof=1 deler numpy på $n$ og gir en litt for liten varians.
- Vektorisering. Skriv -mu * np.log(1 - u) på hele vektoren u på én gang — ikke en for-løkke over enkeltelementer. Vektorisert kode er kortere, raskere og det sensor forventer.

Funksjonsmønsteret som går igjen:

def simuler(n):
    u = np.random.uniform(size=n)
    return ...   # inversjonsformelen anvendt vektorisert på u

Løkke 3 — Haler og flere fordelinger (~15 min)

Inversjonsmetoden virker for enhver fordeling der du klarer å løse $F(x) = u$. De «eksotiske» fordelingene fra kap. 2.3 gir alle pene inversjonsformler:

- Pareto (hale, minsteverdi $k$, formparameter $\theta$): $F(x) = 1 - (k/x)^{\theta}$ for $x > k$, som gir $x = k(1-u)^{-1/\theta}$.
- Weibull (skala $\beta$, form $\alpha$): $F(x) = 1 - e^{-(x/\beta)^{\alpha}}$, som gir $x = \beta\big(-\ln(1-u)\big)^{1/\alpha}$.
- Rayleigh (skala $\sigma$): $F(x) = 1 - e^{-x^2/(2\sigma^2)}$, som gir $x = \sigma\sqrt{-2\ln(1-u)}$.

Hver av dem er samme tre-stegs oppskrift: skaff $F$, løs for $x$, implementer vektorisert.

Løkke 4 — Kodelesing som ferdighet (~15 min)

Eksamen ber deg ikke bare skrive kode — den ber deg lese en snutt og svare på to spørsmål: (1) hvilken fordeling trekkes det fra, og (2) hva estimeres? Framgangsmåten:

1. Finn inversjonsformelen i koden og kjenn den igjen. -mu*np.log(1-u) er eksponensial; k*(1-u)**(-1/theta) er Pareto; sigma*np.sqrt(-2*np.log(1-u)) er Rayleigh.
2. Les den siste linjen — estimandet. np.mean(x) estimerer $E(X)$; np.var(x, ddof=1) estimerer $\text{Var}(X)$; np.mean(x > a) estimerer $P(X > a)$ (gjennomsnitt av en indikator er en sannsynlighet); np.mean(g(x)) estimerer $E[g(X)]$.

Nøkkelen til punkt 2: en sannsynlighet er en forventet indikator. x > a er en vektor av True/False (1/0), og gjennomsnittet av den er andelen som oppfyller betingelsen — altså et estimat av $P(X > a)$.

Begrepsbank

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