UPPSALA UNIVERSITET Matematiska institutionen M˚ans Thulin
Statistik f¨or ingenj¨orer 1MS008 VT 2011
DATOR ¨ OVNING 2: SKATTNINGAR OCH KONFIDENSINTERVALL
1 Inledning
I den h¨ar dator¨ovningen ser vi hur R kan utnyttjas f¨or att kontrollera modellantaganden och ber¨akna konfidensintervall.
2 Att kontrollera modellantaganden
F¨or att kunna genomf¨ora djupare statistisk analys av data s˚a m˚aste vi g¨ora modellan- taganden. Vi ska h¨ar titta p˚a metoder f¨or att grafiskt unders¨oka antagandet att data kommer fr˚an en normalf¨ordelad slumpvariabel.
Vi s˚ag p˚a f¨orel¨asningen att histogram och sannolikhetspapper kan anv¨ands f¨or att unders¨oka normalf¨ordelningsantagandet. F¨or data fr˚an en normalf¨ordelning b¨or histo- grammet likna normalf¨ordelningens klockformade t¨athetsfunktion och p˚a sannolikhets- pappret s˚a b¨or punkterna ligga l¨angs linjen – framf¨orallt vid linjens mitt.
Vi b¨orjar med att titta p˚a histogram och sannolikhetspapper f¨or simulerade data fr˚an normalf¨ordelningen:
par(mfrow=c(1,2)) # G¨or att man f˚ar tv˚a figurer i samma grafikf¨onster y<-rnorm(50,0,1) # Generar 50 N(0,1)-f¨ordelade observationer
hist(y,freq=FALSE) # Ritar histogram
curve(dnorm(x,mean(y),sd(y)),col=2,add=TRUE) # L¨agger till t¨athetsfunktion qqnorm(y);qqline(y) # Ritar sannolikhetspapper
Ovning. Prova att k¨¨ ora koden ovan n˚agra g˚anger f¨or att se hur figurerna varierar.
J¨amf¨or sedan med motsvarande figurer f¨or simulerade data fr˚an exponentialf¨ordelningen:
par(mfrow=c(1,2)) # G¨or att man f˚ar tv˚a figurer i samma grafikf¨onster y<-rexp(50,1) # Generar 50 Exp(1)-f¨ordelade observationer
hist(y,freq=FALSE) # Ritar histogram
curve(dnorm(x,mean(y),sd(y)),col=2,add=TRUE) # L¨agger till t¨athetsfunktion qqnorm(y);qqline(y) # Ritar sannolikhetspapper
Alternativt kan l˚adagram anv¨andas f¨or att unders¨oka f¨ordelningsantagandet. F¨or nor- malf¨ordelningen s˚a b¨or l˚adagrammet vara symmetriskt, medan det f¨or exponentialf¨ordelningen b¨or ha l¨angre morrh˚ar upp˚at ¨an ned˚at:
# L˚adagram f¨or simulerade N(1,1) och Exp(1)-data:
boxplot(rnorm(50,1,1),rexp(50,1),names=c("N(1,1)","Exp(1)")) 2.1 Centrala gr¨ansv¨ardessatsen
Centrala gr¨ansv¨ardessatsen s¨ager att summan av n stycken slumpvariabler ¨ar approxi- mativt normalf¨ordelad om n ¨ar tillr¨ackligt stort. Vi ska nu anv¨anda v˚ara verktyg f¨or att se om data ¨ar normalf¨ordelade f¨or att unders¨oka om satsen st¨ammer genom att titta p˚a summor av exponentialf¨ordelade slumpvariabler. Nedanst˚aende kod ritar histogram f¨or 1000 observationer av olika summor av exponentialf¨ordelade simulerade slumpvariabler.
K¨or koden och ¨oppna sedan grafikf¨onstret i fullsk¨arm.
expsum<-function(n,B) {
summa<-0 for(i in 1:n) {
summa<-summa+rexp(B) }
return((summa-n)/sqrt(n)) }
par(mfrow=c(3,4)) for(i in seq(1,60,5)) {
y<-expsum(i,1000)
hist(y,freq=FALSE,main=paste("Summan av",i,"Exp(1)-variabler")) curve(dnorm(x,0,1),col=2,add=TRUE)
}
P˚aminner histogrammet mer om normalf¨ordelningens t¨athetsfunktion f¨or summor av fler slumpvariabler?
Vi kan ¨aven rita sannolikhetspapper f¨or summor av exponentialf¨ordelade slumpva- riabler:
par(mfrow=c(3,4)) for(i in seq(1,60,5)) {
y<-expsum(i,100)
qqnorm(y,main=paste("Summan av",i,"Exp(1)-variabler")) qqline(y)
}
Ligger punkterna mer l¨angs linjen ju fler slumpvariabler man summerar?
3 Konfidensintervall
3.1 Ett stickprov (passningstider)
Vi ska h¨ar anknyta till exempel 7.1 i l¨aroboken d¨ar 8 observerade passningstider f¨or mobiltelefoner analyserades. Mata in data f¨or hand och lagra i en vektor kallad x:
x <- c(210,214,195,190,218,202,207,197)
Samma f¨oruts¨attningar om normalf¨ordelning etc. g¨ors som i boken.
Ovning. Anv¨¨ and verktygen fr˚an avsnittet ovan f¨or att unders¨oka om det verkar rimligt att data kommer fr˚an en normalf¨ordelning (men kom ih˚ag att det antagandet alltid ¨ar sv˚art att unders¨oka f¨or sm˚a stickprovsstorlekar).
Vi ska nu med hj¨alp av n˚agra f˚a kommandon i R skapa ett 95% konfidensintervall f¨or den genomsnittliga vikten. Vi har ingen kunskap om f¨ordelningens standardavvikelse och r¨aknar d¨arf¨or med ok¨and standardavvikelse. Fr˚an teorin (s. 67-68 i boken) vet vi att intervallet d˚a ges av
h
¯
x ± tα/2(n − 1) s
√n i Det finns tv˚a s¨att att angripa detta p˚a numeriskt:
1. Skriv in intervallet ovan manuellt
2. Anv¨and en f¨ardig rutin (vid namn t.test) F¨or den f¨orsta metoden, skriv
mv <- mean(x) stad <- sd(x) n <- 8
tkvantil <- qt(0.975,n-1)
mv - tkvantil*stad/sqrt(n) # Nedre intervallgr¨ans mv + tkvantil*stad/sqrt(n) # ¨Ovre intervallgr¨ans (J¨amf¨or med svaret i exempel 7.2!)
Som du s¨akert insett ger funktionen qt med l¨ampliga inargument kvantiler till t- f¨ordelningen. Observera att 1 − α/2 = 0.975 d˚a α = 0.05; R anv¨ander allts˚a 1 − α/2 ist¨allet f¨or α/2 som inparameter i kvantilfunktionen.
P˚a liknande s¨att kan kvantiler f¨or andra vanliga f¨ordelningar erh˚allas. F¨or nor- malf¨ordelningens kvantiler anv¨ands funktionen qnorm.
Ovning. Skriv t.ex. in f¨¨ oljande och j¨amf¨or de v¨arden du f˚ar med den vanliga tabellen p˚a sidan 82:
qnorm(0.975) qnorm(0.95)
F¨or den andra metoden anropas helt enkelt rutinen t.test och man f˚ar d˚a ut, bok- stavligen i ett enda slag, f¨orutom sj¨alva konfidensintervallet en m¨angd ytterligare in- formation. Kommando:
t.test(x)
Ar man en van anv¨¨ andare av R och kan sin statistik anv¨ands med f¨ordel metod 2. ¨Ar man nyb¨orjare och ”vill veta vad man g¨or” kan metod 1 vara s¨akrare. Du noterade v¨al att de gav samma svar?
Konfidensgraden 0.95 ¨ar f¨orinst¨alld vid anropet av t.test, men kan ¨andras. Se hj¨alptexten ?t.test f¨or detaljer, speciellt parametern conf.level.
3.2 Tv˚a stickprov
Rutinen t.test kan anv¨andas ¨aven f¨or att ber¨akna konfidensintervall f¨or skillnader i v¨antev¨arden mellan tv˚a stickprov. H¨ar finns, som vi sett i teorin, tv˚a angreppss¨att: tv˚a oberoende stickprov respektive stickprov i par.
3.2.1 Tv˚a oberoende stickprov
Tryckh˚allfastheten f¨or tv˚a olika betongblandningar, av typen M20 respektive M25 ska j¨amf¨oras. Vi l¨aser in data i R:
M20<-c(35.50, 27.80, 35.80, 30.10, 27.60, 32.45, 30.20, 26.85, 31.10, 19.20, 25.86, 31.20, 25.60, 31.15, 35.80, 27.50, 28.73, 23.20, 18.95, 24.50, 22.45, 29.80, 35.65, 30.80, 24.01, 25.25, 27.55, 30.15, 24.50, 22.60)
M25<-c(31.20, 35.86, 31.00, 39.01, 35.60, 38.00, 29.68, 27.26, 30.88, 35.50, 28.88, 38.50, 27.60, 26.00, 37.10, 30.80, 34.45, 38.00, 33.51, 35.80, 31.20, 36.52, 29.82, 37.80, 35.01, 36.60, 32.25, 31.50, 28.65, 27.55)
Ovning. R¨¨ akna ut medelv¨arde f¨or respektive datamaterial och rita l˚adagram med boxplot f¨or att unders¨oka om det verkar finnas n˚agon skillnad mellan v¨antev¨ardet (µM 20 respektive µM 25) f¨or tryckh˚allfastheten f¨or de tv˚a blandningarna.
Ett 99 % konfidensintervall f¨or differensen µM 25− µM 20 ges av t.test(M25,M20,conf.level=0.99)
3.2.2 Stickprov i par
Som illustration anv¨ander vi R f¨or exempel 7.6 (dragstyrka hos metallst¨anger). H¨ar
¨
ar det fr˚aga om modellen stickprov i par, vilket m˚aste anges f¨or R med inparametern paired=TRUE:
xfore <- c(370,360,380,395,375); xefter <- c(400,396,412,420,410);
t.test(xefter,xfore,paired=TRUE) J¨amf¨or med ber¨akningarna i boken.
3.3 *Konfidensintervall f¨or p i Bin(n, p)
I avsnitt 7.3 i boken beskrivs hur man genom normalapproximation kan konstruera konfidensintervall f¨or parametern p i Bin(n, p)-f¨ordelningen. Konfidensintervallet
h ˆ p − λα/2
r1
np(1 − ˆˆ p), p + λˆ α/2 r1
np(1 − ˆˆ p) i
har approximativt konfidensgrad 1 − α.
Konfidensgraden s¨ager oss hur stor andelen f¨ors¨ok som resulterar i konfidensintervall som inneh˚aller det sanna v¨ardet p˚a p ¨ar om vi genomf¨or ett stort antal f¨ors¨ok. Om approximationen ¨ar bra s˚a borde andelen ligga n¨ara 1 − α.
Vi ska h¨ar studera den faktiska konfidensgraden f¨or konfidensintervall f¨or p genom simulering. Koden nedan g¨or att funktionen binKonf(x,n) ger konfidensintervallet f¨or en binomialf¨ordelad observation x och antalet f¨ors¨ok n.
binKonf<-function(x,n) {
p.hatt<-x/n
konf.int<-c(p.hatt-1.96*sqrt(1/n*p.hatt*(1-p.hatt)), p.hatt+1.96*sqrt(1/n*p.hatt*(1-p.hatt))) return(konf.int)
}
Ovning. F¨¨ or att testa funktionen j¨amf¨or vi med exempel 7.4 p˚a s. 69 i boken. D¨ar ¨ar x = 12 och n = 200. Prova att skriva binKonf(12,200) och kontrollera att det ger samma konfidensintervall som i boken.
Du kan sedan provk¨ora funktionen f¨or simulerade data genom att k¨ora koden nedan n˚agra g˚anger. Prova g¨arna att ¨andra v¨ardet p˚a n och p och se hur m˚anga g˚anger som v¨ardet p˚a p ligger i konfidensintervallet.
n<-10; p<-0.5
x<-rbinom(1,n,p) # Ger en Bin(n,p)-f¨ordelad observation.
binKonf(x,n)
Slutligen kan du k¨ora nedanst˚aende kod n˚agra g˚anger f¨or att kolla hur stor andel av de 10000 simulerade konfidensintervallen som inneh˚aller det korrekta v¨ardet p˚a p:
n<-10; p<-0.5
antal<-0 # Variabel f¨or simuleringen for(i in 1:10000)
{
x<-rbinom(1,n,p) # Ger en Bin(n,p)-f¨ordelad observation.
konf<-binKonf(x,n)
# Kolla om p ligger i konfidensintervallet:
if(konf[1]<=p && konf[2]>=p){antal<-antal+1}
}
# Skriv resultatet p˚a sk¨armen:
cat(paste("Andel konfidensintervall som inneh˚aller p:",antal/10000,"\n")) F¨or vilka v¨arden p˚a n och p ligger andelen n¨ara 0.95?
Ovning. Anv¨¨ and koden f¨or att unders¨oka tumregeln ”approximationen ¨ar bra om n · p · (1 − p) ≥ 5”!
4 *Skattningar ¨ ar slumpvariabler
Att en skattning ¨ar en slumpvariabel inneb¨ar att vi kan studera den precis som andra slumpvariabler och ber¨akna exempelvis dess v¨antev¨arde och standardavvikelse. P˚a s˚a vis kan vi teoretiskt j¨amf¨ora olika skattningar f¨or att avg¨ora vilken som ¨ar b¨ast. F¨or det mesta vill vi att skattningen ska vara v¨antev¨ardesriktig (s˚a att den ”i genomsnitt” ger
det r¨atta v¨ardet) och att dess standardavvikelse skall vara s˚a liten som m¨ojligt (s˚a att skattningen f¨orhoppningsvis inte avviker s˚a mycket fr˚an det sanna parameterv¨ardet).
F¨or m˚anga skattningar kan man relativt ”enkelt” r¨akna ut v¨antev¨arde och standard- avvikelse. Exempelvis g¨aller det att om X1, . . . , Xn ¨ar oberoende N (µ, σ2)-f¨ordelade slumpvariabler s˚a ¨ar stickprovsmedelv¨ardet ¯X ∼ N (µ, σ2/n) - denna anv¨ands som bekant f¨or att skatta v¨antev¨ardet µ.
Stickprovsmedelv¨ardet ¨ar inte den enda t¨ankbara skattningen av µ. En normalf¨ordelad slumpvariabel med v¨antev¨arde µ har ¨aven median µ, s˚a en t¨ankbar estimator ¨ar stick- provsmedianen ˆX. F¨ordelningen f¨or ˜X ¨ar, liksom v¨antev¨arde och standardavvikelse, betydligt sv˚arare att r¨akna ut ¨an motsvarande egenskaper f¨or stickprovsmedelv¨ardet.
H¨ar kommer R till v˚ar unds¨attning!
Antag att vi har 10 observationer fr˚an N (µ, 1)-f¨ordelningen. Vi vill veta om den b¨asta skattningen ¨ar stickprovsmedelv¨ardet ¯X eller stickprovsmedianen ˜X. Vi vet fr˚an teo- retiska utr¨akningar att E( ¯X) = µ och att V ( ¯X) = 1/10. Genom att simulera ett antal observationer av ˜X kan vi skatta E( ¯X) och V ( ¯X).
Vi vet inte hur v¨antev¨ardet och variansen f¨or ˜X beror p˚a µ, men vi kan prova att stoppa in olika v¨arden p˚a µ f¨or att se om skattningen ¨ar v¨antev¨ardesriktig och om variansen ¨andras d˚a µ ¨andras. Vi provar att s¨atta µ lika med 0, 1 och 5 och att simulera 1000 stickprov f¨or varje v¨antev¨arde. Vi r¨aknar ut medianen i varje stickprov och g¨or d¨armed 1000 simuleringar vardera av ˜X f¨or de olika v¨ardena p˚a µ. Vi anv¨ander sedan dessa f¨or att skatta E( ˜X) och V ( ˜X) i de olika fallen.
med0<-med1<-med5<-NA # Variabler f¨or simuleringen for(i in 1:1000)
{
med0[i]<-median(rnorm(10,0,1)) med1[i]<-median(rnorm(10,1,1)) med5[i]<-median(rnorm(10,5,1)) }
# Skriv resultatet p˚a sk¨armen:
cat(paste(" my=0:\nE(medianen) =",mean(med0),"\nV(medianen) =",var(med0),"\n",
"my=1:\nE(medianen) =",mean(med1),"\nV(medianen) =",var(med1),"\n",
"my=5:\n--- E(medianen) =",mean(med5),"\nV(medianen) =",var(med5),"\n")) Ligger v¨antev¨ardet n¨ara µ? Beror variansen f¨or medianen p˚a v¨ardet p˚a µ? ¨Ar variansen
mindre ¨an V ( ¯X) = 1/10? Vilken av skattningarna tycker du att man ska anv¨anda?
