💻 Simulering i NumPy

Simulering er ein viktig del av læreplanen i matematikk for samfunnsfag (S1 + S2) i LK20. Mange lærebøker er opptatt av løkker i simuleringa si. Det tar tid å køyra, kan vera meir tungt for elevane å forstå og kan by på meir feilmeldingar mtp. syntaks. Her kjem eit (kanskje av fleire) innlegg med tips til simulering med NumPy.

Random generator

For å simulera lagar me ein random generator. Denne kallar me for rng (random number generator). Me kan tilpassa default_rng() med å gje han ulike argument (feks. seed for å få reproduserbare tal) men i praksis i S1 og S2 er det ikkje noko me treng å tenka på.

import numpy as np
rng = np.random.default_rng()

Denne generatoren kan me bruka til mange ulike simuleringar. Om me vil kasta ein terning ein gong kan me trekka eit heiltal. Det gjer me med funksjonen integers(low, high) som ligg i rng-en vår.

terning = rng.integers(1, 7)
print(terning)

Pass på

Her vil terning = rng.integers(1, 7) gje oss

\[\text{terning} \in \{1, 2, 3, 4, 5, 6\}\]

Med andre ord \[1 \leq \text{terning} < 7\]

Om me vil kasta fleire terningar kan me i tillegg gje inn size som argument til rng.integers

terningar = rng.integers(1, 7, 10)
print(terningar)

[1 1 4 2 5 1 6 5 3 5]

Vidare kan me rekna med resultata. Desse er lagra i ein array (ein variabeltype som minnar om lister).

Me kan t.d. sjekka kor mange av dei 10 terningane me kasta som vart seksarar. Først kan me sjekka kva verdiar som er 6:

terningar == 6

array([False, False, False, False, False, False,  True, False, False,
       False])

Dette kan me bruka til å rekna med. Sidan False vert tolka som 0 og True vert tolka som 1 kan me finna summen av arrayen med sanningsvariablane:

sum(terningar == 6)

Døme 🎲

Me trillar terning.

\(A\): terningen viser 5 eller 6.

Finn \(P(A)\) ved hjelp av simulering.

import numpy as np
rng = np.random.default_rng()

# tal simuleringar
N = 10000000

# triller terning N gongar
terningar = rng.integers(1, 7, N)

# finn antall femmarar og seksarar
gunstige = sum(terningar >= 5)

# finn og skriv ut relativ frekvens
rel_frek = gunstige / N
print(f"P(A) = {rel_frek:.4f}")

P(A) = 0.3332

Ut frå 10 000 000 simuleringar ser me at \(P(A)=0,3333\). Dvs. det er \(33,33\%\) sannsynleg å trilla 5 eller 6 på terningen.

Dømet over er eit problem som ikkje er vanskeleg å løysa analytisk, og dermed ikkje veldig spennande å bruka simulering på. Men om problema vert vanskelege å rekna på er dei ofte ikkje så vanskelege å simulera. Eit litt meir avansert problem (som er lett å simulera) finn me i S1-eksamen frå våren 2022:

I eit spel kastar du tre terningar. Du multipliserer saman augetalet på terningane. Dersom dette produktet er større enn 100, vinn du. Bestem sannsynet for å vinne ved å køyre simuleringar. Hugs å vise korleis du kjem fram til svaret. UDIR (2022) - Eksamen S1-V22

For å svare på denne oppgåva treng me litt meir kjennskap til korleis me kan jobba med arrays. Arrays er ein form for lister eller vektorar. Me kan utføra ulike rekneoperasjonar med dei. Me kan t.d. legga dei saman og multiplisera dei:

a = np.array([1, 2, 3, 4, 5, 6])
b = np.array([6, 5, 4, 3, 2, 1])

print(a+b)
print(a*b)

[7 7 7 7 7 7]
[ 6 10 12 12 10  6]

Då har me alt me treng for å løysa eksamensoppgåva:

import numpy as np
rng = np.random.default_rng()

# antall simuleringar
N = 1000000

# triller 3 terningar N gongar
t1 = rng.integers(1, 7, N)
t2 = rng.integers(1, 7, N)
t3 = rng.integers(1, 7, N)

# finn produktet
prod = t1* t2 * t3

# tel opp kor mange av produkta som er over 100
gunstige = sum(prod > 100)

# finn og skriv ut relativ frekvens
rel_frek = gunstige / N
print(f"P(produkt over 100) = {rel_frek:.4f}")

P(produkt over 100) = 0.0928

Det er med andre ord ca. \(9,34\%\) sannsynleg å vinna spelet.

Binomisk sannsynsmodell

I rng har me òg mogleiken for simulering av binomiske fordelingar.

rng.binomial(n, p) der n er talet på delforsøk og p er sannsynet for suksess i kvart delforsøk.

Døme på bruk

Me plantar 150 frø. Sannsynet for at eit frø spirer er 85%.

Kva er sannsynet for at minst 130 frø spirer?

\(X\): antall frø som spirer.

import numpy as np
rng = np.random.default_rng()

N = 1000000     # antall simuleringar

n = 150         # antall delforsøk/frø
p = 0.85        # sannsynet for at eitt frø spirer

# ser kor mange av frøa som spirer kvar simulering
spirer = rng.binomial(n, p, N)

# tell opp dei simuleringane som gjev 130+ frø som spirer
gunstige = sum(spirer >= 130)

rel_frek = gunstige / N
print(f"P(X >= 130) = {rel_frek:.4f}")

P(X >= 130) = 0.3311

Ser at \(P(X\geq 130) = 0,3321\).

Hypergeometrisk sannsynsmodell

På same måte kan me simulera hypergeometrisk fordeling.

rng.hypergeometric(ngood, nbad, nsample) der ngood og nbad er dei to delane av populasjonen vår og nsample er storleiken på utvalet me trekk.

Døme på bruk

Det ligg 4 daim og 7 banantwist i ei skål. Du trekk 3 bitar. Kva er sannsynet for at akkurat to av desse er daim?

import numpy as np
rng = np.random.default_rng()

# antall simuleringar
N = 1000000

ngood = 4      # daim i populasjonen
nbad = 7       # banan i populasjonen
nsample = 3    # antall bitar me trekk

# simulerer
twist = rng.hypergeometric(ngood, nbad, nsample, size=N)

# tel opp kor mange gonagar me fekk 2 daim
gunstige = sum(twist == 2)

rel_frek = gunstige / N
print(f"P(X=2) = {rel_frek:.4f}")

P(X=2) = 0.2539

Det er \(25,5\%\) sannsynleg at det er akkurat to daimtwist.

Oppsummert

Simulering med numpy er altså fleksibelt og raskt med mindre komplisert syntaks enn store løkker.

I tillegg er det mange andre kjekke funksjonar som mellom anna rng.choice() og rng.normal() som kan vera kjekke. Meir om dei ein annan dag 😄

Gjenbruk

CC BY 4.0