Download Report

Lunds univrsitet
Matematikcentrum
Matematisk statistik
Biostatistisk grundkurs, MASB11
Laboration 3 VT-2015, 150217
Fördelningsanpassning och Centrala Gränsvärdes Satsen
Introduktion
Syftet med laborationen är dels att vi skall bekanta oss med normalfördelningen, s.v. X ∈ N (µ, σ 2 ),
vad som sker med sannolikhetsfördelning när man bildar summor av slumpvariabler (Centrala
gränsvärdessatsen) och dels med några metoder och funktioner som finns i R för att kunna skatta,
studera och jämföra fördelningen hos ett eller flera stickprov, empirisk fördelningsfunktion, F ∗ (x)
och olika typer av fördelningsplottar. Datamaterialen och script finns på kursens hemsida under
namnet www.maths.lth.se/matstat/kurser/masb11/vtlp3/
• Laborationen redovisas med en skriftlig rapport som ska vara inlämnad senast 150303. Rapporten baseras på korta svar från respektive Uppgift 1-9, vid behov kan också lämplig figur
bifogas. Förberedelseuppgifterna ska vara handskrivna och kan lämnas in tillsammans med
rapporten.
Förberedelseuppgifter
1. Serum kolesterol hos amerikanska kvinnor i åldrarna 20 år och uppåt är (enligt National
Health and Nutrition Examination Survey) normalfördelad med väntevärde 206 mg/dl och
standardavvikelse 44.7 mg/dl.
(a) Hur stor andel av de amerikanska kvinnor som är 20 år eller äldre har en serum kolesterol
som understiger 230 mg/dl?
(b) Hur stor andel av de amerikanska kvinnor som är 20 är eller äldre har en serum kolesterol
mellan 150 och 250 mg/dl?
2. Vid läkemedelstillverkning är det viktigt att mängden aktiv substans hålls så konstant som
möjligt. En viss variation är emellertid omöjlig att undvika, oftast kan den beskrivas med hjälp
av en normalfördelning. För ett visst läkemedel gäller att mängden aktiv substans (mg) i en
tablett antas vara N 2, 0.12 . Mängden aktiv substans i en tablett påverkas inte av de andra
tabletternas tillverkning. Tabletterna ordineras patienterna i förpackningar som innehåller
20 tabletter.
(a) Vad är sannolikheten av mängden aktiv substans i en tablett i förpackningen understiger
1.975 mg?
(b) Låt Xi vara mängden aktiv substans i tablett nr i. Hur kan då Y = ”mängden aktiv
substans i hela förpackningen” tecknas?
(c) Vad är fördelningen för Y ?
(d) Vad är sannolikheten av mängden aktiv substans i hela förpackningen understiger 39.5 mg?
Simulering av slumpvariabler i R
Simulering i R görs genom färdiga funktioner unika för respektive fördelning. Exempel på några av
dessa funktioner finns i efterföljande tabell.
Fördelning
Binomial
Poisson
Normal
Likformig (Rektangel)
Exponential
Funktion i R
rbinom(antal,n,p)
rpois(antal,mean)
rnorm(antal,mean,stddev)
runif(antal,min,max)
rexp(antal,scale) (yr)
Exempel i R
rbinom(100,10,0.5)
rpois(100,2)
rnorm(100,10,15)
runif(100,10,20)
rexp(100,1)
mean = väntevärdet, µ, i fördelningen
stddev = standardavvikelsen, σ, i fördelningen
scale = 1/väntevärdet i exponentialfördelningen
1
Normalfördelningen
Om man till exempel vill ha ett antal slumptal (50 stycken) från en normalfördelning med väntevärdet
(populationsmedelvärdet) 100 och standardavvikelsen 15 skriver man följande kommando:
a1 <- rnorm(50,100,15)
I variabeln a1 lagras då de genererade värdena. Vill man kolla på värdena skriver man bara a1 så
får man en lista på värdena. Vill man beräkna basstatistiken för datamaterialet kan man skriva:
a1sum <- summary(a1)
En boxplot får man genoma att skriva boxplot(a1) och ett histogram genom att skriva hist(a1).
Prova dessa kommandon så du ser vad som händer. Vill du veta mer om ett kommando kan du
skriva kommandot med ett frågetecken framför t.ex: ?hist, ?boxplot, etz.
Uppgift 1
Skapa nu ett stickprov normal om 10 slumptal från en normalfördelning med väntevärdet 10 och
standardavvikelsen 2. Bilda sedan ett andra stickprov som heter uni och innehåller 10 slumptal från
en likformig fördelning mellan 10 och 20. Ledning: rnorm(10,10,2) och runif(10,10,20).
Observera att slumptalen kan ses som stickprov om 10 observationer från två kända populationer. Kontrollera nu med hjälp av kommandot hist(normal) och hist(uni) hur väl stickproven
överensstämmer med populationerna. För ett stort stickprov bör de se ut som den teoretiska
fördelningen tex så här:
Hur väl stämmer stickproven överens med populationerna? Upprepa nu förfarandet för stickprovsstorlek n=50 observationer samt n=500. Bör överensstämmelsen bli bättre eller sämre?
Uppgift 2
Skapa nu tre stickprov om n=1000 observationer från följande fördelningar (kalla dem t.ex. norm1,
uni1 resp. exp1):
• Normal(10, 2).
• Rektangelfördelning – Uniform(10, 20).
• Exponential med väntevärdet 1 (rexp(1000,1)).
Kontrollera med histogrammet hur fördelningarna ser ut.
2
1500
600
500
1000
400
300
500
200
100
0
0
5
10
15
0
10
20
12
14
16
18
20
Figur 1: Normal(10,4) samt Uniform(10,20).
2
QQ-plot och qqnorm
Om man vill kontrollera hur pass nära ett stickprov är en viss teoretisk fördelning kan man använda
olika grafiska metoder. En sådan metod är en s.k. Q-Q plot (Q=Quantile). I en Q-Q plot jämför man
de verkliga värdena i stickprovet med det man kunde förvänta sig från en viss teoretisk fördelning.
Om de observerade värdena överensstämmer med de förväntade så kommer punkterna i en Q-Q
plot att följa en rät linje.
Jämför nu de tre stickproven ovan med vad vi kunde förvänta oss från en normalfördelning. I R finns
det en standardfunktion qqnorm(stickprovsnamn) där man jämför kvantilerna i ett stickprov
med normalfördelningen. I kommandofönstret:
qqnorm(norm1)
qqnorm(uni1)
qqnorm(exp1)
Dina figurer bör se ut ungefär som i Figur 2 nedan. Notera att olika avvikelser från normalfördelning
resulterar i olika former på kurvan.
Uppgift 3
I de tidigare uppgifterna har vi simulerat vad som händer om vi tar stickprov av olika storlekar
från olika kända fördelningar. Vi skall nu gå vidare och undersöka vad som händer om vi bildar
olika storheter i stickprovet. Vilka egenskaper får då dessa storheter? Det är framförallt väntevärdet
E[X], variansen V [X] och fördelningen FX (x) som vi intresserar oss för.
Vi börjar med att undersöka vilken fördelningen summan av två observationer från en normalfördelning
med väntevärde 10 och standardavvikelse 2 har. Börja med generera två nya stickprov om 1000
observationer norm1 och norm2:
norm1 <- rnorm(1000,10,2)
norm2 <- rnorm(1000,10,2)
När vi kör dessa kommandon kommer det att bildas två nya variabler som heter norm1 och norm2
och som innehåller 1000 slumptal var.
Bilda nu summan (sum12) av de två kolumnerna norm1 och norm2. Undersök vilken fördelning
summan har genom att göra ett histogram och en Q-Q plot.
sum12 <- norm1+norm2
Vilken
fördelning har summan? Vad bör väntevärdet bli? Standardavvikelsen? (Använd gärna R,
√
x beräknas med sqrt(x)) Beräkna också medelvärdet av sum12 med mean(sum12) och stickprovsstandardavvikelsen sd(sum12). Hur passar de med de teoretiska värdena?
3
Normal Q−Q Plot − uni1
18
16
14
12
Sample Quantiles
12
10
8
10
6
Sample Quantiles
14
20
Normal Q−Q Plot − norm1
−2
−1
0
1
2
−2
Theoretical Quantiles
−1
0
1
2
Theoretical Quantiles
5
4
3
2
0
1
Sample Quantiles
6
7
Normal Q−Q Plot − exp1
−2
−1
0
1
2
Theoretical Quantiles
Figur 2: QQ-plot för Normal-, Likformig- samt Exp-fördelning.
Centrala gränsvärdessatsen
Lägger man ihop, adderar, (eller beräknar medelvärdet) av flera oberoende normalfördelade slumpvariabler är summan också normalfördelad. Men vad händer om man lägger ihop flera variabler som
alla är rektangelfördelade? Vilken fördelning fås om man adderar exponentialfördelade variabler?
Centrala gränsvärdessatsen säger att om man adderar ett stort antal oberoende variabler från en
godtycklig fördelning blir summan (eller medelvärdet) normalfördelad. Detta märkliga faktum ska
du i denna uppgift undersöka med hjälp av den interaktiva rutinen cgs().
Konkret kan vi tänka oss att du gör ett antal mätningar av en intressant (bio)variabel, du bildar
summan av mätningarna (eller medelvärdet). Det du ska undersöka är hur summan kommer att
variera från mätserie till mätserie? Beror det på ursprungsfördelningen hos den uppmätta variabeln?
Så här använder du rutinen cgs i RStudio
• När du skriver cgs() får du möjlighet att välja mellan ett antal fördelningar med givna parametrar eller kan du konstruera en egen diskret sannolikhetsfördelning. Välj ett av alternativen
4
genom att mata in tillhörande siffra.
• Du får en figur med täthetsfunktion eller sannolikhetsfunktion för din valda fördelning. Välj
nu hur många mätningar du ska göra från denna fördelning och mata in detta antal.
• I kommandofönstret visas resultatet av din mätningar (de 10 första om du valt ett stort antal),
d.v.s. R har hämtat slumptal från din valda fördelning. Summan av alla mätningarna skrivs
ut. I din figur markeras mätningarna med kryss.
• Antag nu att du gör upprepade serier med det antal mätningar, n, som du valt. För varje serie
beräknas summan av dina mätningar. Hur varierar då summan? Mera matematiskt beskrivet:
Om X1 , X2 , . . . , Xn är oberoende med den fördelning du valt, vad är då fördelningen för
summan X1 + X2 + . . . + Xn ?
• Undersök detta genom att simulera N serier med det antal mätningar (n) du valt. Rutinen
ritar sedan upp ett histogram för summan. Ange alltså ett värde på N, tänk på att välja N
tillräckligt stort så att du kan få en uppfattning av fördelningen i histogrammet.
• Centrala gränsvärdessatsen säger att om du valt ett tillräckligt stort antal mätningar kommer
fördelningen för summan att bli ungefär normalfördelning. Rutinen ger dig möjlighet att
anpassa en normalfördelning till data. Du kan låta R sköta om det och din uppgift blir då att
undersöka grafiskt om du tycker att approximationen verkar bra. Till din hjälp har du också
en Q-Q plot där du kan se om summan verkar passa till en normalfördelning.
• Du kan också anpassa normalfördelningen själv och måste då fundera på vilka värden på
väntevärde och standardavvikelse som gäller (prova gärna detta själv som en utmärkt övning!).
• Om du vill köra rutinen igen kan du undvika den interaktiva fasen genom att direkt skriva in
dina val i anropet. Exempelvis ger cgs(2,10,1000,1) att 1000 serier med vardera 10 mätningar
slumpas från en likformig fördelning, R(0, 4). Histogramet för de 1000 summorna plottas,
normalfördelning anpassas och en Q-Q plot ritas.
Uppgift 4
1. Välj rektangelfördelning, antal=2 i R-funktionen cgs(). Vilka värden kan summan av två
mätningar ligga mellan? Verkar histogrammet rimligt?
2. Öka antalet mätningar i rektangelfördelningen. Vad händer om du tar antalet mätningar till
5? Eller ökar till 10?
3. Försök anpassa ”rätt” normalfördelning till histogrammet, d.v.s. tänk ut värdena på väntevärde
och standardavvikelse.
(b−a)2
Ledning: Om s.v. X är uniform(a,b): ⇒ E[X] = a+b
2 och att V [X] =
12 .
4. Exponentialfördelning: Gör nu motsvarande för exponentialfördelningen. Hur många mätningar
behöver ni ta innan ni tycker att summan är ungefär normalfördelat? Verkar fördelningen
gå snabbare eller långsammare mot en normalfördelning än det gjorde för den likformiga
fördelningen. Vad beror detta på?
5. Normalfördelning: Vad händer om ni tar antal=2? Kan du förklara detta?
6. Undersök gärna på motsvarande sätt vad som händer då man bildar summor från binomialeller poissonfördelningen.
7. Testa gärna med en egen diskret fördelning, tex binomialfördelningen.
8. Du har tittat på vad som händer med summor av variabler. Vad händer om man i stället tar
medelvärdet av variablerna (mätningarna)?
5
3
Fördelningsanpassning
3.1
Empiriska fördelningsfunktionen F ∗ (x) – normalfördelningspapper
Grafiska metoder används främst för tre ändamäl: skattning av parametrar, validering av fördelning
samt skattning av kvantiler. Den grafiska tekniken bygger kort på att, vid givet slumpmässigt
stickprov x1 , x2 ,...,xn :
1. Först ordnas stickprovet, betecknas x(1) , x(2) ,...,x(n) .
2. Man skattar fördelningsfunktionen F (x) med det vi kallar för den empiriska fördelningsfunktionen F ∗ (x). Den definieras som:

; x < x(1)
 0
F ∗ (x) =
i/n ; x(i) ≤ x < x(i+1)

1
; x(n) ≤ x
3. Därefter plottas de n stycken talparen (x(i) , ( ni )). Plottningspositionen ni som vi använder för
den empiriska fördelningsfunktionen har en del fördelar men också vissa nackdelar, tex. att
x(n) kommer att vara den position som svarar mot 1 hos fördelningsfunktionen. Andra val av
i
) eller (x(i) , ( i−1/2
plottningspositioner förekommer: exempelvis, (x(i) , ( n+1
n ), se Holmquist B.
Matematisk statistik för M och V, Kompletteringar och tillämpningar, 1996.
I R kan F ∗ (x) ritas med hjälp av funktionen plot(...,type="s",...). Nedansåtende kommandorader exemplifierar teckniken med hjälp av 100 observationer från en s.v. X ∈ N (2, 1).
>
>
>
>
>
X<-rnorm(100,2,1)
sortX<-sort(X)
Fn<-seq(1,length(sortX),1)/(length(sortX)+1) # plotposition foer
plot(sortX,Fn,type="s",col="blue")
grid()
F_n
• På y-axeln har vi F ∗ (x). Använd denna för att skatta medelvärdet, kvartilerna samt medianen
i fördelningen.
• Eftersom vi känner µ och σ i det här fallet kan vi komplettera figuren med den riktiga fördelningsfunktionen, FX (x). Gör det, glöm inte att använda points(...) istället för
plot(...) innan du plottar ovanpå F ∗ (x).
3.2
Normplot i R-package nsRFA
Om vi vet eller misstänker att stickprovet kommer från en normalfördelning kan vi istället plotta det
ordnade stickprovet i ett normalfördelningspapper. Skalan på y-axeln i ett normalfördelningspapper
är anpassad så att observationerna kommer att följa en rät linje om de är normalfördelade. Jämför
teknken med qqnorm. Om vi får någon kurvatur indikerar detta alltså att observationerna inte är normalfördelade. Om man har installerat ett packages som heter nsRFA i R kan man direkt plotta ett
eller flera stickprov i ett normalfördelningspapper med kommandot normplot() och normpoints().
Använd ?normplot för att komma underfund med funktionen. Plotta därefter stickprovet X i ett
normalfördelningspapper.
> ?normplot
> normplot(X)
Uppgift 5
Skatta nu medelvärdet µ och standardavvikelsen σ i normalfördelningsplotten, skattningstekniken
är en direkt tillämpning av moment 2 i Lektionsblad 4. Stämmer skattningarna med det använda
6
stickprovet?
Om R-packages nsRFA inte redan finns installerat i din dator kan man enkelt göra det i RStudio
genom att välja flicken packages i nedre högra fönstret och därefter kryssa för de packages man
vill installera. Det finns en hel del olika fördelnnigspapper man kan ha nytta av i nsRFA.
3.3
Jordbävningar
Vi ska nu studera ett datamaterial där data insamlats under perioden den 16 december 1902 t.o.m.
den 4 mars 1977. Det rör sig om tidsintervall, mätt i dagar, mellan kraftiga jordbävningar världen
runt. Jordbävningar med en magnitud på åtminstone 7,5 på Richterskalan finns representerade,
alternativt jordbävningar med över 1000 dödsoffer.
Datamaterialet finns på kursens material-hemsida under namnet Quakes.RData. Läs in filen genom
kommandot load("Quakes.RData"). De numeriska värdena finns lagrade i en vektor med namn
quakeper. Använd length för att finna antalet tidsperioder. Som tidigare ritar vi histogram och
beräknar diverse läges- och spridningsmått:
>
>
>
>
>
>
hist(quakeper,freq=FALSE) # histogrammets totala area blir 1, taethetsfuktion#
m<-mean(quakeper)
med<- median(quakeper)
s<-sd(quakeper)
s2<- var(quakeper)
myrange<-range(quakeper)
Uppgift 6
Använd data och fundera: verkar det troligt att det kan gå längre period än 5 år mellan kraftiga
jordbävningar? Hur ofta inträffar de?
I R finns en del användbara villkorssatser som gör det enkelt att skapa nya vektorer och matriser
med hjälp av lämpliga bivillkor, Man kan alltså på så sätt i en given vektor eller matris finna
element som uppfyller ett aller annat intressant villkor. För att exmpelvis finna de perioder mellan
jordbävningar med längd kortare än 1000 dagar (c:a 3 år) kan man skriva:
> less1000 <- quakeper[quakeper < 1000];
> length(less1000)
Första kommandot skapar en vektor som vi kan kalla vad som helst, tex. less1000. Den innehåller
de element i ursprungsvektorn quakeper vilka uppfyller villkoret. För att få reda på hur många
element som uppfyller villkoret använder vi helt enkelt length(less1000).
Uppgift 7
Vi vill uppskatta sannolikheten för att en period mellan jordbävningar är kortare än 1000 dagar
genom att beräkna motsvarande andel i datamaterialet Vi har i själva verket beräknat täljaren
fall i kommandoserien ovan, och nämnaren ges helt enkelt av length(quakeper). Beräkna nu den
intressanta kvoten och notera ditt svar. Hur stor är sannolikheten att det dröer mer än 200 dagar
mellan två stora jordbäningar?
Anmärkning. Den storhet som beräknades som mean(quakeper) benämnes ibland återkomst-tid
(engelska: return period), beteckna den tex. med Tr , detta är egentligen en skattning av vänteärdet.
Intensiteten av de händelser som studeras kan beräknas som 1/T och studeras ofta i statistisk
riskanalys, brukar betecknas med λ.
7
3.4
Anpassning till standardfördelning
Enligt gängse statistisk teori är tidsavståndet för två händelser som uppträder slumpmässigt i tiden
efter varandra exponentialfördelade, s.v. T ∈ Exp(λ). (Se sid 98 i kurboken). Dess Fördelningsfunktion
kan då skrivas som FT (t) = 1 − e−λ·t , t > 0.
Uppgift 8
Avgör nu om tidsavståndet mellan två efterföljande jordbävningar är exponentialfördelad genom
att jämföra den empiriska fördelningsfunktionen, F ∗ (t) , med den teoretiska, FT (t) R-tips:
>
>
>
>
>
>
>
sortX<-sort(quakeper)
Fn<-seq(1,length(sortX),1)/(length(sortX)+1) # plotposition foer
plot(sortX,Fn,type="s",col="blue")
grid()
taxis<-seg(0.01,max(quakeper,1))
FT<- 1-exp(-\lambda*taxis)
points(taxis,FT,type="l",col="red")
F_n
Uppgift 9
Centrala gränsvärdessatsen i praktiken: På 35 patienter med Hodgkins sjukdom mätte man antalet
T4 celler i blodet (antal/mm3 ). Samtidigt mätte man motsvarande antal hos 35 patienter som
hade andra sjukdomar (Non-Hodgkins). Data ligger i filen Hodgkindata.RData som du hittar på
kursens hemsida. Läs in data via Workspace-fönstrets öppna-ikon. Du har nu fått två nya variabler
Hodgkin och NonHodgkin. Undersök om antalet celler i blodet är normalfördelat för de båda
grupperna.
För att avgöra om två datamaterial, stickprov, är lika ska de ha samma fördelning de ska ju vara
stickprov från samma population. Avgör nu om de två stickproven i datamaterialet Hodgkindata.RData
har samma fördelning med hjälp av dess empiriska fördelningsfunktioner och normalfördelningspapperet,
normplot(). Tips små storheters.v., X, som är positiva men med värden nära noll kan ofta modelleras med√hjälp av en lognormalfördelning, s.v. log X ∈ N (µ, σ 2 ). Ett annat alternativ är att
använda en X transformation. Stämmer detta?
Det är också möjligt jämföra grupperna genom att bilda differensen mellan de två gruppmedelvärdena. Kan du använda dig av centrala gränsvärdessatsen i detta fall? Kan du säga något
om vilken fördelning differensen i medelvärden har? är det ett stort problem att variabeln inte
är normalfördelad i de båda grupperna från början? Kan man åtgärda detta på något sätt? Vad
händer med fördelningen
√ för stickprovsmedelvärdet om man istället använder en transformation av
värdena, till exempel X eller log X.
8