Afsnit ML.6: Statistik
ML.6.1 Simulere stokastiske variable
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å MATLAB 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
MATLAB.
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
gange og se, hvor stor en andel
af disse hvor summen af antal øjne er større end eller lig med 9.
De
kast af den ene terning genereres med
unidrnd(6,1,10000), hvor "6"
angiver, at de mulige udfald er
svarende til en
terning, og "10000" siger, at der skal laves en vektor med 10000 tilfældige
udfald.
>> x=unidrnd(6,1,10000); % 1.terning kastes 10000 gange
>> y=unidrnd(6,1,10000); % 2.terning kastes 10000 gange
>> z=x+y; % sum af antal øjne i de 10000 kast
>> 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
af en cylinder
med radius
og højde
,
.
Radius er målt til
med en standard error på
og højden er målt til
med en standard error på
For at simulere en spredning på det målte volumen
vil vi trække nye
målinger af radius fra en
-fordeling og nye målinger af
højden fra en
-fordeling
>> nSim=10000;
>> vhat=pi*(10^2)*2;
>> rtilde=normrnd(10,0.2,1,nSim);
>> htilde=normrnd(2,0.1,1,nSim);
>> vtilde=pi*rtilde.*rtilde.*htilde
>> s=struct('Skoen',vhat,'StandardError',sqrt(sum((vtilde-vhat).^2)/nSim))
ML.6.2 Fordelingsfunktioner og fraktiler
I MATLAB 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
inv 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.
>> normcdf(1.96,0,1)
ans = 0.9750
>> normcdf(1.96)
ans = 0.9750
>> chi2cdf(5.99,2)
ans = 0.9500
>> chi2inv(0.95,2)
ans = 5.9915
ML.6.3 Beskrivende statistik
For en vektor
med observerede dataværdier beregner man ofte
forskellige beskrivende mål. Tabellen nedenfor gengiver de mest
almindelige.
For at få en direkte fornemmelse af fordelingen kan man lave et
histogram.
I MATLAB gøres dette med funktionen
histogram.
Denne funktion laver en histogramfigur og danner en struktur, hvorfra antallene
i de forskellige intervaller kan findes.
Hvis data ligger i en vektor
bliver det grundlæggende kald til
funktionen
histogram(x). His man ønsker en anden intervalinddeling, end den
histogram laver,
kan dette gøres ved at tilføje
endePkt til kaldet, hvor
endePkt er en
vektor med de ønskede endepunkter for intervallerne.
Som default laver
histogram et antalhistogram. Man kan få et
tæthedshistogram ved at lave tilføjelsen
'Normalization','pdf' til kaldet.
I optællingen af antal observationer i et interval inkluderer
histogram
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:
h=histogram(x,endePkt); antal=h.Values;
hvor antallene ligger i vektoren
antal. Hvis man laver tilføjelsen
'Normalization','pdf' 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
plot
ML.6.4 Goodness of fit test
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
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
>> x =[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=length(x); % antal observationer
>> phat=sum(x)/(50*n); % skøn over p
>> endePkt=linspace(24,44,11)-0.5; % intervalendepunkter i histogram
>> h=histogram(x,endePkt); % figur med antalshistogram
>> antal=h.Values; % tæller antal observationer i hvert interval
% beregne forventede fra binomialfordeling
>> endePkt0=endePkt(2:10);
>> pr0=binocdf(endePkt0,50,phat); % sandsynlighed til venstre for endepunkter
>> pr=[pr0,1]-[0,pr0]; # intervalsandsynligheder
% første interval fra 0 til og med 25, sidste fra 42 til 50
>> ex=n*pr; % forventede
>> [antal;ex]'
ans =
1.0000 0.2865
0 1.1492
6.0000 3.9734
10.0000 10.0288
16.0000 18.3151
22.0000 23.8725
22.0000 21.7748
15.0000 13.5182
7.0000 5.4920
1.0000 1.5895
% slår intervaller sammen for at få forventede >= 5
>> antal1=[sum(antal(1:3)),antal(4:8),antal(9)+antal(10)];
>> ex1=[sum(ex(1:3)),ex(4:8),ex(9)+ex(10)];
>> [antal1;ex1]'
ans =
7.0000 5.4090
10.0000 10.0288
16.0000 18.3151
22.0000 23.8725
22.0000 21.7748
15.0000 13.5182
8.0000 7.0815
>> gTest=2*sum(antal1.*log(antal1./ex1));
>> df=length(antal1)-1-1;
>> pval=1-chi2cdf(gTest,df);
>> s=struct('G',gTest,'Df',df,'Pvaerdi',pval)
s = struct with fields:
G: 1.1580
Df: 5
Pvaerdi: 0.9488
ML.6.5 Test
I dette afsnit nævnes, hvordan nogle af de klassiske statistiske tests
kan laves i MATLAB.
Jeg starter med
-tests for normalfordelte data med et eller to
observationssæt.
Hvis
er en vektor med dataværdier, kan man lave et
-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:
[h,p,ci,stats]=ttest(x,mu0)
Hvis
mu0 udelades, laves et test for, om middelværdien er nul.
Output fra funktionen er:
Værdierne i
stats tilgås ved at skrive
stats.navn.
De tre værdier er
-teststørrelsen i
tstat, frihedsgraderne
for
-testet i
df og empirisk spredning af data i
sd.
Resultatet af testet i
afhænger af et niveau
, der er sat
til
, det vil sige, at
, hvis
-værdien er under
.
Konfidensintervallet i
ci er et
%
konfidensinterval. Hvis man ønsker
i stedet, kan man tilføje
'Alpha',0.01 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.
Hvis man ønsker et test, hvor det kun er store positive værdier, der er
kritiske, skal man tilføje
'Tail','right' i kaldet
til
ttest.
Ovenfor blev
-testet for en enkelt observationsrække omtalt. Lad os nu betragte
-testet for samme middelværdi i to observationsrækker.
Dette svarer til modellen i afsnit
6.1.
Antag, at vektoren
indeholder dataværdierne for den ene
observationsrække og
dataværdierne for den anden
observationsrække. Den relevante kommando er
[h,p,ci,stats]=ttest2(x,y,'Vartype','equal')
som laver
-testet under antagelse om, at varianserne i de to grupper
er ens. Hvis varianserne ikke er ens, skal man ændre
'equal' til
'unequal'. Output 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
Spredningsskøn og antallet af frihedsgrader er for det fælles skøn
for de to grupper.
I tilfældet med to grupper af observationer kan man lave et test for,
at de to varianser er ens ved
[h,p,ci,stats]=vartest2(x,y)
Output ligner output ovenfor, bortset fra at
tstat og
df
er blevet til
fstat,
df1 og
df2.
Skøn og konfidensinterval er for forholdet mellem de to spredninger.
ML.6.6 Regression
Den lineære regressionsmodel fra afsnit
7.1
analyseres i MATLAB ved hjælp af funktionen
fitlm.
Input er en
modelformel på formen
'xt',
hvor
er responsvariablen, og
er den forklarende variabel.
Desuden skal data være en del af input i form af en
tabel (
datatabel). Formelt bliver kaldet af funktionen på formen
fitlm(datatabel,modelformel)
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
ML.2
øverste højre hjørne og nederste delfigur.
>> a=[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=[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=table(a',h','VariableNames',{'afstand','hastighed'});
>> plot(a,h,'o')
>> xlabel('Afstand')
>> ylabel('Hastighed')
>> title('Regression')
% model analyseres og resultat placeres i lmUD
>> lmUD=fitlm(ahdata,'hastighed~afstand')
lmUD =
Linear regression model:
hastighed ~ 1 + afstand
Estimated Coefficients:
Estimate SE tStat pValue
________ ______ ________ __________
(Intercept) -40.784 83.439 -0.48878 0.62983
afstand 454.16 75.237 6.0364 4.4775e-06
Number of observations: 24, Error degrees of freedom: 22
Root Mean Squared Error: 233
R-squared: 0.624, Adjusted R-Squared 0.606
F-statistic vs. constant model: 36.4, p-value = 4.48e-06
>> par=lmUD.Coefficients.Estimate;
>> refline(par(2),par(1))
>> lmUD.DFE % frihedsgrader df(M)
ans = 22
>> lmUD.RMSE % spredningsskøn
ans = 232.9107
>> plot(a,lmUD.Residuals.Raw,'o')
>> title('Residualer'); refline(0,0);
Funktionen
fitlm 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.
>> logFlu=log([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=repelem(["a","M"],25);
>> tid=repelem(["T1","T2","T3","T4","T5","T1","T2","T3","T4","T5"],5);
>> datadrug=table(logFlu',prodrug',tid','VariableNames',{'logFlu','prodrug','tid'});
>> lmUD=fitlm(datadrug,'logFlu~prodrug*tid')
% model hvor hver kombination af prodrug og tid har sin egen middelværdi
Linear regression model:
logFlu ~ 1 + prodrug*tid
Estimated Coefficients:
Estimate SE tStat pValue
_________ _______ ________ __________
(Intercept) 5.3198 0.23878 22.279 3.6821e-24
prodrug_M -1.463 0.33768 -4.3324 9.6452e-05
tid_T2 1.474 0.33768 4.365 8.7212e-05
tid_T3 2.3244 0.33768 6.8834 2.7264e-08
tid_T4 2.4653 0.33768 7.3006 7.1788e-09
tid_T5 2.8768 0.33768 8.5193 1.585e-10
prodrug_M:tid_T2 -0.079595 0.47756 -0.16667 0.86847
prodrug_M:tid_T3 0.30995 0.47756 0.64904 0.52002
prodrug_M:tid_T4 0.26366 0.47756 0.5521 0.58395
prodrug_M:tid_T5 0.48821 0.47756 1.0223 0.31278
Number of observations: 50, Error degrees of freedom: 40
Root Mean Squared Error: 0.534
R-squared: 0.878, Adjusted R-Squared 0.85
F-statistic vs. constant model: 31.9, p-value = 1.58e-15
>> lmUD1=fitlm(datadrug,'logFlu~prodrug+tid')
% model med additiv struktur af prodrug og tid
Linear regression model:
logFlu ~ 1 + prodrug + tid
Estimated Coefficients:
Estimate SE tStat pValue
________ _______ _______ __________
(Intercept) 5.2216 0.1805 28.928 2.861e-30
prodrug_M -1.2665 0.14738 -8.5938 5.758e-11
tid_T2 1.4342 0.23303 6.1546 1.9993e-07
tid_T3 2.4794 0.23303 10.64 9.5266e-14
tid_T4 2.5971 0.23303 11.145 2.1248e-14
tid_T5 3.121 0.23303 13.393 3.984e-17
Number of observations: 50, Error degrees of freedom: 44
Root Mean Squared Error: 0.521
R-squared: 0.872, Adjusted R-Squared 0.857
F-statistic vs. constant model: 60, p-value = 1.57e-18
% Lave F-test for additiv struktur
>> Ftest=((lmUD1.SSE-lmUD.SSE)/(lmUD1.DFE-lmUD.DFE))/lmUD.MSE;
>> pval=1-fcdf(Ftest,lmUD1.DFE-lmUD.DFE,lmUD.DFE);
>> [Ftest,pval]
ans = 0.4764 0.7528
Foregående