Milan Bernauer, Bohumil Bernauer, Petr Št astný. dat. Ústav Anorganické Technologie. Obsah

ˇ ’astn´ Milan Bernauer, Bohumil Bernauer, Petr St y

Statistick´ e zpracov´ an´ı namˇ eˇ ren´ ych dat ´ ´ Technologie Ustav Anorganicke

Obsah 1 Z´ akladn´ı pojmy 1.1 Výbˇerové charakteristiky . . . . . . . . . . . . . . 1.2 Aplikace výbˇerových charakteristik - I . . . . . . . . 1.3 Aplikace výbˇerových charakteristik - II . . . . . . .

1 2 3 4

2 Z´ akon ˇ s´ıˇ ren´ı chyb

5

3 Regresn´ı anal´ yza 3.1 Lineárn´ı regrese . . . . . . . . . . . . . . 3.2 Váˇzená lineárn´ı regrese . . . . . . . . . . 3.3 Urˇcen´ı poˇctu parametr˚ u modelové rovnice 2 3.4 χ test konsistence sd a s(yi) . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

7 8 11 12 13

´ Ustav anorganick´ e technologie: Statistické zpracován´ı namˇeˇren´ ych dat

1

Z´ akladn´ı pojmy

Chyby vyskytuj´ıc´ı se pˇri experimentáln´ı práci m˚ uˇzeme rozdˇelit do dvou kategori´ı: N´ ahodné chyby Systematické chyby.

Opakovan´ ym mˇeˇren´ım za identick´ ych podm´ınek na stejném zaˇr´ızen´ı neobdrˇz´ıme vˇzdy identické v´ ysledky, coˇz je d˚ usledkem fluktuac´ı veliˇcin bˇehem experimentu (zmˇeny teploty, tlaku...), resp. vlivem náhodn´ ych chyb. Náhodné chyby m˚ uˇzeme popsat vcelku snadno a rigoróznˇe vztahy uveden´ ymi v následuj´ıc´ı sekci 1.1. Systematické chyby jsou m´ırou odlehlosti zjiˇstˇené veliˇciny od nˇejaké známé ˇci standardn´ı hodnoty. Systematické chyby m˚ uˇzeme roztˇr´ıdit do následuj´ıc´ıch skupin podle jejich p˚ uvodu Instrument´ aln´ı chyby (nejistoty v kalibraci mˇeˇr´ıc´ıch pˇr´ıstroj˚ u, neˇcistoty ve vzorc´ıch) Chyby experimentáln´ı metody (nedokonalost pouˇzit´ ych model˚ u popisuj´ıc´ıch danou metodu) Chyby experiment´ atora

Zjistit vliv systematick´ ych chyb na v´ ysledky mˇeˇren´ı je obt´ıˇzné pokud neznáme v´ yslednou (správnou) hodnotu, at’ jiˇz standardu, nebo napˇr´ıklad hodnotu z literatury. Pokud je tato hodnota nedostupná, je tu moˇznost provést mˇeˇren´ı identické veliˇciny na jiné aparatuˇre, nejlépe zaloˇzené na jiném principu neˇz ta p˚ uvodn´ı. Napˇr´ıklad mˇeˇrit teplotu rtut’ov´ ym a alternativnˇe odporov´ ym teplomˇerem. Dále se setkáváme s pojmy vyjadˇruj´ıc´ımi “kvalitu” experimentálnˇe z´ıskan´ ych dat; pˇresnost a správnost. Uvedené pojmy by mˇely slouˇzit k objektivn´ımu popisu dat, jak z pohledu jejich vzájemné vnitˇrn´ı konzistence a reprodukovatelnosti (pˇresnost) tak jejich vztahu k nˇejaké známé, správné hodnotˇe (správnost). Pˇresnost (precision) a správnost ˇ (accuracy) jsou definovány na základˇe normy CSN ISO 57251 a jejich v´ yznam je patrn´ yz obrázku 1. Pˇresnost mˇeˇren´ı je ovlivnˇena náhodn´ ymi chybami a m˚ uˇzeme ji vyhodnotit statistick´ ymi postupy, které budou uvedeny v následuj´ıc´ım odd´ıle 1.1. Správnost dat m˚ uˇzeme vyhodnotit pouze známe-li správnou (referenˇcn´ı) hodnotu. Hodnota experimentálnˇe stanovené veliˇciny, uvedená bez odhadu jej´ı nejistoty, je ochuzena o podstatnou ˇca´st své vˇedecké informace. V´ ysledek uveden´ y ve formˇe 25.0 ± 0.2 je diametrálnˇe odliˇsn´ y od v´ ysledku 25.0 ± 15.5, aˇc jsou absolutn´ı hodnoty identické 2 . Nejistota (uncertainty) se vˇetˇsinou udává jako odhad standardn´ı smˇerodatné odchylky dané veliˇciny 3 . ˇ CSN ISO 5725:Pˇresnost (spr´ avnost a shodnost) metod a v´ ysledk˚ u mˇeˇren´ı Malanowski, E.R. A Computer Program for Calculating Standard Deviations from Standard Deviations J. Chem. Educ. 1995, 72 , 1079 − 1082. 3 http://physics.nist.gov/cuu/Uncertainty/index.html 1 2

1 ˇ ’astn´ Last change 3.12.2009 Milan Bernauer & Bohumil Bernauer & Petr St y


Obrázek 1: V´ yznam pojm˚ u pˇresnost a správnost. Oba obrázky obsahuj´ı stejná data, pouze je zmˇenˇen zp˚ usob jejich vynesen´ı do grafu. ∆ = xi −µ, kde xi je zmˇeˇrená (zjiˇstˇená) hodnota a µ je správná (referenˇcn´ı) hodnota. - pˇresná a správná data; - nepˇresná a správná data; - pˇresná a nesprávná data (dle teorie poznán´ı plynouc´ı z filozofie externismu Járy Cimrmana: ”Nev´ıme nic, ale zato to v´ıme pˇresnˇe ”); - nepˇresná a nesprávná data.

@

E

1.1

p

u

V´ ybˇ erov´ e charakteristiky

Mˇejme náhodnou veliˇcinu Y , kterou z´ıskáváme opakovan´ ym mˇeˇren´ım, napˇr´ıklad nˇekteré fyzikáln´ı veliˇciny. Obvykle se z tˇechto opakovan´ ych mˇeˇren´ı vyhodnocuj´ı tyto veliˇciny: pr˚ umˇerná hodnota a odhad rozptylu, resp. smˇerodatná odchylka, které poskytnou odhad o stˇredn´ı hodnotˇe a pˇresnosti, resp. nejistotˇe veliˇciny Y . Slovo odhad je zde pouˇzito z d˚ uvodu koneˇcného poˇctu mˇeˇren´ı veliˇciny Y . Aby bylo moˇzno stanovit stˇredn´ı hodnotu µ a rozptyl σ 2 veliˇciny Y (maj´ıc´ı Gaussovo rozdˇelen´ı, téˇz naz´ yváno normáln´ı rozdˇelen´ı) muselo by se provést nekoneˇcnˇe mnoho experiment˚ u. Z evidentn´ıch d˚ uvod˚ u tento poˇcet mˇeˇren´ı nem˚ uˇzeme uskuteˇcnit a proto se vyhodnocuj´ı odhady z vybraného (koneˇcného) poˇctu mˇeˇren´ı na základˇe tzv. v´ ybˇerov´ ych funkc´ı4 . Tyto funkce jsou: V´ ybˇerov´ y pr˚ umˇer

n

Y =

1X Yi n i=1

V´ ybˇerov´ y rozptyl

(1)

n

s2 (Y ) =

1 X (Yi − Y )2 n − 1 i=1

(2)

p s2 (Y )

(3)

V´ ybˇerová smˇerodatná odchylka

s(Y ) = 4

V´ ybˇerové funkce jsou definované na v´ ybˇerovém prostoru a jejich rozdˇelen´ı pravdˇepodobnosti je urˇceno pravdˇepodobnostn´ım rozdˇelen´ım z´ akladn´ıho souboru.



Smˇerodatná odchylka v´ ybˇerového pr˚ umˇeru

s(Y ) s(Y ) = √ n

(4)

K v´ ybˇerové smˇerodatné odchylce ze vztahu 3 a smˇerodatné odchylce v´ ybˇerového pr˚ umˇeru ze vztahu 4 existuj´ı relativn´ı veliˇciny, které jsou definovány jako pomˇer v´ ybˇerové smˇerodatné odchylky a absolutn´ı hodnoty veliˇciny Y sr (Y ) =

s(Y ) , Y

(5)

respektive jako pomˇer smˇerodatné odchylky v´ ybˇerového pr˚ umˇeru a pr˚ umˇerné hodnoty Y sr (Y ) =

s(Y ) . Y

(6)

V pˇr´ıpadˇe, kdy známe rozptyl jednotlivého pozorován´ı s2 (Yi ), napˇr´ıklad kdyˇz tato veliˇcina Y je mˇeˇrena s konstantn´ı relativn´ı smˇerodatnou odchylkou (nejistotou) sr (Y ), m˚ uˇzeme vyjádˇrit ”váˇzen´ y” pr˚ umˇer n X Yi 2 s (Yi ) , (7) Y w = i=1 n X 1 s2 (Yi ) i=1 a rozptyl tohoto váˇzeného pr˚ umˇeru vypoˇcteme z s (Y w ) =

n

1

2

n X i=1

1 s2 (Yi )

1 X (Yi − Y w )2 . n − 1 i=1 s2 (Yi )

(8)

V pˇredchoz´ım odd´ıle byly zavedeny pojmy pˇresnost, správnost a nejistota. Vztah mezi smˇerodatnou odchylkou a pˇresnost´ı je ten, ˇze smˇerodatná odchylka je m´ırou pˇresnosti v´ ysledku mˇeˇren´ı. Taktéˇz nejistota je podle definice z pˇredchoz´ıho odd´ılu rovna hodnotˇe smˇerodatné odchylky v´ ybˇerového pr˚ umˇeru. Nejlépe bude v´ yznam v´ ybˇerov´ ych funkc´ı patrn´ y z následuj´ıc´ıch pˇr´ıklad˚ u, ve kter´ ych pro zjednoduˇsen´ı a zkrácen´ı je vypuˇstˇeno slovo ”v´ ybˇerov´ y”.

1.2

Aplikace v´ ybˇ erov´ ych charakteristik - I

Mˇeˇr´ıme plynovou chromatografi´ı koncentrace N-methyl-pyrrolidonu(NMP) ve vodˇe. Jedná se o velmi zˇredˇen´ y roztok (xi < 0.001) a provád´ıme 12 anal´ yz (nástˇrik˚ u) tak, abychom mohli provést statistické vyhodnocen´ı namˇeˇren´ ych dat. Detekce NMP prob´ıhá pomoc´ı plamenovˇe ionizaˇcn´ıho detektoru a v následuj´ıc´ı tabulce 1 jsou vyneseny hodnoty ploch p´ık˚ u 3 ˇ ’astn´ Last change 3.12.2009 Milan Bernauer & Bohumil Bernauer & Petr St y


n A n A

1 26805 7 26876

2 27712 8 27377

3 27587 9 27014

4 27536 10 26484

5 27001 11 26413

6 26721 12 26789

Tabulka 1: Hodnoty ploch p´ık˚ uAu ´mˇern´ ych koncentraci NMP ve vodˇe. Au ´mˇern´ ych koncentraci NMP ve vodˇe. Z tˇechto hodnot budeme cht´ıt zjistit pr˚ umˇernou yzy (resp. smˇerodatnou odchylku) a smˇerodatnou hodnotu A, rozptyl chromatografické anal´ odchylku pr˚ umˇerné hodnoty A. Pr˚ umˇernou hodnotu A vypoˇcteme ze vztahu 1, A = 27026. Smˇerodatnou odchylku chromatografického mˇeˇren´ı, respektive smˇerodatnou odchylku ”jednoho mˇeˇren´ı” vypoˇcteme ze vztah˚ u 2 a 3, s(A) = 433. A smˇerodatnou odchylku vypoˇcteného pr˚ umˇeru A z´ıskáme ze vztahu 4, s(A) = 125. Závˇery, které m˚ uˇzeme uˇcinit na základˇe pˇredchoz´ı statistiky jsou, ˇze relativn´ı nejistota jedné chromatografické anal´ yzy (vztah 5)je pˇribliˇznˇe sr (A) =

s(A) = 433/27026 = 0.016 A

(9)

a z´ıskaná pr˚ umˇerná hodnota plochy p´ıku A (´ umˇerná koncentraci NMP ve vodˇe) byla stanovena s relativn´ı nejistotou (vztah 6) sr (A) =

1.3

s(A) = 125/27026 = 0.005 A

(10)

Aplikace v´ ybˇ erov´ ych charakteristik - II

Metodou diferenciáln´ı destilace byl opakovanˇe stanovován limitn´ı aktivitn´ı koeficient γ1∞ N-methylformamidu (NMF) ve vodˇe5 . Vˇsechna mˇeˇren´ı byla provedena za konstantn´ı teploty 70o C, tlaku a pr˚ utoku stripovac´ıho plynu, s c´ılem stanovit nejistotu (smˇerodatnou odchylku) stanoven´ı γ1∞ touto metodou. Jelikoˇz kaˇzdá hodnota γ1∞ mˇela jinou váhu (rozptyl), plynouc´ıch z rozd´ıln´ ych nejistot pˇri chromatografické anal´ yze z´ıskan´ ych vzork˚ u, byl ∞ ∞ k v´ ypoˇctu pr˚ umˇeru pouˇzit vztah 7. Hodnoty γ1 a jejich nejistoty s(γ1 ) jsou uvedeny y v následuj´ıc´ı tabulce 2. Neváˇzen´ y pr˚ umˇer byl vypoˇc´ıtán z rovnice 1, γ1∞ = 1.52, váˇzen´ ∞ pr˚ umˇer z rovnice 7, γ1 w = 1.42 a smˇerodatná odchylka váˇzeného pr˚ umˇeru (nejistota) ze ∞ vztahu 8, respektive z jeho druhé odmocniny s(γ1 w ) = 0.05. V´ ysledky jsou vyneseny na obrázku 2. 5

Bernauer, M. and Dohnal, V. Temperature Dependence of Air-Water Partitioning of N-Methylated (C1 and C2) Fatty Acid Amides. J. Chem. Eng. Data. 2008, 53 , 2622 − 2631



n γ1∞ s(γ1∞ ) n γ1∞ s(γ1∞ ) n γ1∞ s(γ1∞ )

1 1.42 0.07 10 1.43 0.09 19 1.25 0.07

2 1.57 0.13 11 1.05 0.11 20 1.57 0.11

3 1.30 0.08 12 1.16 0.08 21 1.58 0.09

4 1.21 0.11 13 1.43 0.10 22 1.98 0.13

5 1.89 0.15 14 2.01 0.16 23 1.90 0.15

6 1.37 0.08 15 1.29 0.08 24 1.38 0.08

7 1.78 0.12 16 1.11 0.07 25 1.19 0.07

8 1.72 0.14 17 1.78 0.14 26 1.62 0.13

9 1.69 0.08 18 1.88 0.11

Tabulka 2: Namˇeˇrené hodnoty γ1∞ NMF ve vodˇe spoleˇcnˇe s hodnotami nejistot s(γ1∞ ), stanoven´ ymi na základˇe zákona ˇs´ıˇren´ı chyb.

Obrázek 2: Váˇzen´ y pr˚ umˇer opakovan´ ych mˇeˇren´ı limitn´ıho aktivitn´ıho koeficientu γ1∞ NMF ve vodˇe pˇri 70o C. , experimentáln´ı data γ1∞ ; − neváˇzen´ y pr˚ umˇer γ1∞ (rovnice 1); − − − váˇzen´ y pr˚ umˇer γ1∞ w (rovnice 7); − · − odlehlé hodnoty γ1∞ > γ1∞ ± 3s(γ1∞ ).

@

2

Z´ akon ˇ s´ıˇ ren´ı chyb

Mˇeˇr´ıme-li (pozorujeme) náhodnou veliˇcinu Y v závislosti na jin´ ych, taktéˇz náhodn´ ych veliˇcinách A1 , A2 aˇz An Y = f (Ai ), i = 1, . . . , n, (11)



m˚ uˇzeme odvodit vztah, kter´ y vystihuje, jak se nejistota plynouc´ı z náhodného charakteru veliˇcin Ai prom´ıtá do vypoˇctené hodnoty Y 6 n X n X ∂f ∂f cov(Ai , Aj ), s (Y ) = ∂Ai ∂Aj i=1 j=1 2

(12)

kde cov(Ai , Aj ) je prvek kovarianˇcn´ı matice C, která je definována 

cov(A1 , A1 ) . . . cov(A1 , An ) | {z }

 s2 (A1 )   .. .. ... C= . .   cov(An , A1 ) . . . cov(An , An ) {z } |

    .  

s2 (An )

Kovariance mezi dvˇema vektory Ai a Aj vyjadˇruje vztah cov(Ai , Aj ) = E{[Ai − E(Ai )][Aj − E(Aj )]} = E(Ai Aj ) − E(Ai )E(Aj ),

(13)

kde E() znaˇc´ı stˇredn´ı hodnotu. Z tohoto vztahu plyne pro kovarianci cov(Ai , Ai ) = s2 (Xi ). Pokud jsou Ai a Aj nezávislé, plat´ı cov(Ai , Aj ) = 0.

(14)

Opˇet si vˇse vysvˇetl´ıme na konkrétn´ım pˇr´ıpadu. Chceme stanovit objemov´ y pr˚ utok dus´ıku membránov´ ym modulem pomoc´ı bublinkového pr˚ utokomˇeru. Zároveˇ n chceme zjistit jeho smˇerodatnou odchylku. Pr˚ utok je nastavován nenakalibrovan´ ym hmotnostn´ım pr˚ utokoo mˇerem. Cela je pˇresnˇe temperována na 25 C, ale mˇeˇren´ı pr˚ utoku je provádˇeno na urˇcitém m´ıstˇe v laboratoˇri, kde teplomˇer, mˇeˇr´ıc´ı s pˇresnost´ı na 0.5o C udává teplotu 20o C. Bylo provedeno 12 mˇeˇren´ı pˇri konstantn´ım pr˚ utoku dus´ıku. Objem bublinkového pr˚ utokomˇeru je 10 ml. Tlak v laboratoˇri je 100.0 kPa. Namˇeˇrené hodnoty jsou uvedeny v tabulce 3. Uˇzit´ım n t/s n t/s

1 11.96 7 10.42

2 11.25 8 10.25

3 11.15 9 10.27

4 12.81 10 10.52

5 11.25 11 12.03

6 12.82 12 10.38

Tabulka 3: Namˇeˇrené hodnoty ˇcasu v sekundách pˇri mˇeˇren´ı pr˚ utoku bublinkov´ ym pr˚ utokomˇerem 6

Tento vztah lze odvodit z Taylorova rozvoje v okol´ı stˇredn´ıch hodnot Ai . Detaily tohoto odvozen´ı jsou k nalezen´ı napˇr´ıklad v: Meloun, M., Militk´ y, J., Chemometrie-Zpracován´ı Experimentáln´ıch Dat na IBM-PC, SNTL Praha 1990



vztah˚ u z pˇredchoz´ıho odd´ılu zjist´ıme pr˚ umˇernou hodnotu ˇcasu a smˇerodatnou odchylku pr˚ umˇerného ˇcasu (vztah 4). t = t ± s(t) = 11.26 ± 0.28 s. Pr˚ utok vypoˇcteme z následuj´ıc´ıho vztahu (zanedbáme korekci na objem vodn´ı páry uvolnˇené z náplnˇe bublinkového pr˚ utokomˇeru7 a pˇredpokládáme ideáln´ı chován´ı plynné fáze) VB T2 V˙ (T2 ) = , t T1

(15)

kde V˙ (T2 ) je objemov´ y pr˚ utok za teploty T2 = 298 K, VB je objem pr˚ utokomˇeru, T1 je absolutn´ı teplota, za které bylo provedeno mˇeˇren´ı pr˚ utoku a t je zmˇeˇren´ y ˇcas. Veliˇciny T1 a t nejsou známy s ”absolutn´ı” pˇresnost´ı, ale známe jejich nejistoty s(T1 ) a s(t). Také v´ıme, ˇze tyto dvˇe náhodné veliˇciny jsou mezi sebou nekorelované (neexistuje vztah mezi nimi, viz. vztah 14) a proto m˚ uˇzeme vztah 12 zjednoduˇsit na 2 np X ∂f s2 (Ai ), s (Y ) = ∂A i i=1 2

(16)

a pro náˇs konkrétn´ı pˇr´ıpad s (V˙ 2 ) = 2

2 2 V1 T2 1 V1 T2 1 2 − s (t) + − s2 (T1 ). T1 t2 t T12

(17)

Dosazen´ım ˇc´ıseln´ ych hodnot do pˇredchoz´ıho vztahu a pouˇzit´ım vztahu 3 vypoˇcteme smˇerodatnou odchylku pr˚ utoku celou s(V˙ 2 ) = 1.4 ml min−1 . V´ ysledek zap´ıˇseme V˙ 2 = 54.9 ± 1.4 ml min−1 .

3

Regresn´ı anal´ yza

Tentokráte se budeme zab´ yvat závislost´ı náhodné veliˇciny Y na promˇenné x, která nen´ı náhodnou veliˇcinou. Obecnˇe x m˚ uˇze b´ yt n-rozmˇerné. Závislost´ı Y = f (x) mysl´ıme matematick´ y model, kter´ ym se snaˇz´ıme popsat nˇejak´ y pozorovan´ y jev a z pozorovan´ ych hodnot 7

Korekce pr˚ utoku na objem vypaˇrené vodn´ı páry je V˙ s = V˙

p , (p + psT2 )

kde V˙ s je pr˚ utok suchého dus´ıku, V˙ je celkov´ y pr˚ utok plynu zmˇeˇren´ y bublinkov´ ym pr˚ utokomˇerem, p je celkov´ y tlak a psT2 je tlak nasycené vodn´ı páry pˇri teplotˇe T2 . Tlak nasycené vodn´ı páry pˇri 25o C je 3.16 kPa coˇz je pˇribliˇznˇe 3% z celkového tlaku a nav´ıc z d˚ uvodu konstrukˇcn´ıho uspoˇrádán´ı pr˚ utokomˇeru lze stˇeˇz´ı pˇredpokl´ adat rovnov´ aˇzné nasycen´ı proud´ıc´ıho plynu vodn´ı parou. Proto tuto korekci m˚ uˇzeme zanedbat.



odhadnout parametry tohoto matematického modelu. Napˇr´ıklad hodnotu tlaku nasycen´ ych par nad kapalinou (pozorovaná veliˇcina) v závislosti na teplotˇe dobˇre vystihuje v rozmez´ı od trojného bodu po normáln´ı bod varu známá Antoinova rovnice ln ps = A −

B , C +T

(18)

kde A, B a C jsou nastavitelné parametry, které z´ıskáme (v tomto pˇr´ıpadˇe nelineárn´ı) ˇ e teoreticky by staˇcily k urˇcen´ı tˇechto tˇr´ı parametr˚ regres´ı zmˇeˇren´ ych dvojic dat ps a T . Cistˇ u s pouze tˇri dvojice hodnot p a T . Nen´ı tˇeˇzké si pˇredstavit jak mizerná by byla správnost takto stanovené teplotn´ı závislosti ps . Proto se provád´ı mˇeˇren´ı v´ıce a je obecnˇe znám´ ym pravidlem, ˇze mˇeˇren´ı (n) má b´ yt pˇrinejmenˇs´ım o dvˇe v´ıce neˇz je parametr˚ u (np ) v modelové rovnici n = np + 2. Regresn´ı anal´ yza potom spoˇc´ıvá v nalezen´ı takov´ ych hodnot parametr˚ u A, B, C, tak aby po dosazen´ı do Antoinovy rovnice co ”nejlépe” popsaly zmˇeˇrenou závislost. Slovo ”nejlépe” je velmi d˚ uleˇzité a vede nás k urˇcen´ı nˇejakého kritéria, pomoc´ı kterého objektivnˇe a nestrannˇe urˇc´ıme ty ”nejlepˇs´ı” hodnoty parametr˚ u. Bohuˇzel nestaˇc´ı namˇeˇrená data ps (T ) vynést do grafu a proloˇzit je kˇrivkou odpov´ıdaj´ıc´ı Antoinovˇe rovnici s náhodnˇe zvolen´ ymi parametry, a spokojit se vizuáln´ım ”oko-metrick´ ym” vyhodnocen´ım a s konstatován´ım, ˇze kˇrivka pˇeknˇe ”sed´ı” na namˇeˇren´ ych datech. T´ımto objektivn´ım kritériem je v pˇr´ıpadˇe veliˇcin s normáln´ım rozdˇelen´ım 8 suma ˇctverc˚ u rozd´ıl˚ u mezi zmˇeˇrenou a vypoˇctenou hodnotou Y respektive n X S(A, B, C) = (psi,zmˇeˇrená − psi,vypoˇctená )2 .

(19)

i=1

Tato funkce se naz´ yvá objektivn´ı funkce, nˇekdy téˇz c´ılová nebo kriteriáln´ı funkce. Hodnoty parametr˚ u A, B, C pˇri kter´ ych nab´ yvá S(A, B, C) minima jsou potom jejich nejlepˇs´ım nestrann´ ym odhadem pro dan´ y soubor experimentáln´ıch dat. V´ yˇse zm´ınˇen´ ymi slovy jsou tˇemi ”nejlepˇs´ımi”:-). ´ Ukoly regresn´ı anal´ yzy se daj´ı shrnout do následuj´ıc´ıho seznamu: 1. Stanoven´ı hodnot parametr˚ u matematického modelu 2. Odhad nejistoty vypoˇcten´ ych parametr˚ u 3. Statistické vyhodnocen´ı kvality korelace goodnes-of-fit

3.1

Line´ arn´ı regrese

Nyn´ı se budeme zab´ yvat pˇr´ıpadem, kdy modelová rovnice je lineárn´ı v˚ uˇci parametr˚ um, jako napˇr´ıklad v následuj´ıc´ım polynomu y(x) = a1 + a2 x + a3 x2 +, . . . , +anp xnp −1 , 8

Odvozen´ı je napˇr´ıklad uvedeno v: Bard, Y., Nonlinear Parameter Estimation, Academic Press, New York 1973



y

yk y(xk )

xk

x1

xn

x

u

Obrázek 3: Konstrukce objektivn´ı funkce a princip metody nejmenˇs´ıch ˇctverc˚ u. jsou experimentáln´ı (zmˇeˇrené) body, kter´ ym odpov´ıdá hodnota yk pˇri xk , n je celkov´ y poˇcet pozorován´ı, y(xk ) je vypoˇctená hodnota z modelové rovnice pro xk . Plochy jednotliv´ ych 2 ˇctverc˚ u jsou rovny (yk − y(xk )) . Plná kˇrivka reprezentuje proloˇzen´ı bod˚ u zvolenou modelovou rovnic´ı. nebo v kompaktnˇejˇs´ı formˇe y(x) =

np X

ai xi−1 ,

(20)

i=1

kde np je poˇcet parametr˚ u. Pro tuto rovnici sestav´ıme kriteriáln´ı funkci minima souˇctu ˇctverc˚ u odchylek n X S(~a) = (yk − y(xk ))2 , (21) k=1

nebo Sw (~a) =

n X (yk − y(xk ))2

s2 (yk )

k=1

,

(22)

2

kde s (yk ) je rozptyl k-tého experimentáln´ıho bodu, yk je hodnota náhodné veliˇciny a n je poˇcet bod˚ u. Pro pochopen´ı jednotliv´ ych veliˇcin je na obrázku 3 znázornˇena grafická konstrukce objektivn´ı funkce a jej´ı v´ yznam v metodˇe nejmenˇs´ıch ˇctverc˚ u. Kvalita v´ ysledné korelace je vystiˇzena veliˇcinou naz´ yvanou reziduáln´ı rozptyl s2d =

S(~a) , n − np

nebo

s2d =

Sw (~a) , n − np

(23)

kde ~a je vektor parametr˚ u. Jej´ı odmocnina se naz´ yvá smˇerodatná odchylka korelace q sd = s2d . (24) 9 ˇ ’astn´ Last change 3.12.2009 Milan Bernauer & Bohumil Bernauer & Petr St y


Minimalizac´ı rovnice 21 nebo 22 se vyhodnot´ı parametry ai rovnice 20. Minimalizace spoˇc´ıvá v ˇreˇsen´ı tzv. soustavy normáln´ıch rovnic, které vzniknou derivac´ı S(~a) podle parametr˚ u ai ∂S(~a) = 0, i = 1, . . . , np . (25) ∂ai Pokud zap´ıˇseme vektor pozorovan´ ych hodnot jako ~y = (y1 , . . . , yn )T a vektor parametr˚ u T ~a = (a1 , . . . , anp ) , m˚ uˇzeme zm´ınˇenou soustavu normáln´ıch rovnic zapsat maticovˇe FFT ~a = F~y ,

(26)

kde F je matice prvn´ıch derivac´ı modelové funkce podle parametr˚ u   ∂f1 ∂fn ...  ∂a1 ∂a1   . ..  . . .. F= .  ,  .  ∂f1 ∂fn  ... ∂anp ∂anp kde fk =

np X

ai xki−1 ,

k = 1, . . . , n.

(27)

i=1

Napˇr´ıklad pro rovnici 20 bude se daj´ı vypoˇc´ıtat z9

∂fk ∂a1

= 1, pro vˇsechna k = 1, . . . , n. Neznámé parametry ai ~a = (FFT )−1 F~y .

(28) T

Souˇcin matice F a transponované matice FT (respektive FF ) je v pˇr´ıpadˇe lineárn´ı funkce 20 roven následuj´ıc´ı matici   1 ∂ 2 S(~a) 1 ∂ 2 S(~a) ...  2 ∂a ∂a 2 ∂a1 ∂anp  1 1     . ..     2   1 ∂ S(~a) .. .. T   FF =  . .  2 ∂a ∂a i j     .. .     2 2 1 ∂ S(~a)   1 ∂ S(~a) ... 2 ∂anp ∂a1 2 ∂anp ∂anp coˇz je regulárn´ı matice typu (np , np ). Souˇcin reziduáln´ıho rozptylu a inverzn´ı matice Ca = s2d (FFT )−1

(29)

9

Vˇetˇsinou se pˇri v´ ypoˇctu parametr˚ u ai dává pˇrednost ˇreˇsen´ı rovnice 26 , pˇred ˇreˇsen´ım rovnice 28 spojené s inverz´ı matice FFT , coˇz je v´ ypoˇcetnˇe n´ aroˇcnˇejˇs´ı u ´kon.



se naz´ yvá kovarianˇcn´ı matice parametr˚ u a. Vezmeme-li v u ´vahu pouze lineárn´ı rovnici (np = 2), to znamená ˇze máme pouze dva nastavitelné parametry a1 a a2 , soustava normáln´ıch rovnic bude m´ıt rozmˇery 2 × 2   n X xi   n   i=1 FFT =  X , n n X  2  xi xi i=1

i=1

pˇriˇcemˇz matice F bude m´ıt rozmˇer 2 × n, kde n je poˇcet mˇeˇren´ı (experimentáln´ıch bod˚ u) 1 ... 1 F= . x1 . . . xn V tomto pˇr´ıpadˇe bude rovnice 26 vypadat následovnˇe   n X   n x y i 1   1 . . . 1  ..    a1 i=1 =  X   . . n n X a2 x1 . . . xn  2  xi xi yn i=1

i=1

Prvky v kovarianˇcn´ı matici Ca Ca =

cov(a1 , a1 ) cov(a1 , a2 ) cov(a2 , a1 ) cov(a2 , a2 )

,

poskytuj´ı odhad rozptylu s2 (ai ) = cov(ai , ai ) vypoˇcten´ ych parametr˚ u a1 , a2 a kovarianci mezi parametry cov(a1 , a2 ) = cov(a2 , a1 ). Samozˇrejmˇe v dobˇe, kdy je k dispozici nepˇreberné mnoˇzstv´ı softwarového vybaven´ı (Matlab, Excel, Gnuplot, Origin . . . atd.) nemus´ı se ”ruˇcnˇe” provádˇet v´ ypoˇcty souvisej´ıc´ı s ˇreˇsen´ım soustavy normáln´ıch rovnic. Ale málokter´ y z v´ yˇse uveden´ ych nástroj˚ u poskytuje v´ ypoˇcet kovarianˇcn´ı matice, pomoc´ı které se m˚ uˇze vyhodnotit nejistota ve vypoˇcten´ ych parametrech, a proto ˇcasto ˇreˇsen´ı tohoto problému padá na bedra nebohého uˇzivatele.

3.2

V´ aˇ zen´ a line´ arn´ı regrese

V pˇredchoz´ım odd´ıle se tiˇse pˇredpokládalo, ˇze pozorované veliˇciny y maj´ı shodné rozptyly σ 2 (y1 ) = σ 2 (y2 ) =, . . . , = σ 2 (yn ). Jin´ ymi slovy, ˇze vˇsechny hodnoty y maj´ı stejnou váhu. Pokud bude tˇreba ”penalizovat” urˇcité hodnoty pozorovan´ ych hodnot, tˇreba z d˚ uvodu zhorˇsen´ ych experimentáln´ıch podm´ınek (pˇribl´ıˇzen´ı se k hranici citlivosti aparatury), je tˇreba pouˇz´ıt metodu váˇzené regrese, která je pouh´ ym zobecnˇen´ım postupu popsaného v pˇredchoz´ım odd´ıle. Tam byla pro odhad parametr˚ u lineárn´ıho modelu odvozena rovnice 28 ~a = (FFT )−1 F~y .

(28)



Pro zobecnˇen´ı v pˇr´ıpadˇe nestejn´ ych rozptyl˚ u σ 2 (yi ) je zavedena kovarianˇcn´ı matice Cy = σ 2 (y)W,

(30)

kde W je matice vah jednotliv´ ych pozorovan´ ych hodnot. Séri´ı u ´prav10 lze odvodit vztah pro ˇreˇsen´ı soustavy normáln´ıch rovnic vedouc´ı k v´ ypoˇctu parametr˚ u lineárn´ıho modelu T

~a = (FW−1 F )−1 FW−1 ~y .

(31)

Kovarianˇcn´ı matice Cy parametr˚ u ~a je potom T

Cy = s2d (FW−1 F )−1 ,

(32)

kde s2d je reziduáln´ı rozptyl vypoˇcten´ y podle rovnice 23.

3.3

Urˇ cen´ı poˇ ctu parametr˚ u modelov´ e rovnice

Problém urˇcen´ı poˇctu parametr˚ u rovnice 20 potˇrebn´ ych k popisu experimentáln´ıch dat m˚ uˇze b´ yt vyˇreˇsen pouˇzit´ım statistického F-testu (Fisher˚ uv F-test) o rovnosti rozptyl˚ u dvou základn´ıch soubor˚ u, kdy se testuj´ı reziduáln´ı rozptyly dvojice fit˚ u s r˚ uzn´ ym poˇctem pouˇzit´ ych parametr˚ u. Problém je moˇzno formulovat tak, ˇze se snaˇz´ıme rozhodnout zda reziduáln´ı rozptyl fitu s np,1 parametry je vˇetˇs´ı neˇz reziduáln´ı rozptyl fitu s np,2 parametry na zadané hladinˇe v´ yznamnosti α. V ˇreˇci matematické statistiky vyslov´ıme nulovou hypotézu (H0), ˇze nedoˇslo ke sn´ıˇzen´ı s2d a tuto budeme testovat proti alternativn´ı hypotéze (H1), ˇze doˇslo k sn´ıˇzen´ı s2d H0 : s2d,1 ≤ s2d,2 ; H1 : s2d,1 > s2d,2 . (33) Reziduáln´ı rozptyly s2d,1 a s2d,2 odpov´ıdaj´ıc´ı regresi pomoc´ı rovnice 20 s np,1 respektive np,2 parametry a n1 respektive n2 poˇctem experimentáln´ıch bod˚ u se vypoˇc´ıtaj´ı Pn1 2 i=1 (Yi,1 − Yi,1 ) s2d,1 = , (34) n1 − np,1 Pn2 2 2 i=1 (Yi,2 − Yi,2 ) sd,2 = . (35) n2 − np,2 Testovac´ı statistika je popsána funkc´ı s2d,1 R= 2 . sd,2

(36)

Rozhodnut´ı zda zv´ yˇsen´ım poˇctu parametr˚ u rovnice 20 z np,1 na np,2 z´ıskáme menˇs´ı reziduáln´ı rozptyl prob´ıhá na základˇe nerovnosti R > Fα (n1 − np,1 , n2 − np,2 ) ⇒

zam´ıtáme H0, tud´ıˇz s2d,1 > s2d,2 ,

R < Fα (n1 − np,1 , n2 − np,2 ) ⇒

nem˚ uˇzeme zam´ıtnout H0, tud´ıˇz plat´ı H1 s2d,1 ≤ s2d,2 ,

kde Fα (n1 −np,1 , n2 −np,2 ) je tabelovaná hodnota Fisherova rozdˇelen´ı pro dan´ y poˇcet stupˇ n˚ u volnosti a dané hladinˇe v´ yznamnosti. 10

ˇ Viz napˇr´ıklad skripta: Jaroˇs, F. a spol., Pravdˇepodobnost a Statistika, VSCHT Praha 2002



3.4

χ2 test konsistence sd a s(yi )

Pokud jsou hodnoty smˇerodatn´ ych odchylek experimentáln´ıch veliˇcin s(yi ) statisticky konzistentn´ı s v´ yslednou hodnotou souˇctu ˇctverc˚ u váˇzen´ ych odchylek Sw (~a), mus´ı se hodnota váˇzené smˇerodatné odchylky korelace sd pohybovat v okol´ı jedniˇcky s Sw (~a) sd = ≈ 1, (37) n − np kde n je poˇcet experimentáln´ıch bod˚ u a np je poˇcet parametr˚ u. Toto kritérium je moˇzno povaˇzovat za splnˇené, pokud hodnota Sw (~a) leˇz´ı v intervalu vymezeném kvantily rozdˇelen´ı χ2 pro dan´ y poˇcet stupˇ n˚ u volnosti (n − np ) na zvolené hladinˇe v´ yznamnosti α χ2α/2 (n − np ) < Sw (~a) < χ21−α/2 (n − np ),

(38)

kdy hladina v´ yznamnosti je obvykle volena α = 0.05. Pokud je hodnota Sw (~a) niˇzˇs´ı neˇz doln´ı kvantil χ2α/2 (n − np ) jsou nejistoty pˇr´ısluˇs´ıc´ı experimentáln´ım bod˚ um ”nadhodnocené” oproti v´ ysledné hodnotˇe sd . Slovy terminologie zavedené v u ´vodu tohoto dokumentu se pravdˇepodobnˇe jedná o správná ale nepˇresná data. um ”podhodPokud Sw (~a) > χ21−α/2 (n − np ) jsou nejistoty pˇr´ısluˇs´ıc´ı experimentáln´ım bod˚ nocené” a data jsou pravdˇepodobnˇe pˇresná ale nesprávná.

Line´ arn´ı regrese - nev´ aˇ zen´ a Rozˇs´ıˇr´ıme pˇr´ıklad z pˇredchoz´ıho odd´ılu (mˇeˇren´ı pr˚ utoku bublinkov´ ym pr˚ utokomˇerem) o kalibraci hmotnostn´ıho regulátoru pr˚ utoku. Budeme mˇeˇrit pr˚ utok plynu (V˙ B ) pro r˚ uzné hodnoty nastavené na regulátoru (VD ) a z tˇechto uspoˇrádan´ ych dvojic se budeme snaˇzit vyhodnotit parametry modelové rovnice. Pouˇzijeme rozvoj rovnice 20 do prvn´ıho stupnˇe, resp. lineárn´ı závislost y(x) = a1 + a2 x, (39) kde x = VD a y = V˙ B . Postup mˇeˇren´ı pr˚ utoku je stejn´ y jako v pˇredchoz´ım pˇr´ıpadˇe. V tabulce 4 jsou uvedeny namˇeˇrené hodnoty a pˇrepoˇctené hodnoty za normáln´ıch podm´ınek (p = 101.325 kPa a T2 = 273.15 K) Pˇrepoˇcet na normáln´ı podm´ınky byl proveden podle rovnice 15 se zapoˇcten´ım vlivu zmˇeny tlaku (z atmosférického p1 na standardn´ı p2 ) V1 (T2 ) T2 p1 . V˙ (T2 ) = t T1 p2

(40)

Smˇerodatná odchylka s(V˙ ) byla vypoˇctena ze vztahu 16 zahrnuj´ıc´ı pˇr´ıspˇevky s(t) s(T ) a s(p1 )11 . Lineárn´ı regresi je moˇzno provést napˇr´ıklad v tabulkovém procesu (MS Excel, OpenOffice Calc, Gnumeric), kter´ y podporuje maticové operace (maticov´ y souˇcin a operaci 11

Atmosférick´ y tlak byl mˇeˇren s pˇresnost´ı 0.1kPa



VD / (ml min−1 ) VB / ml T1 /o C n 1 2 3 4 5 6 7 8 9 10 11 12 t/s s(t)/s V˙ B s(V˙ B )

50 10 23 6.90 7.22 7.21 7.6 7.2 7.31 7.05 7.22 6.92 7.25 7.25 7.28

40 10 23

30 10 25

9.2 9.06 9.16 9.5 9.12 9.03 9.16 9.12 9.4 9.5 9.1 9.12

20 4 25 ˇ t/s Cas 6.13 6.23 6.35 6.04 6.20 6.22 6.13 6.16 6.13 6.16 6.26 6.22

10 2 25

11.50 5.74 11.60 6.08 11.59 5.93 11.46 5.58 11.49 5.62 11.54 5.68 11.8 5.86 11.43 5.52 11.63 5.33 11.78 5.86 11.67 5.49 11.67 5.42 Statistika 7.20 9.26 11.6 6.19 5.68 0.05 0.05 0.04 0.02 0.07 −1 ˙ Pr˚ utok VB /(ml.min ) pˇri 273.15 K a 101.325 kPa 76.90 60.03 47.41 35.55 19.40 0.57 0.31 0.14 0.13 0.22

15 4 25

12 2 25

10.09 10.19 10.31 10.28 10.12 10.28 10.22 10.15 10.18 10.25 10.25 10.19

5.75 5.63 5.62 5.75 5.75 5.68 5.69 5.72 5.59 5.62 5.62 5.68

10.21 0.02

5.68 0.02

21.54 0.04

19.37 0.06

Tabulka 4: Tabulka experimentáln´ıch hodnot mˇeˇren´ı pr˚ utoku bublinkov´ ym pr˚ utokomˇerem. VD je hodnota udávaná na hmotnostn´ım pr˚ utokomˇeru, VB je objem bublinkového pr˚ utokomˇeru, T1 je teplota pˇri které bylo provedeno mˇeˇren´ı pr˚ utoku (V˙ B )

inverze matice). Tento postup je jednoduch´ y ale zároveˇ n neefektivn´ı. Pokud bude zmˇenˇen poˇcet parametr˚ u nebo pozorován´ı bude tˇreba provést ”pˇreprogramován´ı” celého v´ ypoˇctu, protoˇze tyto tabulkové procesory neumoˇzn ˇuj´ı dynamickou alokaci pol´ı. Jinou moˇznost´ı je pouˇz´ıt nˇekterého z programovac´ıch jazyk˚ u, at’ jiˇz kompilovaného - FORTRAN, C, nebo interpretovaného Matlab, Maple. V´ ysledky byly v tomto pˇr´ıpadˇe z´ıskány pomoc´ı programu WLINREG12 napsaném v programovac´ım jazyce FORTRAN. V´ ystup z v´ ypoˇctu je uveden v dodatku. V´ ysledné hodnoty (odhady!) parametr˚ u rovnice 39 jsou a1 = 2.938

a2 = 1.467

Souˇcet ˇctverc˚ u odchylek a smˇerodatná odchylka korelace (rovnice 24) jsou S(~a) = 29.98 12

sd = 2.449

Zdrojov´ y k´ od programu je k dispozici na adrese ...



Smˇerodatné odchylky stanoven´ ych parametr˚ u a1 , a2 se urˇc´ı z kovarianˇcn´ı matice (rovnice 29) 3.6079 −0.1088 Ca = , −0.1088 0.0043 kde prvky na hlavn´ı diagonále jsou rozptyly parametr˚ u, které po odmocnˇen´ı dávaj´ı smˇerodatné odchylky s(a1 ) = 1.899 s(a2 ) = 0.066. Pokud bude tˇreba stanovit nejistotu v pr˚ utoku pˇri nastavené hodnotˇe VD = 30, pouˇzijeme rovnici 12 v podobˇe np np X X ∂Vb ∂Vb 2 s (Vb ) = cov(ai , aj ) ∂a ∂a i j i=1 j=1 = 1 cov(a1 , a1 ) + x2 cov(a2 , a2 ) + 2x cov(a1 , a2 ),

kde za cov(ai , aj ) dosad´ıme hodnoty z vypoˇctené kovarianˇcn´ı matice. Po dosazen´ı obdrˇz´ıme Pro VD = 30 je pr˚ utok Vb = 46.94 ± 0.98 ml min−1 Na obrázku 4 jsou znázornˇeny zmˇeˇrené hodnoty a v´ ysledky lineárn´ı regrese.

Line´ arn´ı regrese - v´ aˇ zen´ a φ Zdánlivá molárn´ı tepelná kapacita Cp,1 N-methyl-pyrrolidonu (NMP) ve vodˇe byla mˇeˇrena na kalorimetru Pickerova typu v závislosti na molalitˇe roztoku pˇri atmosférickém tlaku a teplotˇe 298 K13 . C´ılem tohoto experimentu bylo stanoven´ı zdánlivé molárn´ı tepelné kapaφ.∞ city NMP ve vodˇe v nekoneˇcném zˇredˇen´ı Cp,1 , která je rovna parciáln´ı molárn´ı tepelné ∞ φ modekapacitˇe v nekoneˇcném zˇredˇen´ı C p,1 . Proto bylo potˇreba koncentraˇcn´ı závislost Cp,1 φ φ lovat lineárn´ı funkc´ı Cp,1 (b) = a1 + a2 b a z v´ ysledku urˇcit hodnotu Cp,1 v b = 0 (respektive urˇcit parametr a1 ). V´ ysledky experimentu jsou uvedeny v tabulce 5 Nejistota ve stanovené φ hodnotˇe Cp,1 vzr˚ ustala smˇerem k niˇzˇs´ım hodnotám molality roztoku, z d˚ uvodu niˇzˇs´ı odezvy mˇeˇr´ıc´ı aparatury a t´ım pádem vˇetˇs´ım vlivu ”ˇsumu” na v´ yslednou hodnotu. Proto byla zvolena metoda váˇzené lineárn´ı regrese, která zohledn´ı nejistoty v jednotliv´ ych mˇeˇren´ıch. K proloˇzen´ı dat byl pouˇzit jiˇz zm´ınˇen´ y lineárn´ı model φ (b) = a1 + a2 b, Cp,1

a nástrojem, kter´ ym byla povedena samotná regrese byl opˇet program WLINREG. Matice vah W obsahovala na své hlavn´ı diagonále rozptyly experimentáln´ıch hodnot a zbytek prvk˚ u byl roven nule. V´ ysledky jsou znázornˇeny na obrázku 5. 13 Bernauer, M. and Dohnal, V. Temperature dependences of limiting activity coefficients and Henry’s law constants for N-methylpyrrolidone, pyridine, and piperidine in water. Fluid Phase Equilib. 2009, 282 , 100 − 107



Obrázek 4: Kalibrace pr˚ utokomˇeru, VD u ´daj na pr˚ utokomˇeru, Vx ”reálná” hodnota pr˚ utoku. , experimentáln´ı body; −, kˇrivka lineárn´ı regrese; − − −, 95% konfidenˇcn´ı pásy.

p

n

b/mol kg−1

φ Cp,1 (exptl)

1 2 3 4 5 6 7 8

0.500 0.314 0.194 0.106 0.081 0.047 0.019 0.009

302 305 304 306 297 308 275 308

φ s(Cp,1 (exptl)) −1 J mol K−1 2 2 2 2 5 10 15 15

φ Cp,1 (calc)

304 304 304 305 305 305 305 305

wdev 0.3 −0.5 0.2 −0.7 1.5 −0.3 2.0 −0.2

Tabulka 5: Tabulka experimentáln´ıch hodnot a v´ ysledk˚ u váˇzené lineárn´ı regrese. φ φ φ (exptl), Cp,1 (calc) jsou experimentáln´ı, resp. vypoˇctené hodnoty a wdev = (Cp,1 (calc)− Cp,1 φ φ Cp,1 (exptl))/s(Cp,1 (exptl)) je váˇzená odchylka.



φ Obrázek 5: Váˇzená lineárn´ı regrese dat zdánlivé tepelné molárn´ı kapacity Cp,1 NMP ve vodˇe v závislosti na molalitˇe b. , experimentáln´ı body; −, kˇrivka lineárn´ı regrese; − − −, 95% konfidenˇcn´ı pásy.

p


Milan Bernauer, Bohumil Bernauer, Petr Št astný. dat. Ústav Anorganické Technologie. Obsah

Recommend Documents