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