I statistik vil vi ofte simulere data for at se, hvordan disse opfører sig under idealiserede forhold. Dette kan være noget så simpelt som at se, hvad der sker, når vi kaster en "fair" terning. Som et eksempel kan vi undersøge, hvor ofte summen af antal øjne på to terninger er større end eller lig med 9. Vi vil derfor gerne kunne få python til at kaste terningerne for os. Man taler om at generere tilfældige tal (random number generator). Følgende tabel viser nogle af mulighederne i python. For binom, poisson, multinomial og norm skal disse importeres via kommandoen from scipy.stats import navn. $$\begin{array}{ll} \hline \text{Navn} & \text{Beskrivelse} \\ \hline \text{np.random.choice(6,10,replace=True)} & \text{10 tilfældige hele tal mellem 0 og 5} \\ \text{binom.rvs(20,0.5,0,10)} & \text{10 tilfældige udfald fra en binomial(20,0.5)} \\ \text{poisson.rvs(3,0,10)} & \text{10 tilfældige udfald fra en poisson(3)} \\ \text{multinomial.rvs(3,[0.1,0.6,0.3],10)} & \text{10 tilfældige udfald fra en} \\ & \text{multinom(3,(0.1,0.6,0.3))} \\ \text{norm.rvs(mu,sigma,10)} & \text{10 tilfældige udfald fra en N(mu,sigma)} \\ \text{np.random.permutation(6)} & \text{tilfældig permutation af tallene 0 til 5} \\ \hline \end{array}$$ Som et første eksempel vil jeg svare på det spørgsmål, der blev stillet ovenfor i forbindelse med kast af to terninger. Jeg vil kaste de to terninger 10$\,$000 gange og se, hvor stor en andel af disse hvor summen af antal øjne er større end eller lig med 9. De 10$\,$000 kast af den ene terning genereres med sample(6,10000,replace=T), hvor "6" angiver, at de mulige udfald er $1,2,3,4,5,6$ svarende til en terning, og "10000" siger, at der skal laves en vektor med 10$\,$000 tilfældige udfald.

import numpy as np

x=1+np.random.choice(6,10000,replace=True) # 1.terning kastes 10000 gange
y=1+np.random.choice(6,10000,replace=True) # 2.terning kastes 10000 gange
z=x+y                                    # sum af antal øjne i de 10000 kast
print(sum(z>=9)/10000)                   # andel med sum >=9

Som et andet eksempel vil jeg finde variansen på en transformeret størrelse ved simulering. Lad os betragte volumet $v$ af en cylinder med radius $r$ og højde $h$, $v=\pi r^2 h$. Radius er målt til $10\,\text{cm}$ med en standard error på $0.2\,\text{cm},$ og højden er målt til $2\,\text{cm}$ med en standard error på $0.1\,\text{cm}.$ For at simulere en spredning på det målte volumen $\hat v=\pi\cdot 10^2\cdot 2\,\text{cm}^3,$ vil vi trække nye målinger af radius fra en $N(10,0.2^2)$-fordeling og nye målinger af højden fra en $N(2,0.1^2)$-fordeling

import numpy as np
from scipy.stats import norm

nSim=10000
vhat=np.pi*(10**2)*2
rtilde=norm.rvs(10,0.2,nSim)
htilde=norm.rvs(2,0.1,nSim)
vtilde=np.pi*rtilde*rtilde*htilde
print('Skøn:',vhat,'   Standard Error:',np.sqrt(sum((vtilde-vhat)**2)/nSim))

I python kan man finde fordelingsfunktionen, det vil sige sandsynligheden for at ligge til venstre for et punkt, for alle standard fordelinger. Disse er navngivet ved, at der sættes .cdf efter navnet på fordelingen. De tilsvarende funktioner for at finde fraktiler i en fordeling fremkommer ved at sætte .ppf efter navnet på fordelingen. Ved en standard normalfordeling kan man udelade mu=0 og sigma=1 i kaldet til norm. Følgende tabel viser nogle af mulighederne for fordelingsfunktionen. For at have fordelingen til rådighed skal denne importeres via kommandoen from scipy.stats import navn. $$\begin{array}{ll} \hline \text{Navn} & \text{Beskrivelse} \\ \hline \text{norm.cdf(1.96,mu,sigma)} & N_{\text{cdf}}(1.96,\text{mu,sigma}) \\ \text{binom.cdf(3,10,0.5)} & P(X\leq 3)\enspace X\sim\text{binomial(10,0.5)} \\ \text{poisson.cdf(3,1.5)} & P(X\leq 3),\enspace X\sim\text{poisson(1.5)} \\ \text{chi2.cdf(5.99,2)} & \chi_{\text{cdf}}(5.99,2) \\ \text{t.cdf(4.3,2)} & t_{\text{cdf}}(4.3,2) \\ \text{f.cdf(2.9,7,12)} & F_{\text{cdf}}(2.9,7,12) \\ \hline \end{array}$$

from scipy.stats import norm, chi2

print(norm.cdf(1.96,0,1))
| 0.9750021048517795

print(norm.cdf(1.96))
| 0.9750021048517795

print(chi2.cdf(5.99,2))
| 0.9499633729134137

print(chi2.ppf(0.95,2)
| 5.991464547107979

For en vektor $x$ med observerede dataværdier beregner man ofte forskellige beskrivende mål. Tabellen nedenfor gengiver de mest almindelige, hvor numpy er importeret somnp. $$\begin{array}{ll} \hline \text{Navn} & \text{Beskrivelse} \\ \hline \text{np.mean(x)} & \text{gennemsnit} \\ \text{np.var(x,ddof=1)} & \text{empirisk varians} \\ \text{np.std(x,ddof=1)} & \text{empirisk spredning} \\ \text{np.median(x)} & \text{median} \\ \hline \end{array}$$ For at få en direkte fornemmelse af fordelingen kan man lave et histogram. I python gøres dette med funktionen plt.hist efter komandoen import matplotlib.pyplot as plt. Denne funktion laver en histogramfigur og danner en struktur, hvorfra antallene i de forskellige intervaller kan findes. Hvis data ligger i en vektor $x$ bliver det grundlæggende kald til funktionen plt.hist(x). Den intervalinddeling, som hist laver, er ikke særlig køn, idet den blot laver 10 intervaller fra den mindste til den største værdi i vektoren $x$. Man kan lave sin egen inddeling ved at tilføje bins=endePkt til kaldet, hvor endePkt er en vektor med endepunkter for intervallerne. Som default laver hist et antalhistogram. Man kan få et tæthedshistogram ved at lave tilføjelsen density=True til kaldet.I optællingen af antal observationer i et interval inkluderer hist det venstre endepunkt, men ikke det højre (pånær det sidste interval, som også inkluderer det højre endepunkt). Dette er modsat af definitionen af et histogram i denne webbog. Hvis vi også vil have antallet af observationer i de forskellige intervaller som output bliver det generelle kald:

antal,bi,pa=plt.hist(x,bins=endePkt)

hvor antallene ligger i vektoren antal. Hvis man laver tilføjelsen density=True indeholder vektoren antal tæthedsværdierne i stedet for antallene. Titel for figuren, og titler på akserne, kan sættes på som ved et almideligt kald af plt.plot Et jordstykke er delt op i 100 lige store områder. I hvert af disse områder er der plantet 50 ærter, og der registreres hvor mange af disse, der spirer og vokser op. De målte antal betegnes $x_i,$ $i=1,\ldots,100.$ Nedenfor laves et histogram med intervallængde 2 og efterfølgende et goodness of fit test for, om data kan beskrives med en binomialfordeling. Metoden er omtalt i afsnit 3.4

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import binom, chi2

x=np.array([31,34,31,35,37,35,37,25,32,35,33,40,31,34,30,37,41,34,36,\
31,28,33,40,28,32,31,39,34,38,34,35,36,36,32,36,35,38,40,\
37,42,38,37,33,34,35,36,37,32,36,39,32,37,28,34,33,31,36,\
41,38,32,33,39,34,37,36,40,32,39,39,38,34,34,38,35,31,33,\
32,28,38,37,34,36,35,37,34,37,34,40,39,38,34,28,30,29,32,\
39,36,33,36,31])
 
n=len(x)  # antal observationer
phat=sum(x)/(50*n)  # skøn over p
endePkt=np.linspace(24,44,11)-0.5 # intervalendepunkter i histogram
antal,bi,pa=plt.hist(x,bins=endePkt)  # figur med antalshistogram

# beregne forventede fra binomialfordeling
endePkt0=endePkt[1:10]
pr0=binom.cdf(endePkt0,50,phat) # sandsynlighed til venstre for endepunkt
pr=np.append(pr0,1)-np.append(0,pr0) # intervalsandsynligheder
# første interval fra 0 til og med 25, sidste fra 42 til 50
ex=n*pr  # forventede

print(np.array([antal,np.around(ex,2)]))
| [[ 1.    0.    6.   10.   16.   22.   22.   15.    7.    1.  ]
|  [ 0.29  1.15  3.97 10.03 18.32 23.87 21.77 13.52  5.49  1.59]]

# slår intervaller sammen for at få forventede >= 5
antal1=np.concatenate(([sum(antal[0:3])],antal[3:8],[antal[8]+antal[9]]))
ex1=np.concatenate(([sum(ex[0:3])],ex[3:8],[ex[8]+ex[9]]))
print(np.array([antal1,np.around(ex1,2)]))
| [[ 7.   10.   16.   22.   22.   15.    8.  ]
|  [ 5.41 10.03 18.32 23.87 21.77 13.52  7.08]]

gTest=2*sum(antal1*np.log(antal1/ex1))
df=len(antal1)-1-1
pval=1-chi2.cdf(gTest,df)

print('G:',gTest,'   Df:',df,'   P-værdi:',pval)
| G: 1.1580228811532804    Df: 5    P-værdi: 0.9488405467978349

I dette afsnit nævnes, hvordan nogle af de klassiske statistiske tests kan laves i python. Jeg starter med $t$-tests for normalfordelte data med et eller to observationssæt. Her er de indbyggede funktioner i scipy.stats eller statsmodels ikke optimale, og jeg har i stedet lavet nogle funktioner, der ligger i filen pytFunktioner.py. Hvis $x$ er en vektor med dataværdier, kan man lave et $t$-test for, at middelværdien har en bestemt værdi mu0 med ttest. Dette er gennemgået i afsnit 4.7. Funktionen kaldes på følgende vis:

ttest(x,mu=mu0)

Hvis mu=mu0 udelades, laves et test for, om middelværdien er nul. Output fra ttest er en pandas dataframe. Hvis output placeres i testUD kan man få de forskellige dele af output ved at skrive testUD['navn'][0], hvor mulighederne er: $$\begin{array}{ll} \hline \text{Navn} & \text{Beskrivelse} \\ \hline \text{tstat} & \text{værdien af t-teststørrelsen} \\ \text{df} & \text{antallet af frihedsgrader} \\ \text{p} & p\text{-værdien} \\ \text{est} & \text{skøn over middelværdien} \\ \text{lower} & \text{nedre grænse for 95\%-konfidensinterval} \\ \text{upper} & \text{øvre grænse for 95\%-konfidensinterval} \\ \hline \end{array}$$ Konfidensintervallet er for den ukendte middelværdi, og størrelsen af konfidensintervallet bestemmes af et niveau, der kan vælges ved at tilføje conflevel=0.95 i kaldet til ttest. Funktionen ttest laver et såkaldt tosidet test, det vil sige, at både store positive og store negative værdier er kritiske. Ovenfor blev $t$-testet for en enkelt observationsrække omtalt. Lad os nu betragte $t$-testet for samme middelværdi i to observationsrækker. Dette svarer til modellen i afsnit 6.1. Antag, at vektoren $x$ indeholder dataværdierne for den ene observationsrække og $y$ dataværdierne for den anden observationsrække. Den relevante kommando er

ttest2(x,y,varequal=True)

som laver $t$-testet under antagelse om, at varianserne i de to grupper er ens. Hvis varianserne ikke er ens, skal man ændre True til False. Outputtet fra funktionen indeholder de samme elementer som ovenfor. Skøn og konfidensintervallet er i dette tilfælde for forskel i middelværdi mellem de to grupper $\delta=\mu_1-\mu_2.$I tilfældet med to grupper af observationer kan man lave et test for, at de to varianser er ens ved

vartest2(x,y)

Output ligner output ovenfor, borset fra at tstat og df er blevet til fstat, df1 og df2. Skøn og konfidensinterval er for forholdet mellem de to spredninger. Den lineære regressionsmodel fra afsnit 7.1 analyseres i python ved hjælp af funktionen ols fra modulet statsmodels.formula.api. Input er en modelformel på formen 'x$\sim$t', hvor $x$ er responsvariablen, og $t$ er den forklarende variabel. Desuden skal data være en del af input i form af en dataframe. Formelt bliver kaldet af funktionen på formen

ols(formula=modelformel,data=dataframe).fit()

Den følgende kode viser et eksempel med Hubbles data omtalt i afsnit 7.6. De to figurer der produceres, kan ses i den sidste figur i afsnit Py.2 øverste højre hjørne og nederste delfigur. En meget brugbar efterbehandling af resultatet fra ols er funktionen summary2. Denne giver tre tabeller hvor det især er den miderste parametertabel der er af interesse, hvor der er skøn og konfidensintervaller for parametrene såvel som et $t$-test for, om parameteren kan antages at være nul. I den første tabel optræder Scale, som indeholder variansskønnet $s^2_r$ fra afsnit 7.2.

from statsmodels.formula.api import ols
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pytFunktioner import *

a=np.array([0.032,0.034,0.214,0.263,0.275,0.275,0.450,0.500,\
0.500,0.630,0.800,0.900,0.900,0.900,0.900,1.000,\
1.100,1.100,1.400,1.700,2.000,2.000,2.000,2.000])
h=np.array([170,290,-130,-70,-185,-220,200,290,270,200,300,\
-30,650,150,500,920,450,500,500,960,500,850,800,1090])
ahdata=pd.DataFrame({'a':a,'h':h})

plt.plot(a,h,'o')
plt.xlabel('Afstand')
plt.ylabel('Hastighed')
plt.title('Regression')
plt.show()

# model analyseres og resultat placeres i lmUD
lmUD=ols(formula='h~a',data=ahdata).fit()

summaryLM(lmUD)  # parameterestimater

| Estimated Coefficients:
|              Coef.  Std.Err.       t  P>|t|    [0.025   0.975]
| Intercept -40.7836   83.4389 -0.4888 0.6298 -213.8253 132.2580
| a         454.1584   75.2371  6.0364 0.0000  298.1262 610.1906
|  
| Number of observations: 24  Error degrees of freedom: 22
| Root Mean Squared Error: 232.9
| R-squared: 0.624  Adjusted R-Squared: 0.606
| F-statistic vs. constant model: 36.4  p-value = 4.48e-06

par=lmUD.params

plt.refline(par[1],par[0])  # linje indtegnes

print({'s(M)':np.sqrt(lmUD.mse_resid)})
| {'s(M)': 232.9106701830066}

plt.figure()
plt.plot(a,lmUD.resid,'o')
plt.title('Residualer')
plt.axhline(0)  # nullinje indtegnes
plt.show()

Funktionen ols kan også bruges til den generelle lineære model beskrevet i kapitel 8 og kapitel 9. Følgende eksempel ser på data i afsnit 8.6. Variable med tekststrenge bliver automatisk opfattet som faktorer.

from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
import numpy as np
import pandas as pd
from scipy.stats import f
from pytFunktioner import *

logFlu=np.log(np.array([303,335,86,168,243,1323,568,1165,669,966,1724,3339,
3259,1731,1224,3140,2276,2861,2077,1893,2969,3000,5285,3536,3781,
132,23,105,62,12,229,264,119,81,434,917,266,782,1291,506,622,883,
532,1388,493,1459,1168,898,1774,1772]))
prodrug=np.repeat(["a","M"],25)
tid=np.repeat(["T1","T2","T3","T4","T5","T1","T2","T3","T4","T5"],5)
datadrug=pd.DataFrame({'logFlu':logFlu,'prodrug':prodrug,'tid':tid})

# tosidet variansanalyse
lmUD=ols(data=datadrug,formula='logFlu~prodrug*tid').fit()
summaryLM(lmUD) 
 
| Estimated Coefficients:
|                          Coef.  Std.Err.       t  P>|t|  [0.025  0.975]
| Intercept               3.8569    0.2388 16.1524 0.0000  3.3743  4.3394
| prodrug[T.a]            1.4630    0.3377  4.3324 0.0001  0.7805  2.1455
| tid[T.T2]               1.3944    0.3377  4.1293 0.0002  0.7119  2.0769
| tid[T.T3]               2.6344    0.3377  7.8013 0.0000  1.9519  3.3169
| tid[T.T4]               2.7289    0.3377  8.0814 0.0000  2.0465  3.4114
| tid[T.T5]               3.3651    0.3377  9.9651 0.0000  2.6826  4.0475
| prodrug[T.a]:tid[T.T2]  0.0796    0.4776  0.1667 0.8685 -0.8856  1.0448
| prodrug[T.a]:tid[T.T3] -0.3100    0.4776 -0.6490 0.5200 -1.2751  0.6552
| prodrug[T.a]:tid[T.T4] -0.2637    0.4776 -0.5521 0.5840 -1.2288  0.7015
| prodrug[T.a]:tid[T.T5] -0.4882    0.4776 -1.0223 0.3128 -1.4534  0.4770 
|  
| Number of observations: 50  Error degrees of freedom: 40
| Root Mean Squared Error: 0.5339
| R-squared: 0.878  Adjusted R-Squared: 0.85
| F-statistic vs. constant model: 31.9  p-value = 1.58e-15

# additive model
lmUD1=ols(data=datadrug,formula='logFlu~prodrug+tid').fit()
summaryLM(lmUD1)

| Estimated Coefficients:
|               Coef.  Std.Err.       t  P>|t|  [0.025  0.975]
| Intercept    3.9551    0.1805 21.9115 0.0000  3.5913  4.3189
| prodrug[T.a] 1.2665    0.1474  8.5938 0.0000  0.9695  1.5636
| tid[T.T2]    1.4342    0.2330  6.1546 0.0000  0.9646  1.9038
| tid[T.T3]    2.4794    0.2330 10.6399 0.0000  2.0098  2.9490
| tid[T.T4]    2.5971    0.2330 11.1451 0.0000  2.1275  3.0668
| tid[T.T5]    3.1210    0.2330 13.3931 0.0000  2.6513  3.5906
|   
| Number of observations: 50  Error degrees of freedom: 44
| Root Mean Squared Error: 0.5211
| R-squared: 0.872  Adjusted R-Squared: 0.857
| F-statistic vs. constant model: 60.0  p-value = 1.57e-18

# Lave F-test for additiv struktur
df=lmUD.df_resid; SSD=lmUD.mse_resid*df
df1=lmUD1.df_resid; SSD1=lmUD1.mse_resid*df1
Ftest=((SSD1-SSD)/(df1-df))/(SSD/df)
pval=1-f.cdf(Ftest,df1-df,df)
print([Ftest,pval])

| [0.47644440231630847, 0.7527549308670061]

Foregående