Katedra teoretické a fyzikální chemie, Přírodovědecká fakulta MU, Brno
Nelineární regrese v chemické kinetice
Miroslav Holík
F 457/1999
2
Obsah
Úvod Část 1 - Regresní rovnice 1. Lineární regrese 2. Vícerozměrná lineární regrese 3. Nelineární regrese 4. Statistické testování regresních parametrů Část 2 – Integrované kinetické rovnice 1. Unimolekulární vratná reakce 2. Bimolekulární reakce 3. Enzymově katalyzovaná reakce Část 3 – Výpočtové programy 1. Nelin: exponenciální regrese s Taylorovým rozvojem 2. Simpl: exponenciální regrese s modifikovaným simplexem 3. Rad2a: nelineární regrese s Taylorovým rozvojem 4. Michme: nelineární regrese s Taylorovým rozvojem
3
Úvod V chemické kinetice se často setkáváme s rovnicemi, které, pokud jsou použity pro regresní odhad parametrů, jsou vzhledem k těmto parametrům nelineární a nelze tedy na ně aplikovat postupy lineární regrese. Nejčastější řešení spočívá v transformaci těchto rovnic na lineární tvar buď jejich inverzí nebo zlogaritmováním. To ovšem s sebou přináší problémy, které si běžný uživatel takové zlinearizované rovnice třeba vůbec neuvědomuje. Obecně lze regresní rovnici napsat takto: y = f(x,β) + ε
(1)
kde y jsou naměřené hodnoty tzv. vysvětlované neboli závislé promměnné, x jsou nastavené (zvolené) hodnoty vysvětlující neboli nezávislé proměnné, β jsou neznámé paramentry (koeficienty úměrnosti), které pomocí regrese odhadujeme tak, aby součet druhých mocnin reziduálů ε byl co nejmenší. Reziduály představují náhodné chyby při měření hodnot veličiny y a předpokládá se o nich, že mají konstantní rozptyl a nulovou střední hodnotu. (2)
E( ε2 ) = σ2
E( ε ) = 0
Každá nelineární transformace transformuje také tyto reziduály, což vede k tomu, že předpoklady pro regresi (t.j. rovnice 2) nejsou splněny. To se může projevit na výsledných hodnotách odhadů parametrů β, které nejsou z hlediska původní nelineární regrese optimální a hodnoty ycalc, vypočítané z hodnot x a těchto odhadů β, se mohou značně lišit od hodnot yexp získaných měřením. Proto je vhodnější volit pro odhad parametrů β raději některou z metod nelineární regrese. Některé z častěji používaných kinetických rovnic jsou uvedeny v tomto spisku.
4
ČÁST 1. REGRESNÍ ROVNICE 1. Lineární regrese V praxi se často setkáváme s případy, kdy měřená veličina je lineární funkcí jedné nebo více nastavitelných veličin; nejjednodušší případ nastane tehdy, je-li měřená veličina závislá na jedné veličině nezávislé a úkolem je stanovit míru této závislosti. yi = a + bxi + εi
(3)
Experimentálně tedy naměříme pro řadu hodnot x (i = 1 až n) řadu závislých proměnných y (i = 1 až n) a výpočtem určíme optimální hodnoty parametrů a a b tak, aby součet druhých mocnin reziduálů εi byl minimální. Vzhledem k tomu, že jednotlivá měření jsou zatížena chybou, je třeba vyhodnotit větší počet dvojic měření a to minimálně tři na jeden optimalizovaný parametr. Postup výpočtu: - pro hodnoty xi a yi vypočítáme aritmetické průměry xp a yp (nezaokrouhlujeme!). - hodnoty xi a yi centrujeme, tzn. od každé odečteme odpovídající průměr xic = xi - xp, yic = yi - yp. - vypočítáme součty součinů centrovaných hodnot Σ (xic)2, Σ (yic)2 a Σ (yic*xic) - parametry a a b pak počítáme dle těchto vztahů: b = Σ (yic*xic) / Σ (xic)2
(4) (5)
a
= yp - b*xp
- z parametrů a a b a nezávisle proměnných hodnot xi vypočítáme vyrovnané hodnoty Yi a reziduály εi jako rozdíly y-hodnot naměřených a vypočtených (6)
Yi = a + b*xi
εi = yi - Yi
- úspěšnost regresního modelu se testuje pomocí směrodatné odchylky regrese (standard error of estimate) se a pomocí korelačního koeficientu r; regrese je tím lepší, čím bližší je korelační koeficient číslu 1( resp. –1). se = √ [ (yi - Yi )2 / (n - 2)]
(7) (8)
5
r
= Σ (yic*xic)/ √[Σ (xic)2 * Σ (yic)2]
-1 ≤ r ≤ 1
nebo (9)
r
= Σ (yic*Yic)/ √[Σ (yic)2 * Σ (Yic)2] t
0≤r≤1
= √ [ r2 (n - 2) / (1 -r2)]
Korelační koeficient se musí testovat ve vztahu k počtu experimentálních dat.Obvykle očekáváme, že bude významný alespoň na 0.1%-ní hladině spolehlivosti, t.j. α = 0.001. Je-li hodnota t vypočítaná z r větší než kritická hodnota t-rozdělení pro odpovídající počet stupňů volnosti, t.j. ν = n-2, lze z 99.9%-ní pravděpodobností usoudit, že odpovídající lineární závislost nevznikla náhodou. - ze směrodarné odchylky regrese vypočítáme směrodatné odchylky parametrů aab (10) (11)
sa = sE*√ [ (1/n) + xp2 / Σ (xic)2 ] sb = sE*√ [ 1 / Σ (xic)2 ]
a testujeme podíly a/sa a b/sb pomocí kritických hodnot t-rozdělení. Je-li vypočtená hodnota větší než tabelovaná pro α = 0.05 a odpovídající počet stupňů volnosti, t.j. ν = n-2, pak se odpovídající parametr významně liší od nuly. Tabulka kritických hodnot t-rozdělení -----------------------------------------------------------------------------------------------ν α = 0.05 α = 0.001 ν α = 0.05 α = 0.001 −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− 2 4.303 22.326 12 2.179 3.930 3 3.182 10.213 13 2.160 3.852 4 2.776 7.173 14 2.145 3.787 5 2.571 5.893 15 2.131 3.733 6 2.447 5.208 16 2.120 3.686 7 2.369 4.785 17 2.110 3.646 8 2.306 4.501 18 2.101 3.610 9 2.262 4.297 19 2.093 3.579 10 2.228 4.144 20 2.086 3.552 11 2.201 4.025 ∞ 1.960 3.090 -----------------------------------------------------------------------------------------------ν = počet stupňů volnosti, α = hladina pravděpodobnosti
6
2. Vícerozměrná lineární regrese
Pokud je vysvětlovaná proměnná y lineárně závislá na více vysvětlujících proměnných x, např. (12)
y = b0 + b1x1 + b2x2 + ... + bp-1xp-1 + ε
pak se tento vztah nazývá vícerozměrná lineární regrese. Obecně lze tuto regresi vyjádřit v maticovém zápisu takto (13)
Ymn = Xmp. Bpn + Εmn
kde m je počet měření t.j. regresních vztahů mezi závisle proměnnou y a nezávisle proměnnými x), p je počet hledaných parametrů b regresní rovnice (t.j. počet vysvětlujících proměnných x zvětšených o jednu pro konstantní neboli lokalizační člen), n je počet vysvětlovaných proměnných Y (běžně je n=1, ale regresi několika závislých proměnných vysvětlovaných stejnými nezávislými proměnnými lze spojit v jednu rovnici, t.j. n>1) a E je matice reziduálů t.j. odchylek naměřených y hodnot od hodnot vysvětlených regresní závislostí ycalc. Parametry b se v metodě nejmenších čtverců odhadnou pomocí generalizované inverze matice X (14)
B = (X'X)-1X'Y
V případě silné korelace mezi jednotlivými vysvětlujícími proměnnými x, mohou při výpočtu B nastat problémy s inverzí matice X'X, která v takovém případě bude kvazisingulární a výsledek může být zatížen velkou chybou. Potom je lépe počítat B pomocí t.zv. pseudoinverze matice X. Každou matici (např. X) lze rozložit singulárním rozkladem na tři matice charakteristických vlastností: (15)
Xmp = Ump. Spp. (Vpp)'
kde U a V jsou orthonormální matice, pro které platí VV'=V'V=I a U'U=I a S je diagonální matice singulárních hodnot uspořádaných tak, že s11 > s22 > ... > spp. Pseudoinverze X je pak rovna (16)
(Xmp)+ = Vpp. (Spp) -1. (Ump)'
V případě korelace mezi vektory matice X se matice S upraví před inverzí tak, že se z ní vypustí velmi malé singulární hodnoty, které by po inverzi dávaly 7
příliš velkou váhu málovýznamným vektorům matic U a V. Stejným způsobem se zmenší i odpovídající rozměr u matic U a V. (17)
(X‡mp)+ = Vpq.(Sqq) -1.(Umq)'
q
Potom se odhad matice parametrů B‡ vypočte podle rovnice (18)
B‡pn = (X‡mp)+.Ymn
q
Matici parametrů B‡, získanou z dat Y a X (t.zv. tréninkový soubor), lze použít k předpovědi dat Y0 z nových dat X0 nezahrnutých do X, (t.zv. testovaný soubor). (19)
Y0 = X0.B‡
Pro zjištění vhodnosti modelu daného maticí X k vysvětlení experimentálních hodnot Y, t.j. k vypočtení teoretických hodnot Ycalc není třeba parametry B‡ počítat. Stačí použít t.zv. 'hat' matici H ke transformaci Y na Ycalc. (20)
Ycalc = X (X'X) -1X'.Y = H.Y
nebo
Ycalc = X.X+.Y = H.Y
nebo po dosazení za X ze singulárního rozkladu dostaneme (21)
Ycalc = U.S.V'.(V.S.U'.U.S.V') –1.V.S.U'.Y = = U.S.V'. (V.S-2.V' ) .V.S.U'.Y = U.U'.Y
nebo
Ycalc = X.X+.Y = U.S.V'. V.S-1.U'.Y = U.U'.Y
Ke stejnému výsledku lze dospět i tehdy, když se do regrese použijí místo korelovaných proměnných X orthogonální proměnné U (tento postup je podstatou t.zv. regrese s hlavními komponentami (PCR = principal component regression). (22)
Ycalc = U (U'U) -1U'.Y = U.U'.Y
Vypočítané hodnoty závisle proměnné Ycalc lze tedy získat projekcí experimentálních hodnot Y pomocí projekční matice vytvořené z orthogonálních vektorů po singulárním rozkladu matice X. Přitom se může použít matice Ump celá nebo podobně zredukovaná v rozměru p na q jako při výpočtu pseudoinverze X (rovnice 17).
8
Při metodě PCR se počítá regrese hodnot Y na hlavních komponentách (skórech) C, které jsou orthogonální, ale nejsou normované (normou je odpovídající singulární hodnota) a jsou zredukované v rozměru p na q. (23)
Ycalc = Xmp.Bpn = Ump.Spp.(Vpp)'.Bpn = Cmq.(Vpq)'.Bpn = Cmq.Aqn
Z rovnosti (Vpq)'.Bpn = Aqn plyne, že odhad parametrů B‡ neovlivněných korelací mezi proměnnými X se získá jako (24)
B‡pn = Vpq.(Vpq)'.Bpn = Vpq.Aqn
Projekční matice Vpq.(Vpq)' převádí matici nestabilních regresních koeficientů Bpn v matici stabilných (ale vychýlených) regresních koeficientů B‡pn. Z rovnice 23 také plyne, že parametry regrese X na C (t.j. V) a Y na C (t.j. A) jsou součástí parametrů regrese Y na X (t.j. B). Jak úspěšně lze matici Y vysvětlit proměnnými matice X je možno testovat pomocí součtu čtverců rezidulálů (SSR) (25)
SSR = (Y - UU'Y)'.(Y - UU'Y)
respektive pomocí odpovídající reziduální standardní odchylky (RSD) (26)
9
RSD = √ [ SSR / (n (m-q))]
q
3. Nelineární regrese Jak již bylo uvedeno, v regresní rovnici je závisle proměnná y funkcí nezávisle proměnných veličin x a parametrů β: (1)
yn,1 = f( xn,p, βp,1) + εn,1
kde n je počet měření a p je počet parametrů β. Je-li alespoň jedna parciální derivace této funkce podle některého parametru b1, b2, b3,...bp nekonstantní, t.j. závislá na jiném b, pak je regresní funkce nelineární. např. pro (27) jsou parciální derivace g rovny:
y = b1 + b2. exp(b3.x) g1 = ∂y/∂b1 = 1 g2 = ∂y/∂b2 = exp(b3.x) g3 = ∂y/∂b3 = b2.x.exp(b3.x)
Protože g2 a g3 jsou nekonstantní je uvedená rovnice nelineární. Pro jinou funkci, např. y = b1 + b2.x + b3.x2 g1 = 1, g2 = x a g3 = x2 jsou parciální derivace nezávislé na parametrech b a tato rovnice je z hlediska parametrů lineární. V nelineární regresi nelze vypočítat parametry b přímo,ale jejich optimální hodnoty se hledají iterativně, t.j. postupným vylepšováním původních odhadù b01, b02, b03, ...b0p. Jedna z často používaných iterativních metod využívá linearizaci Taylorovým rozvojem; další, často používaná metoda se nazývá simplex.
LINEARIZACE TAYLOROVÝM ROZVOJEM První krok: navrhne se model y = f(x,b), např. rovnice y = b1 + b2. exp(b3.x), a původní odhad parametrů bi, t.j. b01, b02, ... b0p, získaný třeba výpočtem ze zlogaritmovaného modelu. Druhý krok: 10
z derivací gi (i = 1, 2,..p) funkce y kolem původního odhadu jednotlivých parametrů bi se sestaví Jakobiho matice J. (28)
[ ∂y/∂bi ]b=bo = n[ g1 g2...gp ] = J(n,p),
Třetí krok: funkce y se rozvine do Taylorovy řady s tím, že se expanse omezí jen na první derivaci: (29)
y = y0 + J. (b - b0)
t.j. z původních odhadů parametrů b0 a modelové rovnice (27) se vypočítá y a dosadí se jako y0 do rovnice (29); tamtéž se za y dosadí experimentální hodnoty. Čtvrtý krok: neznámý vektor parametrů b se vypočítá pomocí pseudo-inverze Jakobiho matice jako první zlepšení pùvodního odhadu parametrů b0. (30)
b - b0 = ( J' .J )-1.J'.( y - y0 )
Pátý krok: toto první zlepšení parametrù b se nyní označí jako b0, vypočte se s jeho pomocí podle rovnice (27) nové y0 a nové jednotlivé vektory gi matice J a podle rovnice (30) se získá druhý zlepšený odhad b. Tyto dva kroky se opakují tak dlouho, pokud se nový odhad liší od předcházejícího o více než předem definovanou hodnotu. Jako tato mezní hodnota se může např. zvolit norma vektoru ||b - b0|| = 0.0001. Výpočet inverse matice J'.J může skončit neúspěchem, pokud budou vektory g na sobě lineárně závislé. Proto je vhodné vložit mezi třetí a čtvrtý krok t.zv. vyšetření podmíněnosti matice J a) od Jakobiho matice se vypočítá korelační matice bez centrování dat, t.j. čtvercová matice L, jejíž prvky jsou: (gi)'.gj (30) Li,j = -------------------√[(gi)'.gi].√[(gj)'.gj] b) je-li determinant matice L < 0.001 je matice J špatně podmíněná. Důvodem pro špatnou podmíněnost matice J může být přeurčený model, t.j. v rovnici y = f(x,b) je příliš mnoho parametrù, které nelze samostatně optimalizovat. Bude-li v rovnici (31) alespoň jedno vi rùzné od nuly, je regresní model přeurčený. (31)
11
Σi=1 gi.vi = 0
p
např. pro rovnici (32)
y = b1.exp( b2 + b3.x )
jsou parciální derivace g rovny:
a rovnice
g1 = ∂y/∂b1 = exp( b2 + b3.x ) g2 = ∂y/∂b2 = b1.exp( b2 + b3.x ) g3 = ∂y/∂b3 = b1.x.exp( b2 + b3.x)
Σ gi.vi = ( v1 + v2b1 + v3b1x ).exp( b2 + b3.x ) = 0
je splněna např. když v1 = -b1, v2 = 1, v3 = 0 a model (32) je přeurčený. Přeučenost modelu je možno rovnici (32)] je možno přepsat jako (33)
někdy omezit reparametrizací. Např.
y = exp( a1 + a2.x )
kde a1 = ln b1 + b2, a2 = b3. A pak Σ gi.vi = ( v1 + v2x ).exp( a1 + a2.x ) = 0 neni splněna pro nenulová v1 a v2. Nebo při jiné reparametrizaci (34)
y = a1.exp( a2.x )
kde a1 = b1.exp( b2 ), a2 = b3. Pak rovnost Σ givi = ( v1 + v2a1x ).exp( a2.x ) = 0 platí jen jsou-li v1 i v2 rovny nule a model není přeurčený. V případě, že se reparametrizací modelu špatná podmíněnost matice J'.J nezlepší, je možno matici J'.J upravit podobně jako při tzv. hřebenové (ridge) regresi a počítat inversi matice (J'.J + αI), kde α se volí podle stupně podmíněnosti matice J'.J. Špatná podmíněnost může nastat také v případě, když průměr nezávisle proměnné hodnoty je daleko od nuly. Jednoduchá transformace - centrování těchto hodnot může podmíněnost značně zlepšit. To znamená, že rovnici typu (34) je možno počítat v alternativní formě (35) a hledaná hodnota parametru a1 se pak vypočte z a0, a2 a známé prùměrné hodnoty xp (rovnice 36). (35)
y = a0.exp( a2( x-xp ))
(36)
a1 = a0 / exp( a2.xp ).
12
OPTIMALIZACE SIMPLEXEM Simplex je prostorový útvar, který má pro p optimalizovaných paramentrů p+1 vrcholů pravidelného mnohostěnu; pro dva parametry je to rovnostranný trojúhelník, pro tři parametry je to pravidelný čtyřstěn atd. Souřadnice vrcholů simplexu představují odhady optimalizovaných paramentrů. Ze všech p+1 odhadů se vypočítají odpovídající hodnoty ycalc a zjistí se nejhorší shoda s yexp, např. pomocí standardní odchylky. Nejhorší případ se vyřadí a nový vrchol simplexu se vypočítá z odhadů parametrů získaných reflexí rovinou zbývajících p bodů. Znovu se posoudí všechny ycalc a nejhorší výsledek se opět vyřadí. Tento postup se opakuje tak dlouho,až se nový vrchol simplexu neliší od předcházejícího o více než zadanou toleranci. V případě dvou parametrů se tedy sestrojí trojúhelníkový simplex takto: jako první vrchol se umístí do počátku souřadné soustavy, t.j. bod 0,0 odhad parametrů získaný např. z odpovídající ineární regrese. další vrcholy simplexu se získají vynásobením vhodně zvoleného kroku číslem uvedeným v následující tabulce. Vrchol Počet faktorů simplexu 1 2 3 4 5 ----------------------------------------------------1 0.0 0.0 0.0 0.0 0.0 2 1.0 0.0 0.0 0.0 0.0 3 0.5 0.866 0.0 0.0 0.0 4 0.5 0.289 0.817 0.0 0.0 5 0.5 0.289 0.204 0.791 0.0 6 0.5 0.289 0.204 0.158 0.775 ----------------------------------------------------Trojúhelníkový simplex vypadá tedy např. takto:
13
Je-li nejhorší shoda vypočtených a naměřených hodnot y při parametrech B, provede se nový výpočet za podmínek D1; v případě nejhoršího výsledku při parametrech A je nový výpočet nasměrován do polohy D2. Jestliže se přesáhne reálný rozsah optimalizovaného parametru, může se vzdálenost bodu po reflexi zmenšit. Zkracování, prodlužování nebo zkracování s inverzí vzdálenosti po reflexi (t.j. modifikovaných simplex) se používají k rychlejšímu získání optimalizovaných hodnot parametrů.
14
4. Statistické testování regresních parametrů Pro posouzení kvality proložení vypočítáné funkce eperimentálními body (goodnes of fit) se používá korelační koeficient, r, a standardní odchylka regrese (standard error of estimate), se. Jinou, důležitou charakteristikou při regresi je standardní odchylka měřených dat, t.j. sy. Jednak umožňuje posoudit, jak jsou data rozptýlena od průměru (lze předpokládat, že regrese založená na větším rozsahu experimentálních dat bude významnější než ta, která využívá jen dat velmi blízkých průměru) a jednak napomáhá při rozhodování, zda některou experimentální hodnotu lze považovat za vychýlenou nebo ne. Známe-li kterékoliv dvě z těchto tří hodnot, tu třetí lze snadno dopočítat ze vztahu (37) pro korelační koeficient (samozřejmě je třeba znát i počet měření, n, a optimalizovaných parametrů, p).
(37)
r=
1−
se2 n − p = s 2y n − 1
1−
SR SY
Při hodnocení korelačního koeficientu (t.j. jeho statistické odlišnosti od nuly) se počítá jeho F test (38), který se někdy považuje za test spolehlivosti regrese jako takové; to nemusí být pravda pro malé směrnice regrese, t.j. při malé hodnotě sy. (38)
r 2 n − p SY − SR n − p F= = ≥ F(0.05 , p −1,n − p ) SR p −1 1− r 2 p −1
Ke zjištění statistické významnosti jednotlivých paramentrů b ve více-rozměrné lineární nebo nelineární regresi je třeba znát jejich standardní odchylky sb. Ty se vypočítají ze standardní odchylky regrese, se, a z odmocniny odpovídajících diagonálních elementů matice V. (39)
V = (Z' Z)-1
kde Z = [ g1 g2 g3 ....gp ]
a g jsou derivace funkce y podle jednotlivých paramentrů b. Pro lineární regresi se matice Z rovná matici X. (40)
sb = se √ ( diag (V))
Statistická významnost se pak zjistí Studentovým t-testem tak, že se vypočtená hodnota t (41) porovná s hodnotou kritickou, tabelovanou (viz kapitola 2) pro hladinu spolehlivosti (většinou 0.05, t.j. 95%) a počet stupňů volnosti (n-p): je-li splněna nerovnost (41), pak testovaný rozdíl (b - β) je statisticky významný.
15
t=
(41)
b−β ≥ t (0.05, n − p ) sb
Je-li β = 0, pak se testuje statistická významnost odpovídající proměnné (resp. odchylka příslušného paramentru od nuly), v případě, že β = 1, lze takto testovat významnost odchylky odpovídající směrnice od jedné. Statistická významnost přidání další proměnné do regresní závislosti se dá také testovat pomocí parciálního F-testu (42).
F=
(42)
SR p − SR p +1 SR p +1
(n − p ) ≥ F(0.05,1,n− p )
kde SRp je součet čtverců reziduálů z regrese s p paramentry a SRp+1 je totéž pro regresi s p+1 parametry. Vypočítaná hodnota F se srovnává s kritickou hodnotou tabelovanou pro zvolenou hladinu statistické významnosti a odpovídající počet stupňů volnosti: je-li splněna nerovnost (42), pak přidání další proměnné je statisticky významné. Jiným testem pro ocenění kvality v případě lineárního regresního modelu je analýza reziduálů, ε = yexp-ycalc. K zvýraznění reziduálů, které mohou patřit podezřelému odlehlému měření, se reziduály převádějí na t.zv. 'jackknife' reziduály (43). Tento přepočet simuluje výpočet reziduálů s vynecháním konkretního měření, t.j. jisté měření yexp se vynechá, regresní parametry se vypočítají ze zbylých měření y, z těchto parametrů a z odpovídajících hodnot vysvětlujících proměnných, x, se vypočítá ycalc a porovná s experimentální hodnotou: εjk = yexp – ycalc. (43)
ε jk = ε n
n − p −1 n − p − ε 2n
kde
εn =
se
ε 1 − hi
a hi je odpovídající diagonální element transformační matice H=X(X'X)-1X'. Jeli εjk > 0.2 (n-p), pak se yi považuje za měření podezřelé z odlehlosti a regrese by se měla zopakovat bez něho. Zjistí-li se, že nové paramentry se neliší od dřívějších o více než jednu standardní odchylku sb, nebo že se při nové regresi zvýšila standardní odchylka experimentálních dat od průměru, sy (t.j. v důsledku snížení počtu stupňů volnosti ve jmenovateli), pak můžeme podezření z odlehlosti opustit.
16
ČÁST 2. INTEGROVANÉ RYCHLOSTNÍ ROVNICE 1. Unimolekulární vratná reakce
kA
A
B
kB Rychlost této reakce υ je možno vyjádřit např. jako úbytek reaktantu A s časem:
υ =
dA t = −k A A t + k B B t dt
Z výpočetních důvodů je vhodné, nahradit dvě proměnné veličiny, At a Bt, t.j. koncentrace reaktantu A a reaktantu B v čase t jednou pomocnou proměnnou x, která představuje zreagované množství v čase t a nazývá se rozsah reakce (angl. conversion variable) x = A0 - At = Bt - B0 potom
−
dA t dx = = −(k A + k B )x + k A A 0 − k B B 0 = − kx + m dt dt t dx m − kx dt = ⇒ ln = −kt ∫ ∫ m kx − m 0 0
x
nechť
m =F k
potom po odlogaritmování
F−x = exp( −kt ) F
a po dosazení za x dostaneme rovnice pro pokles koncentrace A a nárust koncentrace B s časem t.
A t = A 0 − F + F exp( −kt )
(1)
B t = B 0 + F − F exp( − kt )
Rovnovážná konstanta reakce je (čas dosažení rovnováhy je označen ∞ ): (2)
17
K=
B∞ k A B0 + F = = A∞ k B A0 − F
protože
F=
KA 0 − B 0 K +1
Varianta 1: Pro případ, že K = 1, což nastane např. při racemizaci:
F= a rovnice (1) se změní takto:
A 0 − B0 2
A 0 + B0 A 0 − B0 exp( − kt ) + 2 2 A + B 0 A 0 − B0 Bt = 0 − exp( − kt ) 2 2
At = (3)
Jestliže se racemizace sleduje polarimetricky, pak naměřená optická otáčivost je úměrná rozdílu koncentrací složek A a B, t.j. At – Bt: (4)
A t − B t = (A 0 − B 0 ) exp( −kt ) αt
= α 0 exp( −kt )
Pro nelineární regresi je vhodné doplnit rovnici (4) o konstantní parametr a0, který představuje systematickou chybu měření (nepřesné vynulování přístroje) a o ε, což je reziduál neboli náhodná chyba měření. (5)
α t = a 0 + α 0 exp( − kt ) + ε
Varianta 2: Pro případ, kdy počáteční koncentrace reagentu B je nulová, B0 = 0, t.j. vychází se z čistého reagentu A, není třeba zavádět do výpočtu pomocnou proměnnou x.
dA t = −k A A t + k B B t dt Pak se za Bt dosadí z rovnice
At + Bt = A0
⇒
Bt = A0 - At
dA t = −(k A + k B )A t + k B A 0 = − kA t + m 0 dt
18
pro
m0 = F0 k
F0 − A t = exp( − kt ) F0 − A 0
po odlogaritmování
Po dosazení z rovnice At=A0-Bt
A t = F0 − ( F0 − A 0 ) exp( − kt )
(6)
B t = A 0 − F0 + ( F0 − A 0 ) exp( − kt )
protože F0 =
k BA 0 A A − F0 je rovnovážná konstanta K = 0 = 0 kA + kB K +1 F0
Je-li v tomto případě K = 1 pak
At =
F=
A0 A0 + exp( −kt ) 2 2
A0 2 Bt =
A0 A0 − exp( −kt ) 2 2
a podobně jako v předcházející variantě
A t − B t = A 0 exp( −kt ) resp. α t = α 0 exp( −kt ) Varianta 3: Jiný případ nastane, když rychlostní konstanta zpětné reakce, t.j. kB, je nulová, nebo vzhledem ke kA zanedbatelně malá. Potom reakční rychlost υ je dána
υ =− At
dA t dx = = k AAt = k A ( A0 − x ) dt dt
t dA t = −k ∫ dt ∫ A t A0 0
⇒
ln A t − ln A 0 = −kt
t.j. (7)
A t = A 0 exp( − kt ) resp. x = A 0 ( 1 − exp( −kt ))
kde x je přírustek koncentrace produktu B (t.j. Bt) a rovnice (7) jsou shodné s rovnicemi (6) pro F0 = 0.
19
2. Bimolekulární reakce
P
A + B
Při této reakci se mohou výchozí koncentrace A0 a B0 lišit, nebo mohou být stejné; podle toho dostáváme odpovídající kinetické rovnice. Varianta 1: Pro A0 > B0
a
At = A0 – x,
Bt = B0 = x
dx = kA t B t = k (A 0 − x )(B 0 − x ) dt x x M t dx M −1 =∫ −∫ = k ∫ dt ; kde M = ∫ A 0 − B0 0 (A 0 − x )(B 0 − x ) 0 A 0 − x 0 B 0 − x 0
x
po integraci:
−
A0 B0 1 1 ln + ln = kt A 0 − B0 A 0 − x A 0 − B0 B0 − x
B A −x = (A 0 − B 0 ) kt ln 0 . 0 A 0 B0 − x
V
nechť
V=
B0 A0
A0 − x = exp ((A 0 − B 0 ) kt ) = E B0 − x
pak
V A 0 − V x = E B 0 − E x ⇒ x (V − E ) = B 0 (1 − E )
(8)
20
A 0 − B0 1− V x = B 0 1 + = B 1 + 0 V − E B − A exp ( ( A − B ) kt ) 0 0 0 0
Varianta 2: Pro A0 = B0
2A −
P
1 dA t = kA t 2 2 dt
⇒ −
dA t At
2
= 2kdt
po integraci
1 1 − = 2kt ⇒ A t A0
(9)
At =
1 1
A0
+ 2kt
A0 = 1 + A 0 2kt At
⇒ At = a0 +
1 a1 + a 2 t
Je vhodné rozšířit rovnici o člen a0 a po výpočtu testovat jeho statistickou významnost. 3. Enzymově katalyzovaná reakce
k1 E+A
k2 M
P+E
k-1 Rychlost reakce, t.j. rychlost vzniku produktu P závisí na koncentraci meziproduktu M a rychlostní konstantě k2.
υ=
dPt = k 2M t dt
Koncentrace meziproduktu v čase t je dána
dM t = k1A t E t − k −1M t − k 2 M t = 0 dt
21
Předpokládá se, že v ustáleném stavu je koncentrace meziproduktu malá a prakticky s časem neměnná; proto je možno položit dMt/dt = 0. Celková (t.j. původní) koncentrace enzymu je v čase t rozdělena mezi enzym volný Et a enzym v meziproduktu Mt.
E0 = Et + Mt Poslední dvě rovnice se vyjádří maticově takto:
1 k A 1 t
E t E 0 = − (k −1 + k 2 ) M t 0 1
Koncentraci meziproduktu Mt potřebnou pro výpočet reakční rychlosti vypočítáme pomocí inverze čtvercové matice: Její determinant
D = −(k −1 + k 2 ) − k1A t E t 1 − (k −1 + k 2 ) − 1 E 0 M = D − k A 1 0 t 1 t
Mt =
− k 1A t E 0 D
⇒ υ=
k 1k 2 A t E 0 (k −1 + k 2 ) + k 1A t
Po dosazení KM = (k-1+k2)/k1 a VM = k2 E0 dostaneme rovnici MichaelisovuMentenové:
υ=
VM A t KM + At
kde VM je maximální rychlost reakce a KM je Michaelisova konstanta (t.j. koncentrace substrátu, při níž dosahuje reakční rychlost polovinu rychlosti maximální). Pro nelineární regresi je vhodné tuto rovnici upravit do tvaru:
y = a0 +
22
a 1x a2 + x
ČÁST 3. VÝPOČTOVÉ PROGRAMY 1. Nelin: exponenciální regrese s Taylorovým rozvojem
% nelin.m clc disp('program nelin.m verze 27.6.1999') if exist('Bh'), clear Bh; end if exist('x1'), clear x1; end if exist('y1'), clear y1; end disp('pro vypocet parametru nelinearni regrese') disp('1.model: y = b1*exp(b2*x) ') disp('2.model: y = b1 + b2*exp(b3*x) ') wm=input('vyber model wm = '); x1=input('zaved sloupcovy vektor x1 = '); y1=input('zaved vektor amplitud y1 = '); un=ones(size(y1)); md=length(x1); X=[un x1]; if y1(1)
23
b00=b0; if y1(1)
Z=[g1 g2 g3]; y0=b0(1)+b0(2)*g2; end RSS=(y1-y0)'*(y1-y0); disp([' linearni odhad a1 = ' num2str(a1)]); disp([' linearni odhad a2 = ' num2str(a2) '; RSS = ' num2str(RSSa)]); if wm==1 disp(['nelinearni odhad b1 = ' num2str(b0(1))]); disp(['nelinearni odhad b2 = ' num2str(b0(2)) '; RSS = ' num2str(RSS)]); elseif wm==2 disp(['nelinearni odhad b1 = ' num2str(b0(1))]); disp(['nelinearni odhad b2 = ' num2str(b0(2))]); disp(['nelinearni odhad b3 = ' num2str(b0(3)) '; RSS = ' num2str(RSS)]); end % vyšetření podmíněnosti poslední matice LL=korr(Z); del=det(LL); if del<pod disp(['posledni model je spatne podmineny: det(LL) = ' num2str(del)]) else disp(['posledni model neni spatne podmineny: det(LL) = ' num2str(del)]) end se2=RSS/ny; se=sqrt(se2); disp(['se = ' num2str(se)]) FII=iZ*iZ'; sb=sqrt(se2*diag(FII)); bsb=abs(b0./sb); disp(' b0 sb bsb'); format long disp([b0 sb bsb]); format short st=student(ny); disp(['je-li bsb vetsi nez ' num2str(st) ', zamita se nulova hypoteza']); % end function y = korr(x) [m,n] = size(x); if (m == 1) m = n; x = x.'; end c = x' * x / m; 26
d = diag(c); y = c./sqrt(d*d');
function y = student(x) z4=4.549/x; z3=(0.621+z4)/x; z2=(3.226+z3)/x; z1=(2.350+z2)/x; y=1.96+z1;
27
Protokol z výpočtu programem Nelin.m
Data uložená pod názvem "bradam" pocházejí z polarimetrického sledování rotační enantiomerace 3-brom-2,4,6-trimethylfenyladamantyl ketonu (dosud nepublikováno). t/s α/mdeg t/s α/mdeg t/s α/mdeg ----------------------------------------------------------0 464 1740 267 4260 125 120 444 1980 247 4620 113 240 425 2220 231 4980 102 360 412 2460 212 5340 93 480 395 2700 200 5700 84 660 375 2940 186 6060 77 900 357 3180 172 6420 68 1260 311 3540 155 6780 62 1500 288 3900 139
----------------------------------------------------Program nelin.m verze z 27.6.1999 pro vypocet parametru nelinearni regrese 1.model: y = b1*exp(b2*x) 2.model: y = b1 + b2*exp(b3*x) vyber model wm = 2 zaved sloupcovy vektor x = t zaved vektor amplitud y = alfa puvodni model neni spatne podmineny: det(LL) = 0.011556 linearni odhad a1 = 452.5912 linearni odhad a2 = -0.00029732; RSS = 488.1518 nelinearni odhad b1 = 13.7205 nelinearni odhad b2 = 448.3981 nelinearni odhad b3 = -0.00032631; RSS = 116.1667 posledni model neni spatne podmineny: det(LL) = 0.016037 se = 2.2474 b0 sb bsb ----------------------------------13.7205 2.6391 5.1989 448.3981 2.3461 191.1255 -0.0003 0.0000 72.4645
je-li bsb vetsi nez 2.0683, zamita se nulova hypoteza program nelin.m pro vypocet parametru nelinearni regrese 1.model: y = b1*exp(b2*x) 2.model: y = b1 + b2*exp(b3*x) vyber model wm = 1 28
verze z 27.6.1999
zaved sloupcovy vektor x = t zaved vektor amplitud y = alfa puvodni model neni spatne podmineny: det(LL) = 0.56568 linearni odhad a1 = 452.5912 linearni odhad a2 = -0.00029732; RSS = 488.1518 nelinearni odhad b1 = 459.4156 nelinearni odhad b2 = -0.00030541; RSS = 238.7313 posledni model neni spatne podmineny: det(LL) = 0.56844 se = 3.1539 b0 sb bsb ---------------------------------------459.4156 1.4270 321.9511 -000.0003 0.0000 161.0068 je-li bsb vetsi nez 2.0636, zamita se nulova hypoteza ------------------------------------------------------------------------------------------------lineární odhad
----------------------------------------------------------------------------------------------------------
nelineární odhad 29
2. Simpl: exponenciální regrese s modifikovaným simplexem
% simpl.m clc disp('program simpl.m verse 24.3.1999') if exist('b0'),clear b0; end global x1 y1 md ny wm b0 RSS smo disp('simplex (NELDER.M) na exponencialni rovnici ') disp('1.model: y=exp(b1+b2*x)') disp('2.model: y=b1+exp(b2+b3*x)') wm=input('vyber model wm = '); x1=input('zaved vektor x = '); y1=input('zaved vektor y = '); [md nd]=size(x1); un=ones(md,1); X=[un x1]; a=inv(X'*X)*X'*log(y1); if wm==1 b00=a; b0=b00; y0=exp(b0(1)+b0(2)*x1); elseif wm==2 b00=[0;a]; b0=b00; y0=b0(1)+exp(b0(2)+b0(3)*x1); end plot(x1,y1,'ko',x1,y0,'k-') wk=b0; wk=nelder('simpla',wk,0.001); disp(' '); disp(b00') disp(b0') disp(['RSS = ' num2str(RSS)]) disp(['smo = ' num2str(smo)]) function fn=simpla(wk) global x1 y1 md ny wm b0 RSS smo b0=wk; if wm==1 yy=exp(b0(1)+b0(2)*x1); ny=md-2; elseif wm==2 yy=b0(1)+exp(b0(2)+b0(3)*x1); 30
ny=md-3; end RSS=(y1-yy)'*(y1-yy); smo=sqrt(RSS/ny); plot(x1,y1,'ko',x1,yy,'k-'),title(['RSS = ' num2str(RSS)]) pause(0.1) fn=smo; function x = nelder('simpla',x,tol) % Nelderův-Meadův simplexový algoritmus pro minimalizaci % nelineární funkce několika proměnných. % 'simpla' je jméno funkce, která má být minimalizována, % x je vektor, který se optimalizuje na nový vektor % minimalizující funkci 'simpla'. % tol je tolerance na zastavení optimalizace (normálně 0.001) [n,m] = size(x); if m > n x = x'; n = m; end if nargin < 3, tol = 1.e-3; end v = 0.9*x; f = feval(F,v); for j = 1:n y = x; if y(j) ~= 0 y(j) = 1.1*y(j); else y(j) = 0.1; end v = [v y]; f = [f feval(F,y)]; end [f,j] = sort(f); v = v(:,j); while 1 test = 0; for j = 2:n+1, test = max(test,norm(v(:,j)-v(:,1),1)); end if test <= tol, break, end [v,f,how] = neldstep(F,v,f); end x = v(:,1); 31
function [v,f,how] = neldstep(F,v,f) % funkce používaná simplexovým algoritmem NELDER; % v je simplex (je to n,n+1 matice); hodnoty funkce jsou: % f(j) = fun(V(:,j)), j = 1:n+1, with f(1) <= ... <= f(n+1) % F = 'simpla'; výsledkem je nový simplex v % how popisuje krok, který se právě koná alpha = 1; beta = 1/2; gamma = 2; [n,np1] = size(v); vbar = mean(v(:,1:n)')'; vr = (1 + alpha)*vbar - alpha*v(:,n+1); fr = feval(F,vr); vk = vr; fk = fr; how = 'reflect '; if fr < f(n) if fr < f(1) ve = gamma*vr + (1-gamma)*vbar; fe = feval(F,ve); if fe < f(1) vk = ve; fk = fe; how = 'expand '; end end else vt = v(:,n+1); ft = f(n+1); if fr < ft vt = vr; ft = fr; end vc = beta*vt + (1-beta)*vbar; fc = feval(F,vc); if fc < f(n) vk = vc; fk = fc; how = 'contract'; else for j = 2:n v(:,j) = (v(:,1) + v(:,j))/2; f(:,j) = feval(F,v(:,j)); end vk = (v(:,1) + v(:,n+1))/2; fk = feval(F,vk); how = 'shrink '; end end v(:,n+1) = vk; f(:,n+1) = fk; [f,j] = sort(f); v = v(:,j);
32
% konfi.m verze 4.7.1999 clc disp('program konfi.m-pocita interval spolehlivosti, t.j. ± c') disp('pro parametry nelinearni rovnice') disp('1.model: y = b1 + b2*exp(b3*x)') disp('2.model: y = b1 + exp(b2+b3*x)') wm=input('vyber model wm = '); x0=input('zaved vektor promennych x0 = '); b0=input('zaved vektor parametru [b0] = '); smo=input('zaved res.smer.odchylku smo = '); x0=x0(:); md=length(x0); g1=ones(size(x0)); if wm==1 g2=exp(b0(3)*x0); g3=b0(2)*(x0.*g2); elseif wm==2 g2=exp(b0(2)+b0(3)*x0); g3=x0.*exp(b0(2)+b0(3)*x0); end Z=[g1 g2 g3]; IZ=inv(Z'*Z); VV=(diag(IZ)*(smo^2)); sb=VV.^(1/2); tkrit=student(md-3); c=tkrit*sb; disp(' ') disp(' b sb |b|/sb ± c ') disp('------------------------------------------') disp([b0 sb abs(b0)./sb c]); disp(['tkrit = ' num2str(tkrit)]); % end
33
Protokol z výpočtu programem Simpl.m Data pro závislost rychlostní konstanty na teplotě: J.Sandström, Dynamic NMR Spectroscopy, Academic Press, 1982, str. 155. T/K iT/10-3K-1 ka/s-1 ----------------------------------------------------------7.3 301.55 3.3162 10.5 306.15 3.2664 19.5 311.55 3.2098 28.5 315.95 3.1651 48.0 321.15 3.1138 71.0 327.95 3.0492 125.0 333.95 2.9945 145.0 336.35 2.9731 230.0 341.95 2.9244 260.0 344.75 2.9007 355.0 347.95 2.8740 430.0 351.75 2.8429 525.0 354.35 2.8221 600.0 355.55 2.8125 950.0 362.15 2.7613 1300.0 367.35 2.7222 1700.0 372.75 2.6828 2300.0 377.75 2.6473 -----------------------------------------------------------
Program simpl.m simplex (NELDER.M) na exponencialni rovnici 1.model: y=exp(b1+b2*x) 2.model: y=b1+exp(b2+b3*x) vyber model wm = 1 zaved vektor x = iT zaved vektor y = ka linearni odhad b1 a b2: 30.6 -8631.7 nelinearni odhad b1 a b2: 29.6 -8241.8 RSS = 8497.0592 smo = 23.0449
verse 2.8.1999
Program simpl.m simplex (NELDER.M) na exponencialni rovnici 1.model: y=exp(b1+b2*x) 2.model: y=b1+exp(b2+b3*x) vyber model wm = 2 zaved vektor x = iT zaved vektor y = ka
verse 2.8.1999
34
linearni odhad b2 a b3: nelinearni odhad b1, b2 a b3: -12.9 RSS = 7670.5379 smo = 22.6135
30.6 -8631.7 29.2 -8100.9
program konfi.m - pocita interval spolehlivosti, t.j. +-c pro parametry nelinearni rovnice 1.model: y = b1 + b2*exp(b3*x) 2.model: y = b1 + exp(b2+b3*x) vyber model wm = 2 zaved vektor promennych x0 = iT zaved vektor parametru [b0] = b0 zaved res.smer.odchylku smo = smo b sb |b|/sb +-c -----------------------------------------1.0e+003 * -0.0129 0.0103 0.0013 0.0220 0.0292 0.0004 0.0725 0.0009 -8.1009 0.1517 0.0534 0.3232
tkrit = 2.1313 Hodnocení: parametr b1 není statisticky významný -----------------------------------------------------------------------------------------Model 1, bez konstantního členu
35
3. Rad2a: nelineární regrese s Taylorovým rozvojem
% rad2a.m clc % verse 27.9.1999 if exist('Bh'), clear Bh; end if exist('x'), clear x;end if exist('y'), clear y;end disp('program rad2a.m pro reseni rovnice (1/A)=(1/Ao)+2kt') disp('jako nelinearni rovnice A=a1+(a2+a3*t)^(-1)'); t =input('zaved casy mereni t = '); A =input('zaved odezvy A v case t A = '); x=t; y=A; md=length(x); un=ones(size(y)); ya=un./y; X=[un x]; a=inv(X'*X)*X'*ya; disp(['linearni odhad A0 = ' num2str(1/a(1))]); disp(['linearni odhad ka = ' num2str(a(2)/2)]); b0=[0; a(1); a(2)]; b01=b0; g1=un; gu=b0(2)+b0(3)*x; g2=(-1)*gu.^(-2); g3=g2.*t; Z=[g1 g2 g3]; y0=b0(1)+gu.^(-1); RSSa=(y-y0)'*(y-y0); pod=0.001; LL=korr(Z); del=det(LL); if del<pod disp(['puvodni model je spatne podmineny: det(LL) = ' num2str(del)]) else disp(['puvodni model neni spatne podmineny: det(LL) = ' num2str(del)]) end Bn=100; h=0; RSS=RSSa; while Bn>0.000001 h=h+1; disp(h) 36
if h>200; break; end figure(1) plot(x,y,'ko',x,y0),grid,pure,title(['RSS = ' num2str(RSS)]); pause(1e-6) g1=un; gu=b0(2)+b0(3)*x; g2=(-1)*gu.^(-2); g3=g2.*t; Z=[g1 g2 g3]; y0=b0(1)+gu.^(-1); RSS=(y-y0)'*(y-y0); [u s v]=svd(Z,0); iZ=v*inv(s)*u'; B=iZ*(y-y0); Bn=norm(B); Bh(h,:)=[B' b0']; b0=b0+B; end g1=un; gu=b0(2)+b0(3)*x; g2=(-1)*gu.^(-2); g3=g2.*t; Z=[g1 g2 g3]; y0=b0(1)+gu.^(-1); RSS=(y-y0)'*(y-y0); b2=b0(2); b3=b0(3); disp(['nelinearni odhad A0 = ' num2str(1/b2)]); disp(['nelinearni odhad ka = ' num2str(b3/2) '; RSS = ' num2str(RSS)]); LL=korr(Z); del=det(LL); if del<pod disp(['posledni model je spatne podmineny: det(LL) = ' num2str(del)]) else disp(['posledni model neni spatne podmineny: det(LL) = ' num2str(del)]) end se2=RSS/(md-3); se=sqrt(se2); disp(['se = ' num2str(se)]) FII=iZ*iZ'; sb=sqrt(se2*diag(FII)); bsb=abs(b0./sb); disp(' b0 sb bsb'); 37
format long disp([b0 sb bsb]); format short st=student(md-3); disp(['je-li bsb vetsi nez ' num2str(st) ', zamita se nulova hypoteza']); % end
38
Protokol z výpočtu programem Rad2a.m
Data pro závislost vodivosti vodného roztoku NaOH na postupu hydrolýzy octanu ethylnatého: E.A.Guggenheim, J.E.Prue, Physicochemical Calculations, Interscience Publishers Inc., New York, 1955, str.432 (ruský překlad 1958). t/s κ(rel) -------------------------------------------------0 1.560 5 1.315 7 1.247 9 1.193 11 1.146 13 1.107 15 1.064 18 1.020 20 0.994 25 0.945 27 0.923 --------------------------------------------------
Program rad2a.m pro reseni rovnice (1/A)=(1/Ao)+2kt jako nelinearni rovnice A=a1+(a2+a3*t)^(-1) zaved casy mereni t = t zaved odezvy A v case t A = kp linearni odhad A0 = 1.4625 linearni odhad ka = 0.0078805 puvodni model je spatne podmineny: det(LL) = 0.00018164 nelinearni odhad A0 = 1.001 nelinearni odhad ka = 0.032245; RSS = 5.0405e-005 posledni model neni spatne podmineny: det(LL) = 0.0046891 se = 0.0025101 b0 sb bsb ------------------------------------------0.5589 0.0085 65.729 0.9990 0.0079 127.112 0.0645 0.0017 36.944 je-li bsb vetsi nez 2.3065, zamita se nulova hypoteza
39
4. Michme: nelineární regrese Taylorovým rozvojem
% michme.m verze 6.10.1999 clc if exist('x1'), clear x1; end if exist('y1'), clear y1; end if exist('Bh'), clear Bh; end disp('program michme.m') disp('nelinearni vypocet parametru rovnice Michaelisovy-Mentenove') disp('1.model: y = b1*x/(b2+x) '); disp('2.model: y = b1 + b2*x/(b3+x) '); wm=input('vyber model wm = '); x1=input('zaved sloupcovy vektor x1 = '); y1=input('zaved vektor amplitud y1 = '); [md nd]=size(x1); un=ones(md,1); iv=un./y1; is=un./x1; X=[un is]; a=inv(X'*X)*X'*iv; yy=X*a; rra=corr([iv yy]); ra=rra(1,2); RSSa=(iv-yy)'*(iv-yy); RSS=RSSa; figure(1) plot(is,iv,'k*',is,yy) title(['RSSa = ' num2str(RSSa)]); xlabel('linearni prolozeni') pause figure(2) plot(x,y,'*',x,un./yy) title(['RSSa = ' num2str(RSSa)]); xlabel('nelinearni prolozeni') pause a1=a(1); a2=a(2); if wm==1 b0=[1/a1 a2/a1]'; b00=b0; g1=x1./(b0(2)+x1); g2=-(b0(1)*x1)./((b0(2)+x1).^2); Z=[g1 g2]; 40
y0=b0(1)*g1; elseif wm==2 b0=[0 1/a1 a2/a1]'; b00=b0; g1=un; g2=x1./(b0(3)+x1); g3=-(b0(2)*x1)./((b0(3)+x1).^2); Z=[g1 g2 g3]; y0=b0(1)+b0(2)*g2; end % vysetreni podminenosti puvodni matice pod=0.001; LL=korr(Z); del=det(LL); if del<pod disp(['puvodni model je spatne podmineny: det(LL) = ' num2str(del)]) else disp(['puvodni model neni spatne podmineny: det(LL) = ' num2str(del)]) end Bn=100; h=0; while Bn>0.000001 h=h+1; figure(3) plot(x1,y1,'k*',x1,y0,x1,un./yy,'o') title(['RSS = ' num2str(RSS)]); if wm==1 g1=x1./(b0(2)+x1); g2=-(b0(1)*x1)./((b0(2)+x1).^2); Z=[g1 g2]; y0=b0(1)*g1; elseif wm==2 g1=un; g2=x1./(b0(3)+x1); g3=-(b0(2)*x1)./((b0(3)+x1).^2); Z=[g1 g2 g3]; y0=b0(1)+b0(2)*g2; end RSS=(y1-y0)'*(y1-y0); B=inv(Z'*Z)*Z'*(y1-y0); Bn=norm(B); Bh(h,:)=[B' b0']; b0=b0+B; 41
end if wm==1 g1=x1./(b0(2)+x1); g2=-(b0(1)*x1)./((b0(2)+x1).^2); Z=[g1 g2]; y0=b0(1)*g1; elseif wm==2 g1=un; g2=x1./(b0(3)+x1); g3=-(b0(2)*x1)./((b0(3)+x1).^2); Z=[g1 g2 g3]; y0=b0(1)+b0(2)*g2; end RSS=(y1-y0)'*(y1-y0); rrb=corr([y1 y0]); rb=rrb(1,2); disp(['linearni odhad V (1/a1) = ' num2str(1/a1) '; ra = ' num2str(ra)]); disp(['linearni odhad K (a2/a1) = ' num2str(a2/a1) '; RSS ='num2str(RSSa)]); if wm==1 disp(['nelinearni odhad V (b1) = ' num2str(b0(1)) '; rb = ' num2str(rb)]); disp(['nelinearni odhad K (b2) = ' num2str(b0(2)) '; RSS = ' num2str(RSS)]); ny=md-2; elseif wm==2 disp(['nelinearni odhad (b1) = ' num2str(b0(1)) '; rb = ' num2str(rb)]); disp(['nelinearni odhad V (b2) = ' num2str(b0(2)) '; rb = ' num2str(rb)]); disp(['nelinearni odhad K (b3) = ' num2str(b0(3)) '; RSS = ' num2str(RSS)]); ny=md-3; end % vysetreni podminenosti posledni matice LL=korr(Z); del=det(LL); if del<pod disp(['posledni model je spatne podmineny: det(LL) = ' num2str(del)]) else disp(['posledni model neni spatne podmineny: det(LL) = ' num2str(del)]) end se2=RSS/ny; se=sqrt(se2); disp(['se = ' num2str(se)]) izz=inv(Z'*Z); sb=sqrt(se2*diag(izz)); bsb=abs(b0./sb); disp(' b0 sb bsb'); 42
format long disp([b0 sb bsb]); format short st=student(ny); disp(['je-li bsb vetsi nez ' num2str(st) ', zamita se nulova hypoteza']);
43
Protokol z výpočtu programem Michme.m
Data nasimulovaná a zatížená náhodnou chybou. program michme.m nelinearni vypocet parametru rovnice Michaelisovy-Mentenove 1.model: y = b1*x/(b2+x) 2.model: y = b1 + b2*x/(b3+x) vyber model wm = 2 zaved sloupcovy vektor x1 = xi zaved vektor amplitud y1 = yr xi yr xi yr -------------------------------------------------0.2 0.1492 2.2 0.3334 0.4 0.2168 2.4 0.3408 0.6 0.2532 2.6 0.3537 0.8 0.2795 2.8 0.3471 1.0 0.2915 3.0 0.3556 1.2 0.3162 4.0 0.3663 1.4 0.3177 5.0 0.3665 1.6 0.3340 6.0 0.3638 1.8 0.3326 7.0 0.3719 2.0 0.3407 8.0 0.3833 -------------------------------------------------puvodni model je spatne podmineny: det(LL) = 0.00098105 linearni odhad V (1/a1) = 0.39231; ra = 0.99903 linearni odhad K (a2/a1) = 0.32542; RSS = 0.032584 nelinearni odhad (b1) = -0.007914; rb = 0.99672 nelinearni odhad V (b2) = 0.3989; rb = 0.99672 nelinearni odhad K (b3) = 0.30952; RSS = 0.00042092 posledni model je spatne podmineny: det(LL) = 0.00089643 se = 0.0049759 b0 sb bsb -----------------------------0.0079 0.0220 0.3592 0.3989 0.0203 19.6366 0.3095 0.0347 8.9143 ---------------------------------------------------je-li bsb vetsi nez 2.1096, zamita se nulova hypoteza
44
program michme.m nelinearni vypocet parametru rovnice Michaelisovy-Mentenove 1.model: y = b1*x/(b2+x) 2.model: y = b1 + b2*x/(b3+x) vyber model wm = 1 zaved sloupcovy vektor x1 = xi zaved vektor amplitud y1 = yr puvodni model neni spatne podmineny: det(LL) = 0.34423 linearni odhad V (1/a1) = 0.39231; ra = 0.99903 linearni odhad K (a2/a1) = 0.32542; RSS = 0.032584 nelinearni odhad V (b1) = 0.39166; rb = 0.99669 nelinearni odhad K (b2) = 0.32181; RSS = 0.00042429 posledni model neni spatne podmineny: det(LL) = 0.34585 se = 0.0048551 b0 sb bsb ---------------------------------0.3917 0.0022 176.3955 0.3218 0.0106 30.3028 ------------------------------------------------------------je-li bsb vetsi nez 2.1007, zamita se nulova hypoteza
45
