CHEMOMETRIKA a STATISTIKA Prozatímní učební text – vybrané příklady (srpen 2012) Miloslav Suchánek
Úkol č. 1 Maticové operace s využitím EXCELu V EXCELu jsou dvě důležité maticové operace, které nám pomohou při řešení dalších úloh. V maticových operacích pracujeme s tzv. maticovými vzorci, které uplatníme i při řešení úloh lineárních a kvadratických regresí. Matice je definována číselnými hodnotami a počtem řádků a sloupců, např. matice: A=[
1 2 2
1 3] 4
má tří řádky a dva sloupce, značíme ji A3x4 nebo A(3,4). Součin matic, C, je definován vztahem, např. pro matice A, B: C = A.B Obecně neplatí A.B = B.A! Součin dvou matic je definován jen pro matice, které mají společné „vnitřní“ indexy, např. A(nxp).B(pxm) = C(nxm) V tomto případě součin B.A není definován. Přesvědčíme se o tom v EXCELu. Součin matic naleznete ve funkcích, f(x), které si vyhledáte kliknutím na fx pod matematickými funkcemi a označením SOUCIN.MATIC. Pokud nemáte na liště tuto funkci, naleznete ji v Nástroje-Vlastní-Vložit. Postup při výpočtu součinu dvou matic: myší vyhledáte buňku s dostatečným prostorem pro výslednou matici C (nxm buněk); vyvoláte funkci SOUCIN.MATIC; označíte nejprve matici A, potom B, přesně podle návodu v okénku, které se Vám objeví (označení polí dělejte kliknutím myší na malou červenou šipku, potom myší označte pole, znovu klikněte na malou červenou šipku) ; OK; ve vyhledané buňce se objeví číslo; označíte myší prostor buněk nxm; kliknete F2, přičemž se objeví maticový vzorec v uvedené buňce; stisknete CTRL+levý SHIFT a kliknete ENTER; v dříve označených buňkách (nxm) se objeví další prvky matice. Pozorně zkontrolujte vyznačený prostor, zda odpovídá rozměru matice C. Součet dvou matic, např. A + B je definován vždy, Tedy např. 2 3
3 1
(A + B) =
5 6
A=
B=
3 3
5 3
8 4
Součet matic si sami naprogramujete v excelovském sešitu. Inverze matice je definována pouze pro čtvercové matice, tedy matice o stejném počtu řádků a -1 sloupců. V maticovém zápisu píšeme inverzi matice A jako A . -1 Platí: A.A = I (jednotková matice). V EXCELU naleznete inverzi matice ve funkci fx pod názvem INVERZE. Postup při výpočtu inverze matice:
myší vyhledáte buňku, kam umístíte prvek inverzní matice (1,1), okolo musíte mít dostatečný prostor pro výslednou inverzní matici (má stejný rozměr jako původní matice); označíte myší pole buněk o rozměru nxn; vyvoláte funkci INVERZE; podle pokynů v okénku označíte matici, kterou chcete invertovat; OK; ve vyhledané buňce se objeví číslo;kliknete F2; potom stisknete současně CTRL+levý SHIFT a kliknete ENTER; ve vyznačeném poli se objeví invertovaná matice. Transpozici matice, v maticovém zápisu A nebo A´ můžete vytvořit z matice pomocí následujících operací: myší označíte matici, kterou chcete transponovat; CTRL C; přejdete na jinou buňku a vyvoláte ÚpravyT Vložit jinak-Hodnoty-Transponovat-OK a dostanete transponovanou matici A . Symetrická matice je matice, která má stejnolehlé prvky mimo hlavní diagonálu shodné, např. T
2
7
7
3
A= Symetrická matice se transponováním nemění, tedy A = A. T
Diagonální matice je matice, která má, kromě hlavní diagonály, všechny ostatní prvky nulové. Zvláštním případem diagonální matice je matice jednotková, která má na hlavní diagonále jedničky. Determinant matice je číslo. V EXCELU naleznete tuto funkci pod názvem DETERMINANT, opět v matematických funkcích vyvoláním fx. Platí, že determinant transponované matice je stejný, jako T determinant matice původní, tedy A = A. Pro čtvercové matice stejného řádu (stejný počet řádků) je determinant součinu dvou matic stejný, jako součin determinantů těchto matic, tedy A.B= A.B. Postup při výpočtu determinantu matice: myší vyhledáte buňku, kam umístíte hodnotu determinantu; vyvoláte funkci DETERMINANT; označíte matici podle návodu v okénku; OK; v buňce se objeví číslo, což je výsledná hodnota determinantu. Stopa čtvercové matice se dána součtem jejích diagonálních prvků. Stopa transponované matice se rovná stopě původní matice. Čtvercová matice C se nazývá ortogonální, jestliže pro ni platí, že matice k ní transponovaná se rovná inverzní matici, takže platí rovnice: T
C =C
-1
a
T
C .C = 1.
V přiloženém excelovském sešitu máte definovány úlohy, které vyřešíte.
Úkol č. 2
Lineární regrese Pro řešení úloh z lineární regrese se v EXCELu se naučíme využívat jednak funkce LINREG, jednak pomocí maticových operací si sami sestavíme vlastní program. Lineární regrese je matematicko-statistická metoda, při níž prokládáme experimentální data regresním modelem, kterým může být přímka, rovina nebo nadrovina ve vícerozměrném prostoru. Regresní model vybíráme z nekonečného množství možností takovým způsobem, abychom splnili podmínku minimálního součtu čtverců odchylek experimentálních a regresních hodnot. Teoretický lineární model můžeme formulovat rovnicí:
Y = f (x j ; β p ) ve kterém xj jsou nezávisle proměnné veličiny, βp jsou parametry modelu. Tak např. pro lineární model s jednou nezávisle proměnnou platí rovnice:
Y = β 0 + β1 x , pro kvadratický model s jednou proměnnou rovnice:
Y = β 0 + β1 x + β 2 x2 . Odpovídající regresní modely jsou vlastně odhady teoretických modelů, takže obecně píšeme:
Yreg = f ( x j ; b p ) + e V této rovnici jsou hodnoty bp odhady βp, e náhodná chyba modelu, která má nulovou střední hodnotu a normální rozdělení N(0,σ2). Regresní rovnice musí vyhovovat podmínce:
U = ∑ ( yexp,i − Yreg ,i ) 2 = min , n
kde n je počet pozorování hodnot závisle proměnné. V maticovém vyjádření potom lineární model píšeme ve tvaru y = X.β kde y je sloupcový vektor (nx1), X je matice hodnot nezávisle proměnné rozměru (n x (j+1)), β je vektor neznámých parametrů rozměru ((j+1)x1). Lineární regresní model píšeme ve tvaru yreg = X.b + e ve kterém b je vektor odhadů β, e je vektor n hodnot náhodné chyby rozměru (nx1). Pro složky vektoru e platí E(ei) = 0, takže E(e) = 0n (nulová matice nx1). D(ei) = σ2, takže matice C(e) = σ2.In. Definujme si ještě kovarianční matici C = (XT.X)-1. Pro odhad vektoru β je vektor b odhadnut metodou nejmenších čtverců (viz předchozí výklad). V maticové podobě: přičemž kovarianční matice C(b) má tvar
b = (XT.X)-1XT.y, C(b) = syx2(XTX)-1
syx je nejlepším odhadem σ a počítá se podle vztahů: yreg = X.b
e = yexp - yreg
Q = eTe syx2 = Q/(n-p).
V dalším se dohodneme že vektory píšeme vždy ve sloupcové formě, takže transponovaný vektor je řádkový, dále se dohodneme, že matice a vektory jsou psány tučným písmem, skaláry normálním. Nyní k technice výpočtů. V EXCELU naleznete lineární regresi ve funkci fx pod názvem LINREGRESE. Data musíte před tím seřadit do tabulky po dvojicích jako dva sloupcové vektory (nx1), tedy y1 x1 x2 y2 ............. yn xn Najdete si buňku s místem kolem o rozměru 3x2. Po vyvolání LINREGRESE se Vám objeví tabulka, kterou vyplníte. Do pole B napíšete PRAVDA, do pole Stat rovněž PRAVDA. To Vám zaručí výpočet všech statistik. Po odeslaní OK se ve vybrané buňce objeví číslo. Kolem buňky myší vyznačíte maticový prostor o 3 řádcích a 2 sloupcích. Stisknete F2 a potom stisknete současně CTRL+levý SHIFT a kliknete ENTER; ve vyznačeném poli se objeví matice čísel. Čísla jsou hodnoty podle následujícího schématu: b0 sb0 syx
b1 sb1 R2
kde sb1 a sb0 jsou výběrové směrodatné odchylky parametrů, R2 je koeficient determinace, který nebudeme používat. syx je výběrová směrodatná odchylka závisle proměnné veličiny (viz předchozí text). Pomocí LINREGRESE můžete počítat i kvadratický regresní model. Je třeba ale přeorganizovat tabulku experimentálních hodnot takto: x1
x12
y1
atd. Při vyvolání LINREGRESE označíte pole x jako celou matici x, včetně kvadratických hodnot. Maticový prostor pro výsledky bude ale o rozměru 3x3 a hodnoty budou seřazeny takto:
b2 sb2 R2
b1 sb1 syx
b0 sb0 -
Pro maticové výpočty doplníte vektor hodnot x na matici X tímto způsobem: 1 1
x1 x2
atd. a budete postupovat podle maticových vzorců způsobem, který jsme se již naučili.
Úkol č. 3 Nelineární regrese V této úloze se naučíme řešit úlohy nelineární regrese pomocí ŘEŠITELE, což je velice účinný nástroj EXCELu pro řešení různých matematicko statistických problémů. Nejprve něco o nelineární regresi. Teoretický nelineární model můžeme, stejně jako v předchozím lineárním případu, formulovat rovnicí:
Y = f (x j ; β p ) ve kterém xj jsou nezávisle proměnné veličiny, βp jsou parametry modelu, p je index parametru. Tak např. pro nelineární model s jednou nezávisle proměnnou platí rovnice:
Y = β 0 exp( β 1 x ) , ve kterém p=2. Odpovídající regresní model je odhadem teoretického modelu, takže obecně píšeme:
Yreg = f ( x j ; b p ) + e V této rovnici jsou hodnoty bp odhady βp, e náhodná chyba modelu, která má nulovou střední hodnotu a normální rozdělení N(0,σ2). Tedy totožné s lineárním případem. Regresní rovnice musí vyhovovat podmínce:
U = ∑ ( yexp,i − Yreg ,i ) 2 = min , n
kde n je počet pozorování hodnot závisle proměnné. Řešení rovnice pro minimum čtverců odchylek je značně složitější a jde za rámec tohoto předmětu. Navíc, vše za nás udělali autoři ŘEŠITELE, takže my se naučíme využívat výsledky jejich práce. Nicméně trocha teorie bude nutná pro řešení dalších úloh. Přepišme si rovnici pro součet čtverců do jiného tvaru s využitím faktu, že tento součet, stejně jako v případě lineárního modelu, je funkcí parametrů modelu:
U ( b ) = ∑ ( yexp,i − f ( xi , b p )) 2 n
Derivováním U postupně podle všech bp dostáváme (stejně jako v případě lineárního modelu) soustavu tzv. normálních rovnic o p-proměnných. Potíž je v tom, že když se regresní funkce nelineární, jsou nelineární i normální rovnice. Tvar funkce U kolem minim b je p-rozměrný eliptický paraboloid, z čehož (bez důkazu) plyne, že řešení, t.j. nalezení minima funkce U vzhledem k parametrům bp je velice citlivé na počáteční odhady b, což v případě lineární regrese nenastává. Tam počáteční odhady parametrů b ( v lineárním případě směrnice a úsek) nepotřebujeme. Regresní parametry nelineárního vztahu lze jednoznačně odhadnout pouze v případech, kdy jednotlivé parciální derivace
∂f ( x, b ) / ∂b j (j ∈ 1,..p) jsou lineárně nezávislé. V EXCELu naleznete ŘEŠITELE (v angl. versi SOLVER) pod Nástroje-Řešitel. Kliknutím na Řešitel se Vám objeví panel, ve kterém musíte vyplnit některé údaje. Nejprve si ale vysvětlíme některé pojmy. Nastavit buňku: v této buňce bude výraz nebo rovnice, ve které se zobrazuje řešení problému. Pro náš případ, nelineární regrese, to bude suma čtverců odchylek experimentálních hodnot závisle proměnné a regresní proměnné. Musíte tedy znát tvar regresního modelu. Pozor, v našem případě to bude jenom jedna buňka.
Rovno: zde označíte cílové řešení, tedy v našem případě hledáme minimum sumy čtverců, označíte tedy Min Měněné buňky: zde značíte buňky, ve kterých bude řešení, tedy vektor parametrů. Hodnoty v těchto buňkách se během iterací mění, řešení však vyžaduje, abyste před spuštěním zadali počáteční hodnoty parametrů. Jak již bylo řečeno, výsledné řešení je silně závislé na počátečních hodnotách (nulové přiblížení) parametrů. Zkusíte si to při řešení úloh. Pokud kliknete Možnosti, objeví se Vám panel s dalšími variantami řešení, my prozatím ponecháme nastavené hodnoty (default). Po vyplnění hlavního panelu kliknete na Řešit. V listě máte v příslušných měněných buňkách vektor řešení, tj. hodnoty parametrů regresního modelu a sumu čtverců odchylek experimentálních a regresních hodnot závisle proměnné. Naučíme se používat ŘEŠITELE při řešení úlohy lineární regrese, tj. mat_2. Využijeme toho, že máte data přepsána do tvaru, který byl nutný pro lineární regresi. Jako nulového přiblížení použijeme mírně změněných výsledků lineární regrese, řekněme o 10 % vyšších hodnot. Tabulku vstupních dat musíte doplnit o vektor yi,reg s použitím vstupních hodnot. Příklad excelovského uspořádání: Regresní rovnice (viz minulý příklad): Yreg = -0,00161 + 0,09984.c Vstupní hodnoty parametrů umístěte do buněk v listu, kam jste překopírovali tabulku vstupních experimentálních hodnot (y,c), např. do buněk G3:G4 (G3: vstupní parametr pro b0, G4: vstupní parametr pro b1). Tabulku experimentálních hodnot budete mít např. ve sloupcích A a B, uspořádány takto: A c A2: A3:
0,1 0,1
B yexp B2: B3:
0,009 0,007
atd. celkem 22 řádků. Nyní do sloupce C umístíte hodnoty Yreg pro každou hodnotu c tímto způsobem: do buňky C2 napíšete =$G$3+$G$4*A2 a překopírujete do celého sloupce C postupem: CTRL+C (buňka C2); myší označíte buňky C3:C23; CTRL+V. Nyní máte ve sloupci C regresní hodnoty závisle proměnné pro všechny hodnoty nezávisle proměnné (koncentrace c). Pozor! Nezapomeňte před touto operací vložit do buněk G3 a G4 nulová přiblížení parametrů b0 a b1 !!!! Do sloupce D musíte nyní umístit rozdíly yexp - Yreg , např. takto: do buňky D2 napíšete: =B2-C2 a překopírujete stejným způsobem jako před tím do celého sloupce D. Následuje výpočet sumy čtverců tímto postupem: Aktivní buňka bude např. G6 (myš přesunete na tuto buňku a kliknete levým tlačítkem); na liště vyvoláte funkce (funkční tlačítko f(x)) a najdete funkci SUMA.ČTVERCŮ. Kliknete na tuto funkci a na panelu označíte buňky D2:D23. Po odeslání OK se Vám v buňce zobrazí součet čtverců odchylek experimentálních a regresních hodnot (samozřejmě se vstupními hodnotami parametrů). Nyní máte již připravenu tabulku pro použití ŘEŠITELE. Před použitím si ještě nějak označte pole vektorů parametrů a buňku cílové funkce (suma čtverců), abyste se v tom potom vyznali.
Výchozí tabulka vypadá takto: sl.\řád. 1 2 3 atd.
A c 0,01 0,01
B yexp 0,009 0,007
C Yreg =$G$3+$G$4*A2 =$G$3+$G$4*A3
D yexp - Yreg =B2-C2 =B3-C3
V buňce G6 je SUMA.ČTVERCU(D2:D24) a v buňkách G3 a G4 jsou vstupní hodnoty parametrů lineární regrese b0 a b1 (volte nejprve výsledky lineární regrese s 10 %ní změnou). Zůstanete v aktivním listu a kliknete Nástroje-Řešitel v Nastavit buňku označíte G6, Rovno=Min, Měněné buňky: G3:G4 a kliknete Řešit.Téměř okamžitě je řešení skončeno a objeví se panel Výsledky řešení, kde označíte Uchovat řešení a kliknete OK. To je vše, ani to snad nebylo tak obtížné. Zkuste si různé “nástřely“ vstupních hodnot b0 a b1, abyste se přesvědčili o citlivosti výsledku na vstupní parametry.
Úkol č. 4 Vícerozměrná pozorování Základním podkladem pro vícerozměrnou analýzu je datová matice typu (n x p). Řádky odpovídají jednotlivým studovaným objektům (n), sloupce jednotlivým zjišťovaným znakům p (příznakům, pozorovaným proměnným). Tohoto značení budeme v dalším důsledně používat. Datovou experimentální matici budeme označovat X , popř. Y. Prvek datové matice X, xij, je hodnota j-tého příznaku (j=1,2,...p) zjištěná u i-tého prvku (i=1,2,...n). Vektor datové matice X, xi, nazveme obrazem. Tak např. pro datovou matici 3 5 6 7 9
7 6 8 10 9
která je charakterizována dvěma příznaky, je prvním obrazem vektor x1 (3 7), pátým obrazem vektor x5 (9 9). Studované objekty bývají předměty (vzorky), události, instituce (laboratoř), apod. Používání statistických metod vyžaduje, aby byl studován přiměřeně rozsáhlý soubor objektů. Vedle pojmu objekt budeme používat také pojmy jednotka, individuum nebo prvek. Typickým cílem statistické analýzy je poznání vlastností objektů, popř. závislostí mezi těmito vlastnostmi. O úspěchu zjišťování přitom mimo jiné rozhoduje to, jak se podaří vyjádřit měřitelnými znaky jednotlivé vlastnosti, o které se zajímáme. Různorodý charakter zkoumaných proměnných je běžným jevem ve všech aplikacích vícerozměrné statistické analýzy. Jako příklad lze uvést záznam o kontrole jakosti výrobku, protokoly o kontrole stavu životního prostředí, atd. Klasifikace proměnných a klasifikace objektů je velice složitý problém. V dalším se budeme zabývat základy klasifikací proměnných a objektů, které se uplatňují při interpretaci dat hlavně z oblasti životního prostředí. Příznaky, měřené veličiny, můžeme rozdělit podle toho, zda popisují kvalitativní nebo kvantitativní charakteristiky objektů. Nyní se budeme věnovat vícenásobnému porovnávání základních statistik, tj, středních hodnot a rozptylů. Nejprve si zopakujeme porovnávání pro jeden příznak (p=1), které jsme probírali v Chemometrice I. Toto porovnání znáte pod názvem Analýza rozptylu (ANOVA). Analýza rozptylu pro jeden faktor Uvažujme experiment, v němž je vyšetřován vliv jednoho faktoru, např. A, který bude sledován na k úrovních (k>2). Při každé úrovni provedeme stejný počet opakování měření, r, přičemž pro celkový počet pokusů n platí: n = r.k Výsledky pokusů tvoří tzv. experimentální matici, jejíž obecný člen označíme yiν, kde ν je počet opakování měření na i-té úrovni. Pro i-tou úroveň můžeme model pro analýzu rozptylu vyjádřit vztahem: yiν = µ + αi + eiν ve kterém µ je střední hodnota závisle proměnné pro všechny úrovně, neboť experimentální matici si můžeme představit jako náhodný výběr ze základního souboru. αi je vliv faktoru A na i-té úrovni. Definujme si nyní pomocné mezisoučty v experimentální matici tak, jak je v analýze rozptylu zvykem:
Y.. = ∑∑ yiν k
r
Yi . = ∑ yiν r
Odhady parametrů jsou potom vyjádřeny rovnicemi:
µ= 1 n Y..
α= 1 r Yi . − µ eiν = yiν − 1 r Yi . Nyní budeme testovat hypotézu o tom, že vlivy faktoru A na všech úrovních jsou stejné, tedy hypotézu: H0: α1 = α2 = ...= αk proti alternativní hypotéze H: Σαi > 0 Testové kriterium F je potom vyjádřeno vztahem:
SA F=
Sr
( k − 1)
(n − k)
ve kterém
S A = 1r
∑k Yi
S r = ∑∑ yiν 2 − 1 r k r
2 .
− 1 Y..
∑k Yi
2
n 2
.
Hodnotu F srovnáváme s F1-α,(k-1;n-k) Odhad rozptylu měřené veličiny vypočteme z residuální proměnlivosti: s2 = Sr/(n-k)
Vaším úkolem bude sestavit jednoduchý list pro jednofaktoriální analýzu rozptylu tak, aby bylo možno proces vyhodnocování zautomatizovat. List sestavíte v EXCELu pro počet úrovní k = 10 a počet opakování max. 10 pro každou úroveň. List konstruujte tak, aby se dal použít i pro menší počet úrovní i opakování a měl hezkou grafickou úroveň. Váš list zkontrolujte s analýzou rozptylu, která je v EXCELu v Nástroje-Analýza dat-Anova:jeden faktor.
Porovnání rozptylů na více úrovních
Probereme si ještě další způsob porovnání, a to porovnání rozptylů pro více úrovní jednoho faktoru. Potvrzujeme či vyvracíme přitom hypotézu: H0: σ21 = σ22 = σ2k Test, kterým ověřujeme nulovou hypotézu, se nazývá Batlettův test a spočívá ve výpočtu testového kriteria B a porovnání s χ2 kvantilem.
Bartlettův test B=
1
C
k
[(n − k ) ln s 2 − ∑ (nh − 1) ln s h2 ] h =1
s2 =
C = 1+
1 k ∑ (nh − 1)s h2 n − k h =1
k 1 1 1 ((∑ )− ) 3(k − 1) h =1 nh − 1 n − k
V těchto vztazích je nh počet pozorování na h-té úrovni (h = 1,2,...k) Testovací kriterium: B< χ2(1-α)(k-1) V přiloženém excelovském sešitu máte definovány úlohy, které vyřešíte.
