1. No category

¨
L UNDS TEKNISKA H OGSKOLA
M ATEMATIKCENTRUM
M ATEMATISK STATISTIK
D ATORLABORATIONER – DEL II, HT-14
¨ BIO -, KEMI - OCH NANOTEKNIK ; FMS086 & MASB02
M ATEMATISK STATISTIK F OR
Datafiler samt en del special skrivna Matlab-kommandon h¨amtas fr˚an kursens hemsida:
www.maths.lth.se/matstat/kurser/fms086
6
Skattningar av parametrarna i en normalf¨ordelning
¨
Ovningens
syfte:
• Du ska unders¨oka hur skattningar av v¨antev¨arde och varians beror av stickprovsstorleken.
Kurskompendium:
Olbjer kap 7.1
F¨orberedelseuppgifter:
Givet a¨r ett stickprov x1 , . . . , xn fr˚an N(μ, σ2 ) d¨ar μ och σ2 a¨r ok¨anda.
1. Hur skattas μ?
2. Vilken f¨ordelning har skattningen av μ?
3. Hur skattas σ2 ?
Utg˚a fr˚an en normalf¨ordelning N(3, 4). Antag att v¨antev¨ardet 3 och variansen 4 a¨r ok¨anda f¨or oss och
att viP
vill skatta dem genom att ta ett stickprov, x1 , . . . , xn , om n observationer och bilda ¯x respektive s2 =
n
1
x )2 . Hur n¨ara kommer skattningarna de sanna v¨ardena om stickprovsstorleken a¨r 5? om
i=1 (xi − ¯
n−1
den a¨r 25?
• Simulera 1000 stickprov om 5 v¨arden fr˚an N(3, 4) och skatta μ i varje stickprov. G¨or samma sak f¨or
1000 stickprov som alla best˚ar av 25 observationer. Unders¨ok v¨antev¨ardesskattningarna egenskaper,
t ex genom att g¨ora histogram.
>>
>>
>>
>>
>>
>>
>>
>>
X=normrnd(3,2,5,1000);
Z=normrnd(3,2,25,1000);
mx=mean(X);
mz=mean(Z);
subplot(2,1,1)
hist(mx,20)
subplot(2,1,2)
hist(mz,20)
Vann man mycket p˚a att o¨ ka stickprovsstorleken fr˚an 5 observationer till 25? Hur stor a¨r sannolikheten att skattningen avviker mer a¨n 2 enheter fr˚an det sanna v¨ardet μ=3? d˚a du anv¨ander 5
v¨arden i stickprovet, 25 v¨arden i stickprovet? (Anv¨and normcdf eller normspec, men t¨ank f¨orst ut
vilken f¨ordelning skattningen har d˚a n=5? d˚a n=25?)
¨
• G¨or a¨ven histogram f¨or σ2 -skattningarna och jfr de tv˚a stickprovsstorlekarna n=5 och n=25. Ar
det ovanligt att skattningen av σ2 avviker mer a¨n 2 fr˚an det sanna v¨ardet 4 (dvs understiger 2 eller
o¨ verstiger 6) d˚a n=5? d˚a n=25?
8
7
HYPOTESTEST – ANALYS AV KEMILABORATION
Tolkning av konfidensintervall f¨or v¨antev¨ardet
¨
Ovningens
syfte:
• Du ska, genom att simulera konfidensintervall f¨or v¨antev¨ardet μ f˚a en tolkning av begreppet konfidensintervall och speciellt unders¨oka hur intervallen beror av stickprovsstorleken.
Kurskompendium:
Olbjer kap 7.3
F¨orberedelseuppgifter:
Givet a¨r ett stickprov x1 , . . . , xn fr˚an N(μ, σ2 ) d¨ar μ och σ2 a¨r ok¨anda.
1. Hur bildas ett 95% konfidensintervall f¨or μ?
F¨or att g¨ora konfidensintervall f¨or v¨antev¨ardet beh¨ovs kvantilfunktionen f¨or t-f¨ordelningen som i Matlab
heter tinv. Definitionen av kvantilfunktionen a¨r h¨ar det tal som a¨r s˚adant att sannolikheten att f˚a ett utfall
mindre a¨n talet a¨r lika med argumentet. Vill man ha t0.025 (n − 1) t ex, skriver man tinv(0.975,n-1).
F¨ordelen med Matlabs definition a¨r att kvantilfunktionen blir samma som inversen till f¨ordelningsfunktionen. F¨ordelen med att ha definitionen tv¨artom, som i kurskompendiet, a¨r att argumentet blir lika med
risken vid ett ensidigt test att f¨orkasta nollhypotesen trots att den a¨r sann.
• Anv¨and de 2000 simulerade stickproven fr˚an f¨oreg˚aende o¨ vning f¨or att g¨ora 2000 95% konfidensintervall f¨or μ. H¨alften av dem kommer d˚a att vara baserade p˚a 5 observationer medan de o¨ vriga p˚a
25.
Plotta g¨arna ut, i samma diagram men med olika symboler, o¨ vre och undre gr¨anserna i konfidensintervallet f¨or n˚agra stickprov (f¨orslagsvis 150 st).
Matlabtips: Om o¨ vre och undre gr¨anserna i konfidensintervallen ligger i variablerna over respektive
under kan de 150 f¨orsta intervallen illustreras med f¨oljande sekvens (det a¨r de lodr¨ata avst˚andet
mellan ”stj¨arna” och ”ring” som utg¨or konfidensintervallet)
>>
>>
>>
>>
plot(over(1:150),’*’)
hold
plot(under(1:150),’o’)
line([0 150], [3 3])
Alternativt kan du anv¨anda funktionen plotint.
• Hur m˚anga intervall missar det sanna v¨ardet μ=3, d˚a n=5? d˚a n=25? (Ska det vara n˚agon skillnad?)
Hur m˚anga intervall ska enligt teorin i genomsnitt missa μ?
• Hur p˚averkar stickprovsstorleken konfidensintervallen?
8
Hypotestest – Analys av kemilaboration
¨
Ovningens
syfte:
¨
Ovningen
ska illustrera
• hur olika test kan ber¨aknas med hj¨alp av Matlab,
• hur de tre olika metoderna f¨or test a¨r sammankopplade,
• hur viktigt det a¨r att ha klart f¨or sig vilka modellantaganden man g¨or om data n¨ar man anv¨ander ett
test i Matlab,
• vilka slutsatser man kan dra fr˚an testet.
12
8
HYPOTESTEST – ANALYS AV KEMILABORATION
Kurskompendium:
Olbjer kap 7.5
F¨orberedelseuppgifter:
1. Givet a¨r ett stickprov x1 , . . . , xn fr˚an N(μ, σ2 ) d¨ar μ och σ2 a¨r ok¨anda. Hur testas H0 : μ = μ0 mot
H1 : μ 6= μ0
(a) med ett signifikanstest, d¨ar signifikansniv˚an a¨r fix?
(b) med ett signifikanstest d¨ar man har r¨aknat ut ett P-v¨arde?
(c) med hj¨alp av konfidensintervall?
2. Givet a¨r ett stickprov x1 , . . . , xn1 fr˚an N(μ1 , σ2 ) och ett stickprov y1 , . . . , yn2 fr˚an N(μ2 , σ2 ) d¨ar
μ1 , μ2 och σ2 a¨r ok¨anda. Hur kan du unders¨oka om μ1 och μ2 skiljer sig a˚t?
Nu skall ni anv¨anda m¨atv¨arden fr˚an kemilaborationen ”Best¨amning av koppar med atomabsorptionsspektrofotometri med grafitugn”. De som har gjort laben kan anv¨anda sina egna m¨atv¨arden, o¨ vriga kan
l˚ana v¨arden av n˚agon som gjort laben eller be labhandledaren om anvisning.
Antag h¨ar att det inte finns n˚agon os¨akerhet i den kalibreringskurva ni anv¨ant f¨or att f˚a fram kopparkoncentrationerna.
Om du inte har ”¨oversatt” de tv˚a uppm¨atta dataserierna till kopparkoncentrationer med hj¨alp av kalibreringskurvan kan du g¨ora det med f¨oljande kommandon:
L¨agg in de uppm¨atta absorbtionerna i varsin vektor abs3 och abs10
>> abs3 = [ ... ]
>> abs10 = [ ... ]
% Uppm¨
atta absorbtioner i 3-serien
% Uppm¨
atta absorbtioner i 10-serien
L¨agg in v¨ardena fr˚an kalibreringskurvan, k¨anda kopparkoncentrationer i en kolonnvektor kalibx och
motsvarande absorbtioner i kaliby
>> kalibx = [0 0 50 50 100 100 150 150 200 200]’
>> kaliby = [ ... ]’
% eller de v¨
arden ni anv¨
ant
Kalibreringskurvan, y = α + βx kan skattas med (detaljer kommer senare i kursen)
>>
>>
>>
>>
n = length(kaliby)
b = regress(kaliby, [ones(n,1) kalibx])
alpha = b(1)
beta = b(2)
sedan o¨ vers¨atts absorbtionerna till kopparkoncentrationer med x =
y−α
β
>> cu3 = (abs3 - alpha)/beta
>> cu10 = (abs10 - alpha)/beta
I detta avsnitt skall vi inte ta h¨ansyn till os¨akerheten i kalibreringskurvan utan vi antar att v¨ardena i cu3
och cu10 a¨r uppm¨atta direkt.
• Anv¨and kommandot ttest, eller ttest2 (anv¨and hj¨alptexterna f¨or att avg¨ora vilken test som a¨r
l¨amplig), f¨or att g¨ora 95% konfidensintervall f¨or kopparkoncentrationen baserat dels p˚a 3-serien
och p˚a 10-serien. Skriv ned resultaten s˚a att ni kan j¨amf¨ora dem med resultaten d˚a vi tar h¨ansyn
till os¨akerheten i kalibreringskurvan senare i laborationen. Om ni nu hade f˚att reda p˚a den sanna
kopparkoncentrationen skulle ni kunna anv¨anda intervallen f¨or att testa om ni har n˚agot systematiskt fel i era m¨atningar. Hur g¨or ni detta?
13
9
STYRKEFUNKTION
• Testa om det a¨r n˚agon skillnad mellan v¨antev¨arde f¨or de 10 m¨atningarna och v¨antev¨ardet f¨or de 3
m¨atningarna.
• Sl˚a upp formeln f¨or konfidensintervall f¨or σ2 baserat p˚a ett stickprov fr˚an N (μ, σ2 ) i formelsamlingen och ber¨akna ett 95% konfidensintervall f¨or σ2 (eller σ genom att dra roten ur gr¨anserna) dels
baserat p˚a 10-serien och dels p˚a 3-serien. Blir det stor skillnad mellan intervallbredderna?
9
Styrkefunktion
¨
Ovningens
syfte:
• Att illustrera begreppet styrkefunktion hos ett signifikanstest;
• att illustrera vilka faktorer som p˚averkar styrkefunktionen.
Kurskompendium:
Olbjer kap 7.6
F¨orberedelseuppgifter:
1. Givet a¨r ett stickprov x1 , . . . , xn fr˚an N(μ, σ2 ) d¨ar μ och σ2 a¨r ok¨anda. Du testar H0 : μ = μ0 mot
H1 : μ 6= μ0 . Hur definieras styrkefunktionen f¨or testet ovan? Vad kan styrkefunktionen anv¨andas
till?
Antag att du har n observationer fr˚an N(μ, σ2 ) och vill testa att μ=6. F¨or att g¨ora det konkret, anta att
du vid ett laboratorieexperiment vill pr¨ova om pH-v¨ardet p˚a en l¨osning kan vara 6 genom att g¨ora n
best¨amningar. Antag vidare, att man gjort upprepade best¨amningar tidigare med samma instrument och
d¨arf¨or anser att man k¨anner dess variation p˚a denna typ av l¨osning och att σ2 a¨r 0.6.
Vi har allts˚a ett stickprov x1 , . . . , xn fr˚an N(μ, 0.6) och vill testa
H0 : μ = 6
H1 : μ 6= 6
p˚a signifikansniv˚a α. Hur ”bra” a¨r detta test? Intressanta fr˚agor kan t ex vara:
• Om det sanna pH-v¨ardet inte a¨r 6 utan 5.5 (avvikelsen a¨r −0.5 fr˚an nollhypotesens v¨arde), med
vilken sannolikhet kommer jag d˚a att uppt¨acka att H0 a¨r falsk med detta test?
• Hur m˚anga best¨amningar m˚aste jag g¨ora f¨or att med sannolikheten 0.90 uppt¨acka att H0 a¨r falsk d˚a
μ i sj¨alva verket a¨r 7 (dvs avvikelsen fr˚an nollhypotesens v¨arde a¨r 1)?
Denna typ av fr˚agor kan besvaras med hj¨alp av testets styrkefunktion som definieras som
π(μ) = P(H0 f¨orkastas || det sanna parameterv¨ardet a¨r μ)
Just i det h¨ar fallet, d˚a vi har ett tv˚asidigt test med k¨ant σ, blir styrkefunktionen (se avsnitt 7.6.2 i Olbjer)
√
√
(μ − 6) n
(μ − 6) n
π(μ) = Φ(−zα/2 −
) + 1 − Φ(zα/2 −
)
σ
σ
dvs den beror p˚a α – testets signifikansniv˚a, σ2 – f¨ors¨oksfelvariansen (som vi antar k¨and) samt n – stickprovsstorleken.
14
10
¨ REGRESSION
ENKEL LINJAR
• Anv¨and funktionen styrka (se help styrka eller stencilen om Matlabkommandon) f¨or att se
hur styrkefunktionen ser ut d˚a testets signifikansniv˚a α a¨r 0.05, σ2 =0.6 och stickprovsstorleken n a¨r
9.
Observera att i figuren som f˚as fr˚an styrka ritas p˚a x-axeln avvikelsen fr˚an nollhypotesens μ0 , dvs
funktionen som ritas a¨r
P(H0 f¨orkastas || μ avviker fr˚an μ0 med c)
¨ du n¨ojd med funktionen? Hur stor a¨r slh att f¨orkasta H0 om μ a¨r 6.5? Hur stor a¨r slh att uppt¨acka
Ar
att H0 inte a¨r sann, d˚a μ a¨r 5?
• Hur skulle en ”ideal styrkefunktion” se ut i det h¨ar exemplet? Skissera den p˚a papper! Du vill
naturligtvis att slh att f¨orkasta H0 ska vara liten om μ verkligen a¨r 6, men att slh ska vara stor s˚a fort
μ avviker fr˚an 6 (dvs om H0 inte a¨r sann).
F¨or att f¨orb¨attra styrkefunktionen har du olika strategier till ditt f¨orfogande:
¨
¨
¨
– Andra
p˚a felrisken α – Andra
p˚a f¨ors¨oksfelvariansen σ2 – Andra
p˚a stickprovsstorleken n
• Pr¨ova vilken strategi som a¨r b¨ast f¨or att n¨arma sig den ”ideala” styrkefunktionen. Vilken strategi tror
du a¨r enklast att genomf¨ora i praktiken vid t.ex. laboratoriearbete?
10
Enkel linj¨ar regression
¨
Ovningens
syfte:
• Att studera modellen ”enkel linj¨ar regression” och illustrera hur de olika parametrarna i modellen
p˚averkar data;
• att illustrera begreppet residualer;
• att visa p˚a de olika m¨ojligheter som finns i Matlab att analysera en regressionsmodell.
Kurskompendium:
Olbjer kap 10
F¨orberedelseuppgifter:
Antag att givet a¨r talpar (xi , yi ), i=1,. . .,10 d¨ar man anser att sambandet mellan x och y a¨r linj¨art. Modellen
a¨r yi = α + βxi + εi d¨ar εi a¨r oberoende observationer fr˚an N(0, σ2 ).
1. Vad a¨r den grafiska tolkningen av α och β i modellen? Vilka f¨oruts¨attningar har du p˚a x-v¨ardena?
2. Regressionsmodeller beskrivs enklast med hj¨alp av matriser. Hur ser matrisformuleringen ut f¨or
modellen ovan?
3. Residualanalys a¨r ett viktigt instrument vid analys av regressionsmodeller. Hur definieras residualerna i ovanst˚aende modell? Skattningen av σ2 i modellen bygger p˚a residualerna. Hur ser skattningen ut?
10.1
Illustration av modell:
I ett enkelt simuleringsexperiment ska du unders¨oka hur v¨ardet p˚a σ p˚averkar modellen och de slutsatser
man kan dra fr˚an data. (F¨or att ge illustrativa bilder ges fullst¨andiga Matlab kommandon i denna del
av laborationen.) Skapa en vektor x med v¨arden 1, 2, . . ., 10 och en variabel y som erh˚alls genom det
teoretiska linj¨ara sambandet y=α+βx, d¨ar α och β a¨r k¨anda. V¨alj t ex y=10+2x. Till variabeln y adderas tv˚a
upps¨attningar av normalf¨ordelade m¨atfel N(0, σ2 ) med olika v¨arden p˚a σ, f¨orslagsvis σ=1 och σ=5.
15
11
KALIBRERINGSKURVA
>> x=[1:10]’
>> y1=10+2*x+normrnd(0,1,10,1);
>> y2=10+2*x+normrnd(0,5,10,1);
Vektorn y1 best˚ar allts˚a nu av 10 observationer fr˚an N(10+2x, 1) medan y2 best˚ar av 10 observationer fr˚an
N(10+2x, 25).
• Titta p˚a data i samma diagram och j¨amf¨or.
>>
>>
>>
>>
plot(x,10+2*x)
hold on
plot(x,y1,’x’)
plot(x,y2,’o’)
• F¨or att skatta regressionslinjen och titta p˚a residualerna utnyttjar vi den specialskrivna m-filen
reggui.
>> reggui(x,y1)
>> reggui(x,y2)
• Titta p˚a residualerna f¨or de b˚ada linjerna. Hur p˚averkas de av v¨ardet p˚a σ?
• I figurerna som alstras av reggui ges a¨ven skattningar och konfidensintervall f¨or modellens parametrar. J¨amf¨or de erh˚allna intervallen med de sanna v¨ardena p˚a α och β; t¨acker intervallen o¨ ver parametrarna?
10.2
Matlabs egen inbyggda regressionsrutin
I Matlab finns en inbyggd funktion f¨or regressionsanalys, regress, som kan anv¨andas vid multipel linj¨ar
regression (och d¨armed f¨orst˚as a¨ven vid enkel linj¨ar regression). Observera att reggui endast kan anv¨andas
vid enkel linj¨ar regression samt vid polynomregression som a¨r ett specialfall av multipel linj¨ar regression.
Pr¨ova hj¨alpkommandot help regress f¨or att ta reda p˚a hur in- och utargumenten ser ut.
• Anv¨and regress f¨or att skatta en av de tv˚a regressionslinjerna ovan. D˚a m˚aste vi f¨orst bilda matrisen X som a¨r en (10 × 2) matris med f¨orsta kolumnen enbart ettor och andra kolumnen best˚aende
av x-v¨ardena.
>> X=[ones(10,1) x]
>> [b bint r]=regress(y1,X,0.05)
Utargumentet bint ger konfidensintervall f¨or parametrarna α och β (med konfidensgrad 0.95 h¨ar ovan).
• Kontrollera att de erh˚allna skattningarna och intervallen st¨ammer med de du fick fr˚an reggui.
11
Kalibreringskurva
¨
Ovningens
syfte:
• Att studera modellen ”enkel linj¨ar regression” och illustrera hur de olika parametrarna i modellen
p˚averkar data;
• att illustrera begreppet residualer;
• att visa p˚a de olika m¨ojligheter som finns i Matlab att analysera en regressionsmodell.
16
11
KALIBRERINGSKURVA
Kurskompendium:
Olbjer kap 10
F¨orberedelseuppgifter:
Antag att givet a¨r talpar (xi , yi ), i=1,. . .,10 d¨ar man anser att sambandet mellan x och y a¨r linj¨art. Modellen
a¨r yi = α + βxi + εi d¨ar εi a¨r oberoende observationer fr˚an N(0, σ2 ).
1. F¨or ett givet v¨arde p˚a x, x0 , a¨r man ofta intresserad av det f¨orv¨antade v¨ardet f¨or y, μ0 . Ange formeln
f¨or ett 95% konfidensintervall f¨or μ0 i ovanst˚aende modell.
2. F¨or ett givet v¨arde p˚a x, x0 , a¨r man ofta intresserad av det predikterade v¨ardet f¨or y, y(x0 ). Ange
formeln f¨or ett 95% prediktionsintervall f¨or y(x0 ) i ovanst˚aende modell.
3. Vad a¨r skillnaden mellan konfidensintervallet och prediktionsintervallet i de f¨oreg˚aende uppgifterna?
Anv¨and g¨arna ett konkret exempel f¨or att klarg¨ora skillnaden.
4. Vid kalibrering har man det omv¨anda problemet: F¨or ett givet v¨arde p˚a y, y0 , vill man skaffa ett
intervall f¨or motsvarande x0 . Visa grafiskt hur detta kan g¨oras utg˚ande fr˚an ett prediktionsintervall.
Man vill g¨ora en kalibreringskurva f¨or en kalorimetrisk analys av fluorjoner i vatten och m¨ater d¨arf¨or
transmittansen tv˚a oberoende g˚anger f¨or ett antal k¨anda koncentrationer av fluorjoner. Resultat (finns i fil
kalibrer.mat):
Konc F− (μg/ml) x
0.608
1.216
1.824
2.432
3.040
% transmittans y
80.3
80.3
80.5
80.4
80.9
81.0
81.2
81.8
81.6
82.0
Konc F− (μg/ml) x
3.648
4.256
4.864
5.472
6.080
% transmittans y
82.9
82.5
83.0
83.1
83.9
84.0
84.0
84.0
85.0
84.8
Materialet blir nog mer hanterbart om man arbetar med
>> x=[Konc’; Konc’];
>> y=[Trans(1,:)’ ; Trans(2,:)’];
• Pr¨ova att anpassa en enkel linj¨ar regressionsmodell till data med hj¨alp av reggui (observera att
reggui beh¨over ej den inledande kolumnen av ettor)
• Verifiera att modellen a¨r rimlig genom att titta p˚a residualerna.
• Prediktion: Vad a¨r den f¨orv¨antade transmittansen d˚a fluorkoncentrationen a¨r 5.0? Vad a¨r motsvarande
95% prediktionsintervall?
• Kalibrering (invers prediktion): D˚a man i framtiden ska anv¨anda linjen som kalibreringskurva, vill
man till ett v¨arde y best¨amma ett intervall som med 95% sannolikhet t¨acker provets verkliga halt.
Skatta ett 95% intervall f¨or fluorkoncentrationen d˚a man f¨or ett prov med ok¨and koncentration
avl¨ast trans=82.8.
I kemilaborationen ”Best¨amning av koppar med atomabsorptionsspektrofotometri med grafitugn” b¨orjade
ni med att ta upp en kalibreringskurva som ni sedan anv¨ande f¨or att ett- och ett ”¨overs¨atta” era uppm¨atta
17
¨ REGRESSION OCH PROBLEMET MED KOLINJARA
¨
MULTIPEL LINJAR
X-VARIABLER
12
m¨atv¨arden till kalibrerade m¨atv¨arden. N¨ar man sedan anv¨ander sina kalibrerade m¨atv¨arden s˚a tar man allts˚a
inte h¨ansyn till os¨akerheten i kalibreringskurvan. F¨or att g¨ora det kan man g¨ora ett kalibreringsintervall i
st¨allet f¨or det konfidensintervall ni gjorde under punkt 8.
De metoder som tas upp i boken f¨or att g¨ora kalibreringsintervall baseras p˚a ett observerat y-v¨arde, men
era m¨atningar baseras p˚a 3 respektive 10 observerade y-v¨arden, s˚a t.ex. kalibreringsintervallet i formelsamlingen f˚ar modifieras n˚agot till


s
∗
2
s
1 1 (x0 − ¯x ) 
¯y0 − α∗
I x0 = x0∗ ± t α/2 (n − 2) · ∗ ·
+ +
d¨ar x0∗ =
β
k n
Sxx
β∗
d¨ar 1k under kvadratroten tidigare var 1. k a¨r h¨ar antalet observerade y-v¨arden. Anv¨and detta, eller matlabkommandot kalibk f¨or att g¨ora ett kalibreringsintervall baserat dels p˚a 10-serien och dels p˚a 3-serien.
J¨amf¨or resultatet med konfidensintervallen ni fick under punkt 8.
12
Multipel linj¨ar regression och problemet med kolinj¨ara x-variabler
¨
Ovningens
syfte:
• Att ge exempel p˚a en korrelationsmatris;
• att ge exempel p˚a en multipel linj¨ar regressionsmodell;
• att illustrera hur multipel linj¨ar regression kan p˚averkas av kolinj¨ara x-variabler;
• att illustrera vilka kriterier som anv¨ands d˚a man v¨aljer mellan olika modeller.
Kurskompendium:
Olbjer kap 11.1-11.7
F¨orberedelseuppgifter:
Antag att y, responsvariabeln, beror linj¨art av tv˚a oberoende variabler x1 och x2 . Vid 10 olika f¨ors¨ok har
man noterat (x1 , x2 , y). Modellen a¨r nu yi = β0 + β1 x1i + β2 x2i + εi , i = 1, . . . , 10 och d¨ar εi a¨r oberoende
observationer fr˚an N(0, σ2 ) som tidigare.
1. Ange matriserna i matrisformuleringen av modellen.
2. Vad menas med att x-variablerna a¨r kolinj¨ara?
12.1
Cementdata
Detta experiment beskrevs i Industrial and Engineering Chemistry redan 1932. I 13 f¨ors¨ok har man m¨att
v¨armeutvecklingen i stelnande cement som funktion av viktprocenten av n˚agra ing˚aende a¨mnen. I filen
cement.mat finns f¨oljande variabler:
varme
cem1
cem2
cem3
cem4
utvecklad v¨arme i kalorier per gram cement
viktprocent av 3CaO·Al2 O3
viktprocent av 3CaO·SiO2
viktprocent av 4CaO·Al2 O3 ·Fe2 O3 ,
viktprocent av 2CaO·SiO2
Data:
18
13
¨ REGRESSION
POLYNOMREGRESSION – ETT SPECIALFALL AV MULTIPEL LINJAR
cem1
cem2
cem3
cem4
varme
7
1
11
11
7
11
3
1
2
21
1
11
10
26
29
56
31
52
55
71
31
54
47
40
66
68
6
15
8
8
6
9
17
22
18
4
23
9
8
60
52
20
47
33
22
6
44
22
26
34
12
12
78.5
74.3
104.3
87.6
95.9
109.2
102.7
72.5
93.1
115.9
83.8
113.3
109.4
Man vill kunna avg¨ora vilken eller vilka av de fyra cementvariablerna som ska anv¨andas i en modell f¨or att
f¨oruts¨aga v¨armeutvecklingen. Man ans¨atter en multipel linj¨ar regressionsmodell d¨ar varme a¨r responsvariabel och cementvariablerna oberoende variabler. Problemet a¨r att vissa av de fyra cementvariablerna samvarierar kraftigt med varandra vilket kommer att p˚averka regressionsanalysen.
• Unders¨ok om det finns n˚agon samvariation mellan de fem olika variablerna genom att t.ex. ber¨akna
korrelationsmatrisen mellan variablerna (corrcoef). Med kommandot corrcoef([cem1 cem2
cem3 cem4 varme]) f˚ar du en matris d¨ar elementen p˚a t.ex. rad 2 a¨r korrelationskoefficienterna
mellan cem2 och variablerna cem1, cem2, cem3, cem4 respektive varme. Plotta i tv˚adimensionella diagram de variabler mot varandra som verkar samvariera.
• Ber¨akna f¨orst en full regressionsmodell med varme som responsvariabel och samtliga fyra cementvariabler som f¨orklarande variabler (anv¨and regress och kom ih˚ag att ”designmatrisen” X ska
inledas med en kolumn av ettor). Vilka koefficienter a¨r signifikant skilda fr˚an noll? Enligt detta
¨ svaret rimligt?
resultat, vilka variabler ska vara med i modellen? Ar
• Pr¨ova olika kombinationer av f¨orklarande variabler i regressionsmodellen (anv¨and den nyss ber¨aknade korrelationsmatrisen f¨or att ”gissa” vilka variabler som b¨or vara med). J¨amf¨or konfidensintervallen f¨or β -koefficienterna, skattningen av σ och residualerna f¨or de olika modellerna. Vilka
cementvariabler verkar vara de viktigaste f¨or att f¨oruts¨aga v¨armeutvecklingen?
• Funktionen stepwise kan vara till stor hj¨alp vid modellvalet, se help stepwise. Skriv t ex
>>
13
stepwise([cem1 cem2 cem3 cem4],varme)
Polynomregression – ett specialfall av multipel linj¨ar regression
¨
Ovningens
syfte:
• Att illustrera polynomregression.
Kurskompendium:
Olbjer kap 11.9
F¨orberedelseuppgifter:
1. Hur ser modellen ut vid polynomregression;
2. j¨amf¨or modellen med en multipel linj¨ar regressionsmodell!
19
¨
BRODBAK
14
13.1
Kontroll av cyanid
F¨oljande problem presenteras i Meier & Zund: Statistical methods in analytical chemistry. Vid en kemisk
industri ville man utveckla en metod f¨or att unders¨oka m¨angden av cyanid (CN− ) i avloppsvattnet. I
litteraturen fann man en fotometrisk metod som verkade vara l¨amplig och n¨asta steg var att finna en
l¨amplig kalibreringsfunktion. En cyanidl¨osning gjordes och en elektromekanisk automat anv¨andes f¨or att
erh˚alla en sp¨adningsserie med 10, 30, . . ., 250 μg CN− /100 ml. Till 10 ml av kalibreringsl¨osningen
tillsattes 90 ml av den f¨arggivande reagentl¨osningen varefter absorbansen m¨attes (1 cm kyvett). Data finns
i fil cyanid.mat.
• Beskriver en linj¨ar kalibreringskurva sambandet v¨al?
• Ans¨att ist¨allet en kvadratisk kalibreringskurva: y = β0 + β1 x + β2 x 2 + ε. Denna modell a¨r ett
specialfall av en multipel linj¨ar regressionsmodell; vilka a¨r de tv˚a oberoende variablerna?
• Anv¨and kommandot regress f¨or att pr¨ova denna modell p˚a data och j¨amf¨or med den linj¨ara
kalibreringskurvan. Vilken av dem passar b¨ast till data?
Bokens f¨orfattare ger f¨oljande kemiska f¨orklaring: ”Stray light in the photometer dominates at high absorbances, which can be responsible for lower slopes at high concentrations. The chemical workup can
also produce lower yields of the chromophore at higher concentrations.”
13.2
Anpassning av polynom i Matlab
Med hj¨alp av funktionen reggui kan man interaktivt, med hj¨alp av minsta kvadratmetoden, anpassa
polynom av olika gradtal till ett datamaterial.
• Pr¨ova reggui p˚a data fr˚an f¨oreg˚aende exempel f¨or att snabbt j¨amf¨ora de tv˚a olika modellerna ovan.
13.3
Simulering av sj¨okalkning
I en laboration, d¨ar sj¨okalkning skulle simuleras i en tankreaktor, ville man unders¨oka hur uppl¨osningshastigheten f¨or kalkstensmj¨ol kunde beskrivas som funktion av v¨atejonskoncentrationen. I filen kalk.mat
finns f¨or olika pH-v¨arden uppm¨att logaritmerad reaktionshastighet f¨or kornstorlekar 71-90 μm.
• Anv¨and reggui f¨or att unders¨oka sambandet mellan de tv˚a variablerna.
14
Br¨odbak
¨
Ovningens
syfte:
• Att illustrera 22 -f¨ors¨ok;
• att visa p˚a betydelsen av att studera samspel mellan faktorer;
• att varna f¨or f¨ors¨oksuppl¨agget ”variera en faktor i taget”.
Kurskompendium:
Olbjer kap 12.1-12.3
20
15
NITRERINGSPROCESS
F¨orberedelseuppgifter:
1. Hur ser modellen ut vid ett 22 -f¨ors¨ok?
2. Hur tolkas effekterna A, B och (AB) och hur g¨or man konfidensintervall f¨or dem?
3. Vad skiljer ett faktorf¨ors¨ok fr˚an ett ”variera en faktor i taget” f¨ors¨ok. Vilka slutsatser kan du dra fr˚an
respektive f¨ors¨ok?
Vid framst¨allning av vetemj¨ol a¨mnat f¨or br¨odbak ville man unders¨oka hur tillsats av tv˚a mindre ingredienser (A och B) p˚averkade vetemj¨olets bakegenskaper. Vid ett experiment bakades ett antal limpor av
”standardsort” p˚a 4 olika s¨att (lite/mycket tillsats av A, lite/mycket tillsats av B) varefter br¨odets volymitet
((dm)3 /kg) best¨amdes. Marknadsunders¨okningar har visat att medelkonsumenten f¨oredrar h¨oga, fluffiga
br¨od vilket motsvarar h¨og volymitet. Resultat (finns i fil bak.mat):
(1)
(a)
(b)
(ab)
A
B
lite
mycket
lite
mycket
lite
lite
mycket
mycket
volymitet
2.36
2.00
1.70
3.22
2.03
2.34
2.01
2.98
2.31
2.28
1.83
2.76
• Unders¨ok, genom att g¨ora l¨ampliga konfidensintervall, vilka effekter tillsats av A och B har p˚a
br¨odets volymitet.
• Personen K t¨anker l¨agga upp sitt f¨ors¨ok n˚agot annorlunda. Fr˚an ”grundtillst˚andet” (lite A och lite
B) studerar han var som h¨ander d˚a man dels tills¨atter mycket A och dels tills¨atter mycket B. Han
betraktar allts˚a enbart de tre o¨ versta raderna i ovanst˚aende schema. Vilka slutsatser kommer K att
dra om A:s och B:s inverkan p˚a br¨odets volymitet?
• Vilka argument kan du anv¨anda f¨or att o¨ vertyga K om det ol¨ampliga i hans f¨ors¨oksplan?
15
Nitreringsprocess
¨
Ovningens
syfte:
• Att illustrera 23 -f¨ors¨ok;
• att visa p˚a betydelsen av att studera samspel mellan faktorer.
Kurskompendium:
Olbjer kap 12.4
F¨orberedelseuppgifter:
1. Hur ser modellen ut vid ett 23 -f¨ors¨ok?
2. Hur skattas σ2 om du har n replikat i varje niv˚akombination?
3. Hur skattas σ2 om du endast har en observation i varje niv˚akombination?
P˚a ett laboratorium unders¨okte man hur utbytet av en nitreringsprocess p˚averkades av f¨oljande faktorer:
A – tid f¨or salpetersyratillsats,
B – omr¨orningstid,
21
16
RESPONSYTOR
C – bottensatsf¨orekomst.
Faktorn C togs med i laboratorief¨ors¨oket d¨arf¨or att det a¨r vanligt i fabrikstillverkning att man inte reng¨or
processk¨arlet f¨or varje tillverkningsomg˚ang utan l˚ater en bottensats kvarst˚a till n¨asta omg˚ang. Varje faktor
har tv˚a niv˚aer. Utbytet (i procent) vid de olika niv˚akombinationerna anges i nedanst˚aende tabell (Finns i
fil nitrer.mat).
(1)
(a)
(b)
(c)
(ab)
(ac)
(bc)
(abc)
A
2
7
2
2
7
7
2
7
B
0.5
0.5
4.0
0.5
4.0
0.5
4.0
4.0
C
utan
utan
utan
med
utan
med
med
med
utbyte 1
87.1
88.3
81.8
86.5
82.9
89.0
83.3
83.7
utbyte 2
87.3
88.6
82.1
86.6
83.1
89.3
83.6
83.8
• Analysera f¨ors¨oket genom att skatta olika effekter. Hur p˚averkar de tre olika faktorerna processen?
• Man trodde sig vara s˚a pass f¨ortrogen med den kemiska reaktionen ifr˚aga att man kunde anta att
samspelseffekterna var ringa. St¨ammer det?
• I en kommande f¨ors¨oksupps¨attning av nitreringsprocessen vill man spara tid och pengar genom att
f¨or varje niv˚akombination enbart m¨ata ett utbyte. Tycker du detta a¨r l¨ampligt?
16
Responsytor
¨
Ovningens
syfte:
• Att illustrera skattningar av responsytor;
• att illustrera optimering i kemiska system
Kurskompendium:
Olbjer kap 12.10
F¨orberedelseuppgifter:
1. Vad a¨r en responsyta?
2. Hur kan problemet att skatta en responsyta o¨ verf¨oras till ett multipel linj¨art regressionsproblem?
3. Vilka problemst¨allningar a¨r av intresse vid optimering av en responsyta?
F¨or olika typer av kemiska system d¨ar en responsvariabel y beror p˚a variabler x1 och x2 , y = f (x1 , x2 ), kan
responsytan f ofta l¨ampligt beskrivas med hj¨alp av en andra ordningens polynommodell:
Om yi a¨r responsvariabeln vid f¨ors¨ok i (t ex absorbans) och x1 och x2 a¨r faktorer eller experimentella
variabler (t ex koncentrationerna av tv˚a olika a¨mnen) ges modellen av
(1)
2
2
yi = β0 + β1 x1i + β2 x2i + β11 x1i
+ β22 x2i
+ β12 x1i x2i + εi
Modell a¨r allts˚a en multipel linj¨ar regressionsmodell med 5 oberoende variabler, x1 , x2 , x12 , x22 och x1 x2 och
d¨ar slumpvariablerna εi antas vara oberoende och normalf¨ordelade
N(0, σ2 ).
22
16
RESPONSYTOR
100
80
60
40
20
0
100
y
y
N˚agra responsytor med respektive polynom:
80
60
y
40
20
x2
0 0
100
80
60
40
20
0
100
80
20
40
x1
60
40
20
x2
0 0
20
40
x1
60
60
80
80
100
100
80
60
40
20
0
100
80
60
40
20
x2
0 0
20
40
x1
60
80
100
100
Responsyta 1: y = 96.7 − 1.68x1 − 2.18x2 + 0.013x12 + 0.018x22 + 0.01x1 x2
Responsyta 2: y = 3.35 + 1.68x1 + 2.18x2 − 0.013x12 − 0.018x22 − 0.01x1 x2
Responsyta 3: y = 93.3 − 0.37x1 − 1.37x2 − 0.005x12 + 0.005x22 + 0.02x1 x2
Fundera en stund p˚a hur parametrarna β0 , β1 , . . . , β12 ska tolkas. (J¨amf¨or med funktionerna i figurerna
ovan.) Hur p˚averkar storleken p˚a β11 och β22 responsytan? Hur ser ytan ut om b˚ade β11 och β22 a¨r
positiva? eller om b˚ada a¨r negativa?
Parametern β12 a¨r ofta av speciellt intresse eftersom den uttrycker samspelet mellan variablerna x1 och x2 . I
responsyta 3 a¨r parametern β12 signifikant skild fr˚an 0. Vad h¨ander med responsvaribeln y d˚a x1 a¨r fixerad
till 0 men x2 o¨ kar? Vad h¨ander om x1 i st¨allet fixeras till 100?
Ett experiment med vanadin
I en artikel i Journal of Chemical Education 69, 7 (1992), 560-563 beskrivs ett experiment d¨ar man
unders¨oker hur y – absorbans hos en l¨osning av vanadinsulfat p˚averkas av x1 – antal droppar av 1% H2 O2
samt av x2 – antal droppar av 20% H2 SO4 .
Add dropwise the amounts of 1% H2 O2 and 20% H2 SO4 , in that order, that are specified by
the experimental design to 40 drops of stock VOSO4 solution (about 0.1 g dissolved in 250
ml of distilled water). Stir the resulting mixture, and allow it to equilibrate for 5 min after
the addition of H2 SO4 .
Then measure the percent transmittance at 460 nm using a visible spectrophotometer (eg.
Bausch & Lomb Spectronic 20.) Calculate the absorbance. The procedure is a modification
of a vanadium spot test used in organic analysis.
23
16
RESPONSYTOR
F¨ors¨oket a¨r upplagt s˚a att b˚ade x1 och x2 unders¨oks p˚a ”l¨agsta niv˚a”, ”l˚ag niv˚a”, ”mellan niv˚a”, ”h¨og niv˚a”
och ”h¨ogsta niv˚a”. Motsvarande antal droppar a¨r d˚a 8, 10, 15, 20 och 22. (Denna f¨ors¨oksplan kallas
central composite design) Resultat (finns i fil vanadin.mat):
x1 −
x2 −
y−
antal dr H2 O2 antal dr H2 SO4 absorbans
15
10
20
15
15
15
15
8
22
20
10
15
8
10
10
15
15
15
15
15
15
20
20
22
0.397
0.420
0.359
0.334
0.336
0.346
0.323
0.367
0.319
0.330
0.293
0.327
Kemisk kommentar till experimentet
En sur l¨osning av VOSO4 inneh˚aller vanadin(IV) i joner VO2+ (vanadyl(IV)-joner). N¨ar n˚agra droppar
H2 O2 s¨atts till en s˚adan l¨osning, bildas i f¨orsta hand den kraftigt r¨oda l¨osliga f¨oreningen (VO2 )2 (SO4 )3 ,
som a¨r en peroxof¨orening av vanadin(V). H2 O2 har h¨ar fungerat b˚ade som oxidationsmedel och som
ligand:
2VO2+ + 3H2 O2 + 3HSO4− → (VO2 )2 (SO4 )3 + 4H2 O + H +
r¨od
Vanadin(V) kan emellertid koordinera totalt 4 st peroxo- eller oxodigander i blandade oxoperoxovandat(V)joner, som t ex V (O2 )O33− , V (O2 )2 O23− etc. Sedan v¨al allt vanadin oxiderats till femv¨art, kan s˚aledes den
starkt f¨argade (VO2 )2 (SO4 )3 , som troligen ger den maximala absorbansen vid 460 nm, till st¨orre eller mindre del o¨ verg˚a till diverse, ljusare gula, v¨ateoxoperoxovandat(V)-joner – allt beroende p˚a m¨angden tillsatt
H2 O2 och H2 SO4 . Som typexempel kan f¨oljande j¨amvikt tj¨ana:
(VO2 )2 (SO4 )3 + 2H2 O2 + 4H2 O 2H2 V (O2 )2 O2− + 3HSO4− + 5H +
Tillsats av H2 O2 f¨orskjuter j¨amvikten a˚t h¨oger(f¨argen avtar); tillsats av H2 SO4 f¨orskjuter den a˚t v¨anster
(f¨argen djupnar).
Statistisk analys
Data fr˚an experimentet finns i fil vanadin.mat.
• G¨or en tv˚adimensionell plot o¨ ver x1 och x2 f¨or att se vilken f¨ors¨oksuppl¨aggning man anv¨ant sig
av. Kan du komma p˚a n˚agon f¨ordel med denna uppl¨aggning? (T¨ank p˚a avsnittet med kolinj¨ara
x-variabler)
• Skatta en responsyta med hj¨alp av respons (se help respons), skatta parametrarna i modell (1)
och studera dess konfidensintervall. (Eftersom modell (1) a¨r en multipel linj¨ar regressionsmodell
bygger responsyteskattningen p˚a kommandot regress.
• B¨or alla variablerna vara med? Hur ser den slutliga modellen ut?
24
16
RESPONSYTOR
• Rita upp den slutligt skattade responsytan. St¨amde den skattade responsytan med den kemiska
tolkningen?
• Finns det andra v¨arden p˚a x1 och x2 du skulle vilja pr¨ova i experimentet? Vilka i s˚a fall och varf¨or?
25