Sztochasztikus rezonancia a kémiai kinetikában

Izsák Ferenc V. éves matematikus hallgató^{1}

témavezető: Tóth János Budapest, 2000. április 30.

1A szakdolgozat a T025472. számú OTKA pályázat részleges támogatásával készült.

Kivonat

A reakciókinetika tárgya a kémiai reakciók jellemző mennyiségeinek vizsgálata. Ezek leírására többféle matematikai modellt is kidolgoztak, melyek közt az egyik fő különbség az, hogy milyen mennyiségeket vesznek a reakció során állandónak. Mi itt a legegyszerűbb esetre szorítkozunk, amennyiben egyedül az anyagmennyiség, azaz rögzített térfogat esetén a koncentráció változását feltételezzük. A másik fő különbség a lehetséges modellek közt, hogy determinisztikus vagy sztochasztikus módon írják-e le a reakciót. A dolgozat első részében az utóbbi két típushoz tartozó modell összehasonlítására irányítjuk figyelmünket; megmutatjuk, hogy egy konkrét tulajdonságra, a felrobbanásra nézve a két modell lényegében ekvivalens. Ha azonban nem csak egy esemény bekövetkezésére, hanem a reakció viselkedését részletesebben leíró mennyiségekre vagyunk kíváncsiak, akkor kiderül, hogy a sztochasztikus modell erre (legalábbis elvileg) sokkal inkább alkalmas. Ennek konkrét megfogalmazásával zárjuk az első részt. További modellek tárgyalása, azok közti aszimptotikus összefüggések leírása található a később is idézett [1] cikkben. Kémiai reakciók során vagy különböző anyagcsere-, ill. idegrendszeri folyamatoknál a zaj formájában jelentkező „véletlen hatások” is fontos szerephez jutnak. Konstans mennyiségek mérését vagy reakciók bizonyos paramétereinek meghatározását akár pontosabbá is tehetik. Ezt a sztochasztikus rezonanciának nevezett jelenséget példákon keresztül részletesebben megvilágítjuk és megvizsgáljuk a második részben. Először a konstans érték mérését modellezzük, és tisztán matematikai eszközökkel reprodukáljuk a [6] cikkben ismertetett laboratóriumi mérés eredményeit. Az utolsó fejezetben pedig azt vizsgáljuk meg, hogyan javítható zaj hozzáadásával egy determinisztikusan leírt folyamat kezdeti értékére adott becslés. A megfelelő részekhez kapcsolódó három program és három ábra a mellékletben talalálható. A dolgozatban külön hivatkozás nélkül használjuk a reakciókinetika, a valószínűségszámítás és a matematikai statisztika néhány alapfogalmát.

1. fejezet

1.1. A differenciálegyenlet által leírt kémiai rendszer felrobbanásának vizsgálata

A reakciókinetika sztochasztikus és determinisztikus modellezésének összehasonlítását a felrobbanás problémájának vizsgálatával érzékeltetjük. A rendszer felrobbanása azt jelenti, hogy determinisztikus módon, tömeghatás típusú kinetikával modellezve a kémiai reakciót, a neki megfeleltetett differenciálegyenlet megoldása az adott kezdeti feltétellel nem értelmezhető egy egész félegyenesen, mert a megoldás valamelyik komponensének határértéke véges időpontban +∞ (Úgy is fogalmazhatunk, hogy a megoldás értéke véges időn belül végtelenné válik.) Lényegében egyváltozós differenciálegyenletekre fogunk szorítkozni.

1.1.1. A determinisztikus eset

Vizsgáljuk meg, hogy egyváltozós esetben milyen feltételek mellett van a differenciálegyenletnek olyan megoldása, amely nincs értelmezve valamely t_{0} időponttól kezdve. Az általános levezetésnél csak annyit teszünk fel, hogy az egyenlet autonóm, később azonban konkrétan reakciókinetikai egyenletekre alkalmazva megvizsgáljuk a polinomiális differenciálegyenlet esetét. Autonóm esetben az egyenlet a következő alakú:

ẋ(t) = f(x(t)), x(t_{0}) = c,  (1.1)

ahol J ⊂ R egy nyílt intervallum, f ∈ C(J, R), c ∈ J. Ekkor k_{1}, k_{2} ∈ Domx(t) esetén:

k_{2} − k_{1} = ∫_{k_{1}}^{k_{2}} 1dt = ∫_{k_{1}}^{k_{2}} \frac{ẋ(t)}{f(x(t))}dt = ∫_{x(k_{1})}^{x(k_{2})} \frac{1}{f} ds,  (1.2)

ami akkor értelmes, ha f az integrálás határa által kijelölt tartományban nem nulla.

Ha a megoldás korlátos ott, ahol értelmezve van, akkor biztosan értelmezve van az egész [t_{0}, +∞) félegyenesen, hiszen az általános elméletből ismert, hogy a megoldás minden kompakt halmazt elhagy.

Ha a megoldás nem korlátos, akkor k_{2}−k_{1} tetszőlegesen nagy lehet. Ha emellett \frac{1}{f} improprius integrálja valamilyen t > t_{0}, t ∈ R mellett a [t, +∞) félegyenesen létezik, akkor rögzített k_{1}-re x(k_{1}) > t és végtelenhez tartó k_{2} mellett az (1.2) egyenlőség bal oldalának határértéke végtelen, a jobb oldalé viszont korlátos. Ez ellentmondásra vezet, ha az előző gondolatmenet helytálló, azaz, ha a megoldás [t_{0}, +∞) -en értelmezve van és emellett x(k_{2}) > t bármely k_{2}-re. Ha tudjuk, hogy az x(t) megoldás valamilyen t_{1}-re pozitív és t > t_{1}-re monoton növő, akkor mindig pozitív is lesz, az integrál tehát x(t_{1})-től értelmes, és a nemkorlátos tulajdonság miatt ellentmondásra jutunk az előzőhöz hasonlóan. A reakciókinetikában mindig polinomiális differenciálegyenletek szerepelnek, tehát f egy polinom. Vagyis, ha a differenciálegyenlet jobb oldala negatív főegyütthatós polinom, akkor a megoldás nyilván korlátos is, ha pozitív főegyütthatós, akkor pedig elsőfokú esetben exponenciális növekedést tapasztalhatunk, magasabb fokszám esetén pedig egy bizonyos kezdeti feltételtől kezdve felrobban a rendszer. A kritikus kezdeti feltételre nézve nyilván elégséges feltétel az, hogy az f polinom onnan kezdve már csak pozitív értékeket vegyen fel, (így garantálható a megoldás monoton növő tulajdonsága) azaz a kezdeti feltételben szereplő érték nagyobb legyen a polinom legnagyobb gyökénél.

Tehát az előzőek alapján konkrétan így adható eljárás a rendszer felrobbanásának megállapítására egy anyagfajta (de persze belőle alkotott akármilyen komplexek) esetén: Tekintsük a leghosszabb olyan komplexet, amelyből még másik komplex keletkezik (feltehető, hogy ez nem az X komplex maga, mert akkor legfeljebb exponenciális növekedést tapasztalhatunk). Ekkor ha az erre a komplexre vonatkozó sebességi együttható pozitív, akkor a rendszer felrobban, ha negatív, akkor pedig nem. Amennyiben a kapott együttható nulla, akkor az eggyel kisebb komplexet tekintjük. Eredményeinket a következő tételben foglaljuk össze.

1. Tétel. Egy anyagfajtából álló kémiai rendszer esetén egy tömeghatás típusú kinetikával determinisztikus módon modellezett reakció pontosan akkor robban fel, ha a neki megfeleltetett kinetikai differenciálegyenlet legalább másodfokú és pozitív főegyütthatós. Az előző feltétel ekvivalens azzal, hogy a reakcióban résztvevő leghosszabb komplex legalább 3 hosszú, és a reakció során mennyisége nem fogy.

Példa Tekintsük a következő reakciót:

2Y → 3Y Y → 0.

Ennek a reakciónak a kinetikai differenciálegyenlete a következő

ẏ(t) = y^{2}(t) − y(t)

A tétel értelmében itt a másodfokú jobboldal pozitív főegyütthatója miatt felrobbanást fogunk tapasztalni, amit az y(0) = 2 kezdeti feltétel mellett az 1. ábrán szemléltetünk. (A kezdeti feltételt a fentiek értelmében szintén úgy választottuk, hogy teljsesül a felrobbanás elégséges feltétele.)

Az előzőek alapján felmerül a kérdés, azt vajon előre meg tudjuk-e mondani, hogy a megoldás, amely egy rögzített pontból indul, legfeljebb meddig van értelmezve, azaz mikor robban fel a rendszer. A fenti egyenletekből kapjuk, hogy

t − t_{0} = ∫_{x(t_{0})}^{x(t)} \frac{1}{f} dx = ∫_{c}^{x(t)} \frac{1}{f} dx,

vagyis ha több idő telik el, mint \frac{1}{f} -nek a c-től vett improprius integrálja, akkor biztos felrobban a rendszer.

Megjegyzés: A fentiek alapján világos az is, hogy az ẋ(t) = x^{2}(t) egyenletnek nincs olyan pozitív kezdeti értékkel rendelkező megoldása, mely az egész [t_{0}, +∞) félegyenesen értelmezve volna. Ennek alapján néhány speciális, egyszerű típusú, két anyagfajtából álló rendszer esetén is igazolható a felrobbanás.

Vizsgáljuk tehát a következő alakú differenciálegyenletrendszert:

ẋ(t) = ax(t)^{2} + bx(t)y(t) + cy(t)^{2}
ẏ(t) = dx(t)^{2} + ex(t)y(t) + fy(t)^{2}

Ha x(t) ill. y(t) helyett \frac{x(t)}{q} ill. \frac{y(t)}{s} ismeretlen függvényeket tekintjük, a fenti képletben szereplő differenciálegyenleteket ezekkel felírjuk, majd azokat p-vel ill. q-val szorozzuk, akkor azt kapjuk, hogy:

p(qẋ)(t) = paq^{2}x(t)^{2} + sbpqx(t)y(t) + cps^{2}y(t)^{2}
r(sẏ)(t) = rdq^{2}x(t)^{2} + erspx(t)y(t) + frs^{2}y(t)^{2}

Ha ezeket összeadjuk, a bal oldalon pqẋ(t) + rsẏ(t) áll, tehát ha emellett a jobboldalak összege pqx + rsy ennek négyzete, akkor a megoldások pozitív tulajdonságát is figyelembevéve a fenti megjegyzés alapján a rendszer felrobbanására következtethetünk. A jobb oldalak összege:

x^{2}(t)(paq^{2} + rdq^{2}) + x(t)y(t)(pqbs + qsre) + y^{2}(cps^{2} + frs^{2})  (1.3)

Azaz elsőként qs(pb + re) = 2pqrs azaz pb + re = 2pr feltételnek kell teljesülni. Másrészt a következő egyenlőségeknek kell még fennállniuk:

rd + pa = p^{2}
rf + cp = r^{2}

Itt természetesen p ill. r értékeit akarjuk úgy megválasztani, hogy a felírt egyenlőségek valóban teljesüljenek. Összefoglalva az eddigieket azt mondhatjuk, hogy a vizsgált alakú rendszer felrobbanásának elégséges feltétele, hogy a következő (túlhatározott) egyenletrendszer megoldható legyen.

rd + pa = p^{2}
rf + cp = r^{2}
pb + re = 2pr

Ez a feltétel elég kevés esetben teljesül, ezért egy konkrét differenciálegyenletrendszer esetében is igazoljuk a felrobbanást.

2. Tétel. A következő differenciálegyenletrendszer komponensenként pozitív kezdetiértékkel rendelkező megoldása felrobban.

ẋ(t) = y^{2}(t)
ẏ(t) = x^{2}(t)

Bizonyítás: Ha van olyan t_{1} > t_{0}, hogy x(t_{1}) = y(t_{1}), akkor az egyenletben meglevő szimmetria, valamint a megoldások egyértelműsége alapján azt kapjuk, hogy x(t) = y(t) minden t > t_{1} esetén, tehát akkor az x(t_{1}) = c_{1} kezdeti feltétellel a következő egyenlet teljesül: ẋ(t) = x^{2}(t) , ami valóban igazolja a felrobbanást. Ha nem létezik olyan t_{1} > t_{0}, ami a fenti tulajdonsággal rendelkezik, akkor feltehető, hogy x(t) > y(t) teljesül minden t > t_{0} esetén. Ekkor viszont t > t_{0}-ra ẏ(t) = x^{2}(t) = y^{2}(t) + f(t), ahol f(t) > 0. Ez megint a felrobbanást jelenti.

Következmény: A következő típusú differenciálegyenletek által leírt rendszerek is felrobbannak, amennyiben komponensenként pozitív kezdeti feltétellel rendelkeznek.

ẋ(t) = y^{2}(t) + ax(t) + by(t) + cx(t)y(t)
ẏ(t) = x^{2}(t) + dx(t) + ey(t) + fx(t)y(t),

ahol az a, b, c, d, e, f együtthatók nemnegatívak.

A fentiek alapján példát mutatunk olyan reakcióra, amely felrobban.

Példa:

2Y → 2Y + X  2X → 2X + Y
X + 2Y → 2X + 2Y  X + Y → X + 2Y → 2Y

A fenti reakció kinetikai differenciálegyenlete a következő:

ẋ(t) = y^{2}(t)
ẏ(t) = x(t)y(t) + x^{2}

A fejezet első részében levezetett eredményből az is látszik, a rendszer felrobbanásának szükséges feltétele egy dimenzióban, hogy legyen olyan komplex, amelynek hossza legalább három. Ez ugyanígy érvényes két dimenzióban is, amint azt a következő tétel mutatja.

3. Tétel. Két anyagfajtából álló kémiai rendszer esetén egy tömeghatás típusú kinetkával determinisztikus módon modellezett reakció felrobbanásának szükséges feltétele, hogy legalább három hosszú legyen a leghosszabb komplex.

Bizonyítás: Indirekt bizonyítást alkalmazunk. Tegyük fel, hogy csak legfeljebb kettő hosszú komplexek vannak jelen, azaz 0, X, Y, 2X, 2Y, X + Y. Be fogjuk látni, hogy ekkor a rendszer nem robban fel. Foglalkozzunk egyelőre azokkal a reakciókkal, amelyek a kinetikai differenciálegyenletben szereplő pozitív együtthatóval rendelkező nemlineáris tagokat adhatják. Eleve csak olyan reakciókról lehet szó, amelyekben legalább az egyik oldalon szerepel egy kettő hosszú komplex, ezek közül pedig a a következők:

2X ⇀ ↽ 2Y  2X ⇀ ↽ X + Y  2X ⇀ ↽ Y
X ⇀ ↽ 2Y  2Y ⇀ ↽ X + Y

Egyszerű meggondolni, hogy pontosan ezekről a reakciókról van szó. A sebességi állandók legyenek rendre a_{1}, a_{2}, b_{1}, b_{2}, c_{1}, c_{2}, d_{1}, d_{2}, e_{1}, e_{2}, amelyekről ebben az esetben csak a nemnegativitást követeljük meg. A kinetikai differenciálegyenlet(rendszer) ekkor a következő lesz:

ẋ(t) = 2a_{2}y^{2}(t) − 2a_{1}x^{2}(t) − b_{1}x^{2}(t) + b_{2}x(t)y(t) − 2c_{1}x^{2}(t) + 2c_{2}y(t)

−d_{1}x(t) + d_{2}y^{2}(t) + e_{1}y^{2}(t) − e_{2}x(t)y(t)

ẏ(t) = 2a_{1}x^{2}(t) − 2a_{2}y^{2}(t) + b_{1}x^{2}(t) − b_{2}x(t)y(t) + c_{1}x^{2}(t) − c_{2}y(t)

+2d_{1}x(t) − 2d_{2}y^{2}(t) − e_{1}y^{2}(t) + e_{2}x(t)y(t)

Az egyenleteket összeadva a jobb oldalra azt kapjuk, hogy: −c_{1}x^{2}(t) + c_{2}y(t) + d_{1}x(t) − d_{2}y^{2}(t). Ezen kívül még lineáris tagok is lehetnek. Ebből, az előző kifejezésből és az egyenletrendszer bal oldalából kapjuk, hogy: ẋ(t) + ẏ(t) ≤ A(x(t) + y(t)), ahol A egy elég nagy szám. (Itt ismét felhasználjuk azt az általános elméletből ismert tényt, hogy a megoldások mindig pozitívak.) x(t) + y(t) tehát nem robbanhat fel, ami ismét a pozitivitást használva azt jelenti, hogy sem x(t), sem y(t) nem robbanhat fel, amivel a tételt indirekt módon igazoltuk.

Megjegyzés: Többdimenziós esetben a felrobbanás kérdése jóval bonyolultabb.

A nemlineáris rendszerek (nyilván eleve csak ezek érdekesek) elemzésére szolgáló kvalitatív elmélet itt nem alkalmazható közvetlenül, hiszen a fázisdiagrammokról csak azt látjuk ránézésre, hogy a megoldásnak mi a határértéke, az viszont nem derül ki, hogy véges időn belül éri ezt el (ami a felrobbanást jelenti) vagy sem.

1.1.2. A sztochasztikus modellezés vizsgálata a felrobbanás problémájának szempontjából

A továbbiakban az [1] cikk két kitételét pontosítjuk arra az esetre, amit ebben a fejezetben megvizsgálunk. Felidézzük ezért a cikk jelöléseit és fontosabb eredményeit. A reakció alapján felírt kinetikai differenciálegyenlet főegyütthatóját jelöljük c_{d}-vel. A sztochasztikus modell az általunk vizsgált esetben a globális modellt jelenti, vagyis a reakciót egy születési-halálozási folyamatnak tekintjük, és a cikk jelöléseivel:

λ_{k} = L(a_{0} + a_{1} \frac{k}{L} + a_{2} \frac{k(k − 1)}{L^{2}} + ... + a_{p} \frac{k(k − 1)...(k − p + 1)}{L^{p}})

μ_{k} = L(b_{1} \frac{k}{L} + b_{2} \frac{k(k − 1)}{L^{2}} + ... + b_{q} \frac{k(k − 1)...(k − q + 1)}{L^{q}}),

ahol λ_{k} ill. μ_{k} a folyamatban szereplő együtthatók, amelyek azt fejezik ki, hogy milyen valószínűséggel nő ill. csökken a rendszerben levő részecskék száma az idővel arányosan. Ugyanez képlettel:

P (X(t + s) = k + 1|X(s) = k) = λ_{k}t + o(t)
P (X(t + s) = k − 1|X(s) = k) = μ_{k}t + o(t)

A képletben szereplő L érték annak a térrésznek a térfogatát jelöli, ahol a reakció végbemegy. Ezt egységnyinek fogjuk választani, ugyanis ez a későbbi vizsgálatainkban semmiféle szerepet nem játszik. Végül a_{i}, ill. bi jelöli a reakcióban szereplő, az i számú részecskéből álló komplex keletkezési, ill. mgsemmisülési sebességét. Jelölje P_{k}(t) annak a valószínűségét, hogy a t időpillanatban épp k darab részecske van jelen. A második Kolmogorov-egyenletből a fentiek felhasználásával a következő alapegyenlet írható fel: Ṗ_{k}(t) = λ_{k−1}P_{k−1}(t) + μ_{k+1}P_{k+1}(t) − (λ_{k} + μ_{k})P_{k}(t) a következő kezdeti feltétellel együtt:

∑_{k=1}^{∞} P_{k}(0) = 1,

ahol ezek a P_{k}(0) értékek adottak. Ez azt jelenti, hogy ismerjük a kezdetben meglevő tömeg (vagy részecskeszám) eloszlását. Az említett cikkben annak igazolása szerepel, hogy ha determinisztikus modell esetére felírt kinetikai differencálegyenlet főegyütthatója negatív, akkor a sztochasztikus és a determinisztikus modell bizonyos értelemben ekvivalens. Az alábbiakban annak bizonyítását szemléltetjük, hogy a felrobbanás szempontjából a kétfajta modell ekvivalens akármilyen egy anyagfajtából álló rendszer esetében. Mindenekelőtt azonban tisztázni kell, hogy mit jelent a felrobbanás ebben a modellben. Definiáljuk ezt úgy, hogy a részecskék számának várható értéke véges időn belül végtelenné válik, azaz a várható értékre felírt differenciálegyenlet felrobban. Két konkrét példán szemléltetjük, hogy lehet egy ilyen differenciálegyenletet felírni, majd egy általános tételre hivatkozva igazoljuk a fent megfogalmazott állítást.

1. Példa Tekintsük először a X → 2X reakciót (egységnyi nagyságú sebességi állan-

dóval), amelyről tudjuk, hogy kinetikai differenciálegyenlete a determinisztikus modell szerint a következő alakú: ẋ(t) = x(t), vagyis a főegyüttható pozitív, és az egyenlet az exponenciális növekedést jelenti. Az alábbiakban azt fogjuk igazolni, hogy a születési-halálozási folyamattal modellezve sem robban fel ez a rendszer, sőt a várható értékekre kapott differenciálegyenlet épp ugyanaz lesz, mint fent. A fenti reakció esetében a folyamat együtthatóira vonatkozó képlet alapján a kapjuk, hogy λ_{k} = k μ_{k} = 0, vagyis az alapegyenletre a következő adódik: Ṗ_{k}(t) = (k − 1)P_{k−1}(t) − kP_{k}(t). A t időpillanatban jelen levő részecskék számát X(t)-vel jelölve (ez nyilván egy valószínűségi változó) azt írhatjuk, hogy

E(X(t)) = ∑_{k=1}^{∞} kP_{k}(t),

amiből az egyenlet két oldalát deriválva, majd az alapegyenletet felhasználva kapjuk, hogy

Ė(X(t)) = ∑_{k=1}^{∞} kṖ_{k}(t) = ∑_{k=1}^{∞} k(−k)P_{k}(t)+(k+1)kP_{k}(t) = ∑_{k=1}^{∞} kP_{k}(t) = E(X(t)).

Az átalakítás ötletét a [5] könyvből vettem át (144–145. oldal), és a módszert a másik példa esetén is alkalmazni fogjuk. A részecskék számának t időpillanatban vett várható értéke ezek szerint mindenütt értelmezett, tehát valóban nem robban fel a rendszer a sztochasztikus modell szerint sem.

2. Példa Tekintsük most a 2X → 3X reakciót. Most a_{2} = 1, amiből a folyamatot leíró együtthatókat kiszámolva, majd az alapegyenletbe beírva kapjuk, hogy Ṗ_{k}(t) = (k − 2)(k − 1)P_{k−1}(t) − k(k − 1)P_{k}(t).

Írjuk fel az előzőhöz hasonlóan először a várható érték deriváltját, azaz

Ė(X(t)) = ∑_{k=1}^{∞} kṖ_{k}(t) = ∑_{k=1}^{∞} k(−k)(k − 1)P_{k}(t) + (k + 1)k(k − 1)P_{k−1}(t) =

= ∑_{k=2}^{∞} −k^{2}(k−1)P_{k}(t)+(k+1)k(k−1)P_{k}(t) = ∑_{k=2}^{∞} k(k−1)P_{k}(t) = ∑_{k=1}^{∞} k(k−1)P_{k}(t).

Másrészt tudjuk, hogy

E(X^{2}(t)) = ∑_{k=1}^{∞} k^{2}P_{k}(t).

Vagyis az előzőeket összevetve

E(X^{2}(t)) − Ė(X(t)) = ∑_{k=1}^{∞} kP_{k}(t) = E(X(t)),

amiből

E(X^{2}(t)) = D^{2}(X(t)) + E^{2}(X(t)) felhasználásával azt írhatjuk, hogy Ė(X(t)) = D^{2}(X(t)) + E^{2}(X(t)) − E(X(t)), ami viszont egy másodfokú, pozitív főegyütthatós differenciálegyenlet E(X(t))- re. A jobb oldalon szerepel még egy nemnegatív tag is (a szórásnégyzet), de ez a felrobbanás tényén nem változtat. Vagyis a determinisztikus esethez hasonlóan itt is felrobbanást igazoltunk.

A fenti levezetések általánosabban is elvégezhetők a generátorfüggvény alkalmazásával. A generátorfüggvény segítségével a részecskék (konkrét esetben pl. molekulák) számának, azaz X-nek várható értékére egy differenciálegyenlet írható fel, ahogyan azt két konkrét példán szemléltettük. Az így kapott differenciálegyenletek is polinomiálisak lesznek, és a jobboldalon szereplő polinom foka, ill. főegyütthatójának előjele megegyezik a determinisztikus módon tömeghatás típusú kinetikával modellezett reakcióra felírt differenciálegyenlet jobb oldalán szereplő polinom hasonló paramétereivel. Ekkor viszont az első tétel értelmében a két differenciálegyenlet a felrobbanás tulajdonságára nézve ugynúgy viselkedik (olyan értelemben, hogy vagy mindkettő felrobban vagy egyik sem). A két vizsgált modell tehát ilyen értelemben tényleg ekvivalens. A pontos fogalmakat, képleteket, levezetéseket ld. [7] 130–131. oldalán.

1.2. Sebességi állandók pontosabb becslésének lehetősége a sztochasztikus modell esetében

Az előző fejezet alapján úgy tűnhet, hogy nem is érdemes sztochasztikusan modellezni a kémiai reakciókat, hiszen a koncentráció várható értéke hosszú távon nagyjából úgy viselkedik, mint a determinisztikus modell alapján kiszámolt koncentrációérték. Éppen ezért ebben a részben a sztochasztikus modell, azaz a reakció születési-halálozási folyamattal való leírásának egy előnyös tulajdonságára mutatunk rá. Ismét az egyváltozós esetre fogunk szorítkozni, vagyis egyetlen anyagfajta van jelen. A születési-halálozási folyamat fent bevezetett λ_{k}, μ_{k} együtthatóira tett bizonyos feltevések mellett létezik ún. stacionárius eloszlás, amely analóg a determinisztikus modellben szereplő differenciálegyenletek egyensúlyi pontjaival. Egy könnyen ellenőrizhető elégséges feltételt fogalmazunk meg. Pontosabban a következők állnak fenn (a feltételeket illetően ld. a [4][242–247. oldal], a levezetést illetően pedig a [3][236–237. oldal] könyvet). Vezessük be a következő jelöléseket

π_{0} = 1, π_{k} = \frac{λ_{0}λ_{1} … λ_{k−1}}{μ_{1}μ_{2} … μ_{k}} .

Ha teljesülnek a

∑_{k=0}^{∞} \frac{1}{λ_{k}π_{k}} = +∞,  ∑_{k=0}^{∞} π_{k} < +∞

feltételek, továbbá a folyamat által felvett értékek nemnegativitását biztosító

λ_{0} > 0, μ_{0} = 0,
ill. minden pozitív egész k-ra λ_{k} > 0, μ_{k} > 0,

akkor létezik stacionárius eloszlás, tehát egy olyan eloszlás, amelyre fennáll, hogy

lim_{t→∞} p_{n}(t) = lim_{t→∞} p_{n|k}(t) = p_{n}^{∗}.

Az eloszlás egyes tagjait tehát p_{n}^{∗}-al jelöljük, ahol a kitevőben levő csillag arra utal, hogy stacionárius eloszlásról van szó. p_{n|k}(t) pedig annak a valószínűsége, hogy a t időpontban n darab részecske lesz, ha kezdetben (a nulla időpillanatban) k darab volt.

A stacionárius eloszlás kiszámítása analóg a differenciálegyenletek egyensúlyi pontjainak meghatározásával, csak a differenciálegyenlet helyett most az előremutató egyenletben szereplő deriváltat kell nullává tenni. Megoldásként azt kapjuk, hogy

p_{n}^{∗} = p_{0}^{∗} ∏_{k=0}^{n−1} \frac{λ_{k}}{μ_{k}} .

Ha egy kémiai reakció olyan, hogy a sebességi állandóiból alkotott λ_{k} és μ_{k} együtthatók teljesítik a fenti feltételeket, akkor tehát a koncentráció eloszlása egy idő után állandó lesz (pontosabban megfelelő értelemben tetszőlegesen közel kerül egy meghatározott eloszláshoz). Ezt használjuk ki akkor, amikor a stacionárius eloszlás ismeretében megbecsüljük a sebességi állandókat. A stacionárius eloszlás „ismerete” a gyakorlatban azt jelentheti, hogy a részecskeszámra (ill. koncentrációra) vonatkozó mérések alapján a tapasztalati eloszlást ismerjük. Az eddigieket egy konkrét példán szemléltetjük. A reakció legyen a Schlögl-reakció, azaz a következő:

A ⇀ ↽_{k_{2}}^{k_{1}} X A + 2X ⇀ ↽_{k_{4}}^{k_{3}} 3X.

Itt A egy olyan anyagfajta, melynek koncentrációja állandónak vehető, tehát a mérések alapján egy ismert paraméter. Először a születési-halálozási folyamat együtthatóit határozzuk meg.

λ_{k} = L(Ak_{1} + Ak_{3} \frac{k(k − 1)}{L^{2}}) μ_{k} = L(k_{2} \frac{k}{L} + k_{4} \frac{k(k − 1)(k − 2)}{L^{3}})

Megjegyzés: Ismert tény, hogy a fenti típusú (megfordítható és egy anyagfajtából álló) reakcióknak van stabil egyensúlyi pontja (sőt, ennél valamivel több is igaz), vagyis „számíthatunk rá”, hogy a sztochasztikus modellezés szerint is létezik stacionárius eloszlás.

Legyen egyelőre k_{4} az ismeretlen sebességi együttható, L-et pedig, a térrész hosszát válasszuk ismét egységnyinek. A fentieket figyelembevéve a stacionárius eloszlás a következő lesz:

p_{n}^{∗} = p_{0}^{∗} ∏_{k=0}^{n−1} \frac{k_{3}Ak(k − 1) + k_{1}A}{k_{2}(k + 1) + k_{4}(k + 1)k(k − 1)}

Tehát az ismeretlen k_{4} sebességi állandó a következőképpen becsülhető. A tapasztalati eloszlás, azaz a koncentrációra (vagy részecskeszámra) vonatkozó megfelelő számú (sok) mérés ismeretében olyan k_{4} értéket keresünk, amelyhez tartozó stacionárius eloszlás a leginkább hasonlít a megfigyelt eloszláshoz. Hasonló módon becsülhető akár három sebességi állandó is (egyszerre). Akkor volna gond a becsléssel elvileg, ha egyazon eloszláshoz több sebességi együtthatókból álló hármas is tartozhatna. A stacionárius eloszlásra vonatkozó formula alapján ugyanis nyilvánvaló, hogy ha két ilyen megegyezik, akkor az azt előálllító szorzatban szereplő minden törtnek meg kell egyeznie. Ha például k_{3} értéke rögzített, akkor a többi sebességi együtthatót értelemszerűen k_{11}, k_{12}, k_{14}-el, ill. k_{21}, k_{22}, k_{24}-el jelölve a törtek összevetésével, szorzással a kapott polinomok legmagasabb fokú tagjaira azt kapjuk, hogy

k_{14}k_{3}k^{5}, ill. k_{24}k_{3}k^{5},

amiből azonnal látszik, hogy k_{14} = k_{24}. Most a polinomok konstans tagjait felírva

Ak_{14}(k_{21} + k_{3}), ill. Ak_{24}(k_{11} + k_{3}),

amiből az előző eredményt is figyelembevéve kapjuk, hogy k_{11} = k_{21}. Végül innen már könnyen látható, hogy a második sebességi állandók a megfelelő törtek egyenlősége miatt szintén megegyeznek.

Megjegyzés: Amennyiben a determinisztikus modellt használtuk volna a stacionárius pontban végzett eloszlások kiértékelésére, akkor a következő eredményt kaptuk volna. A kinetikai differencálegyenlet a következő lesz:

ẋ(t) = k_{1}A − k_{2}x(t) + k_{3}Ax^{2}(t) − k_{4}x^{3}(t).

Ha itt például k_{1} adott, ezenkívül ismerjük a stacionárius eloszlásnak megfelelő rögzített x(t) függvényértéket, melynél a bal oldal nulla, akkor a megfelelő számértékeket beírva is csak egy

0 = B + Ck_{2} + Dk_{3} + Ek_{4}

alakú egyenletet kapunk, amiből a sebességi együtthatókat nem tudjuk egyértelműen meghatározni. Látható az okoskodásból, hogy még két ismert állandó alapján sem tudjuk a többit meghatározni.

A becslés végrehajtásához pontosan meg kellene fogalmazni, hogy mit értünk azon, hogy a tapasztalati eloszláshoz minél közelebbi eloszlást keresünk, és az ez alapján kidolgozott módszert numerikusan, szimulált mintavétellel tesztelni is kell. Elsőként az az ötlet merül fel, hogy maximum likelihood becslést kellene alkalmazni, azaz olyan sebességi együtthatót választani, amelynek tényleges realizációja esetén a mért tapasztalati eloszlás valószínűsége a lehető legnagyobb. Azonban mivel az eloszlás nem véges tartójú, azonnal látszik a veszély, hogy szorzatként egymáshoz nagyon közeli, kis számokat kapunk majd, amely a mintavételben levő bizonytalansággal együtt a becslést használhatatlanná teszi. Egyszerűbb becslési eljárást kapunk, ha azt a stacionárius eloszlást (ill. a hozzá tartozó sebességi együtthatót) választjuk, amelyiknek maximumnormában vagy L_{2} normában vett eltérése a tapasztalati eloszlástól a lehető legkisebb. Az eljárások numerikus teszteléséhez a következőképpen járunk el. Választunk egy sebességi állandót, és azzal a most ismert stacionárius eloszlást elkészítjük. Némi nehézséget okoz, hogy az első tagot semmilyen formula sem adja meg, de ez a többi adatból a programban látható módon jó közelítéssel (az eloszlás értékeinek exponenciális csökkenését kihasználva) kiszámolható. Most mintát veszünk az eloszlásból, a valódi méréseket szimuláljuk. Ezután úgy tekintjük, hogy csak ez a minta adott, és ennek alapján akarjuk az eredetileg már ismert sebességi állandót megbecsülni. A becsléshez mindhárom fenti változatot kipróbáltam, először csak egy sebességi állandóra, de egyik sem működik jól, illetve csak akkor, ha a minimalizálandó paramétert már eleve a tényleges érték közeléből indítom el. Az említett programok MATLAB-ban íródtak, schloglmax.m ill. schlogl .m néven. Ld. még az ottani kommentárt is.

Megjegyzés: A stacionárius eloszlás maximumának vizsgálata újabb érdekes kapcsolatot tár fel a determinisztikus és a sztochasztikus modellezés között. Erről is szó van mindkét idézett munkában ( [4, 248–249. oldal], ill. [3, 242–245. oldal]).

2. fejezet

2.1. A sztochasztikus rezonancia jelensége

A tudományos kutatásban gyakran találkozunk olyan jelenségekkel, melyek megfigyelése során a jelenségre jellemző bizonyos adatokat akarunk meghatározni. Lássunk néhány konkrét példát:

1. Példa Egy térfogatrészben részecskék mozognak. Meg akarjuk mérni a hőmérsékletet, ami nyilván a részecskék sebességének függvénye.

2. Példa Egy kémiai reakcióban a résztvevő anyagok koncentrációjának vektora a t időpontban az Y_{t} sztochasztikus folyamat t időpillanatban felvett értékével modellezhető. Ennek az egyik sebességi állandó valamilyen (lehet, hogy implicit) függvénye. A feladat ezen sebességi állandó meghatározása. A fentiek alapján általánosabban úgy fogalmazhatunk, hogy az Y_{t} sztochasztikus folyamatból vett megfigyelések alapján (vagyis a trajektória egy részhalmazának ismeretében) meg akarjuk becsülni a folyamat valamilyen paraméterének egy bizonyos függvényét. Jelöljük X-szel a megfigyelésekből álló vektort. Elképzelhető, hogy ha nem magát a folyamatot, hanem annak egy zajjal terhelt módosítását figyeljük, akkor a becslés jobb lesz, sőt megfelelően választott zaj esetén esetleg olyan jelenségről is információt kaphatunk, amelyet egyébként nem is tudnánk mérni. A megfigyelt jel tehát a zavaró hatások által felerősödhet (olyan értelemben is, hogy pontosabban meghatározható). Ezt a jelenséget nevezzük \underline{sztochasztikus rezonanciának}. Gyakran az Y_{t} folyamat valójában determinisztikus, és csak a módosított folyamatot kezeljük sztochasztikus folyamatként. Példa az előzőekre: A vízben élő rákok különösen érzékenyek környezetük periodikus mozgásaira. Gyakran azonban olyan gyenge ingerek is hatást váltottak ki megfigyelésükkor, amelyek érzékelésére a receptorok elvileg képtelenek lennének. A jelenség magyarázata az, hogy egy véletlen zaj hatására az ingerküszöbhöz közeli értékek esetén a zaj és a jel együttesen már meghaladták az érzékelhető küszöböt. Másik példa: az alagútjelenség. Egy tömegpont halad egy „potenciálfal” felé. Amennyiben kevesebb energiával rendelkezik, mint amekkora a falon való áthaladáshoz szükséges, akkor elvileg (a klasszikus mechanikai modell alapján) nem haladhat át rajta. Mivel azonban kísérletileg kimutatható, hogy megjelenik a fal túloldalán is, azért gondolhatjuk itt is úgy, hogy zaj hatására az esetek egy részében mégiscsak sikerül átlépnie. Ezt a zajt a modellbe már be is építették, így az új (kvantummechanikai) modell már (jobban) tükrözi a valóságot. A fentiekben azt említettük, hogy a becslés „jobb lesz”. Értelmezzük ezt úgy, hogy az átlagos veszteség, tehát E(f(X)−c)^{2} minimális legyen, ahol f(X) az X minta (mint valószínűségi változó) függvényében (tehát a megfigyelés alapján) a keresett c értékre adott becslés. A továbbiakban azzal a kérdéssel foglalkozunk, hogy milyen zajt adjunk valamilyen konkrét folyamathoz, hogy a keresett paraméter becslése a fenti értelemben a lehető legjobb legyen.

2.1.1. A sztochasztikus rezonanciával kapcsolatos konkrét modell vizsgálata

A továbbiakban a [6] cikkben szereplő jelenségnek adjuk meg egy egyszerű, lehetséges modelljét, melynek alapján választ kaphatunk a cikkben szereplő jelenségek hátterére is. A szaglás mechanizmusának tanulmányozása során a ténylegesen jelenlevő koncentráció függvényében vizsgálták a zaj mértékét úgy, hogy növelték a koncentrációt. Azt a jelenséget figyelték meg, hogy az érzékelőrendszerben jelen levő zaj kezdetben nagyjából a koncentráció negyedik hatványával arányosan növekszik, azonban ez a növekedés egy idő után lelassul, majd meg is áll, a zaj egy küszöbértéket nem lép túl. Szemléletesen úgy magyarázhatjuk ezt, hogy ha csak a megszokott hatás éri a receptorokat, akkor nincs értelme nagy zajt hozzáadni a tényeges értékhez, hiszen ez esetben akkor is érzékelnénk néha valamit, amikor semmilyen tényleges hatás nem volt. Ha azonban valamilyen átlagos értéket meghaladó hatás mutatkozik, akkor célszerű zajt hozzáadni, hogy érzékelni tudjuk a gyengébb, de ténylegesen meglevő szagot. A fent általánosan megfogalmazott modell tehát a következő. A szagingernek egy sztochasztikus folyamat felel meg, az érzékelésnek, ill. az ebből kialakuló ingernek pedig a folyamat valamely paraméterére (konkrétan majd minden időpillanatban a várható értékére) adott becslés. Abból az alapfeltevésből indulunk ki, hogy a mechanizmus a lehető legjobban működik, azaz olyan zajt ad a jelhez, hogy a hatás érzékelése (azaz a becslés) a lehető legjobb legyen. Jelölje az ingerküszöb értékét z, a tényleges hatás erősségét pedig c. A fentieket képzelhetjük áramerősségértékeknek valamilyen mértékegységben. Az egyszerűség kedvéért tegyük fel, hogy c egy állandó érték (később azonban enyhítünk majd ezen a feltevésen).

Az egyenletes eloszlású zaj esete

Első esetben a tényleges értékhez egyenletes eloszlású zajt adunk hozzá. Jelöljük Y -nal az így kapott valószínűségi változót (itt a sztochasztikus folyamat minden pillanatban egy ilyen Y valószínűségi változó). Világos, hogy Y ∼ c + U(−s, s), ahol s az ismeretlen, keresett optimalizálandó paraméter, U(−s, s) pedig a [−s, s] intervallumon értelmezett egyenletes eloszlás szokásos jelölése. Nyilván Y > z esetében fogjuk érzékelni azt, hogy a jel átlépte a küszöböt. Tudjuk azt is, hogy

P (Y > z) = P (c+U(−s, s) > z) = P (U(−s, s) > z−c) = P (U(−1, 1) > \frac{z − c}{s}) =

= \frac{1 − \frac{z−c}{s}}{2} = \frac{1}{2} − \frac{z − c}{2s}.

Vegyük észre, hogy a bal oldalra vonatkozó torzítatlan becslést ismerünk, ugyanis a megfelelő Y > z esemény indikátorának várható értéke a bal oldal, melyre vonatkozólag a mintaátlag, vagyis a relatív gyakoriság valóban torzítatlan becslés. Az előzőek alapján fel is írhatjuk a c-re vonatkozó becslést, amely a

\frac{K}{N} = \frac{1}{2} − \frac{z − c}{2s}

összefüggésből számolással

z − c = s − \frac{2sK}{N} ,

azaz c értékének becslése

ĉ = z − s + \frac{2sK}{N}

lesz, ahol N a mérések száma, K pedig azon esetek száma, amikor a jel meghaladta az ingerküszöböt. Helyettesítsük be most ezt az eredményt a veszteségfüggvényt megadó formulába. Vegyük ugyanakkor figyelembe, hogy a fenti okoskodás alapján már az indikátorból (és nem a mintaátlagból) álló becslés is torzítatlan lesz (ezt az egyszerű állítást később bizonyítjuk), vagyis a veszteségfüggvény nem más, mint a becslés szórásnégyzete. Amikor pedig a mintaátlagot tekintjük, akkor az N számú független indikátor összegeként írható, ezért

D^{2}(z − s + \frac{2sK}{N}) = D^{2}(∑_{i=1}^{N} \frac{z − s}{N} + \frac{2sχ(Y > z)}{N}) =

= \frac{1}{N} D^{2}(z − s + 2sχ(Y > z)) = \frac{1}{N} D^{2}(s(1 − 2χ(Y > z)) =

= \frac{1}{N} s^{2}(E(1 − 4χ(Y > z) + 4(χ(Y > z))^{2} − (1 − 2c)^{2}) = \frac{1}{N} s^{2}(1 − (1 − 2c)^{2}).

Ezt s-ben minimalizálva azt kapjuk, hogy s-nek minél kisebbnek kell lennie úgy, hogy még érvényes legyen az előző okoskodás, azaz a becslés torzítatlan maradjon. Konkrétan ez azt jelenti, hogy ne legyen nulla annnak a valószínűsége, hogy a mért érték eléri az ingerküszöböt, vagyis mindig igaznak kell lennie a s ≥ c egyenlőtlenségnek. Ez pedig akkor teljesül, ha s = c. Szemléletesen ez azt jelenti, hogy éppen akkora zajt kell mesterségesen hozzáadni, hogy a kapott, zajjal terhelt folyamat lehetséges maximuma éppen elérje a kritikus z értéket. Ez az eredmény azonban „instabil”, mert az optimumot kissé megváltoztatva már nagyon gyenge lesz a becslés, ugyanis az előírtnál kisebb zaj esetén a hatásokat egyáltalán nem lehet érzékelni (pedig éppen ez lett volna az eljárás lényege), mivel akkor a folyamat értéke soha nem éri el az érzékelhető küszöböt. Jogosan merül fel a kérdés, hogy nem azért kaptunk-e ilyen furcsa eredményt, mert nem a lehető legjobb becslést állítottuk elő. Azonban igazolni fogjuk, hogy a fenti becslés torzítatlan, emellett a \frac{K}{N} minimális elégséges statisztika függvénye, így a fent előállított becslés az egyetlen hatásos becslés, azaz a legjobb, amit megfigyeléseink alapján megkaphattunk.

Elsőként a fent már említett torzítatlanságot igazoljuk.

E(ĉ) = E(z − s + \frac{2sK}{N}) = z − s + 2sE(\frac{K}{N}) = z − s + 2sP (Y > z),

amiből átrendezéssel és a P (Y > z) értékére felírt formula behelyettesítésével kapjuk, hogy E(ĉ) = c, amit igazolni akartunk.

Második lépésként megmutatjuk, hogy \frac{K}{N} a c paraméterértékre nézve elégséges statisztika. A mérések összesített eredménye, vagyis a minta felfogható úgy, mint egy nullaegy sorozat, melyben az egyesek száma éppen K. K eloszlása a mérések függetlensége miatt binomiális, az eloszlás rendje N, paramétere pedig P (Y > z). Így egy tetszőleges mintaelem valószínűségét megkaphatjuk a binomiális eloszlás ismert tagjainak függvényében (pontosabban fogalmazva az eloszlás K értékhez tartozó tagját le kell osztani a lehetséges K hosszúságú sorozatok számával). Mivel az így előállított mintaelemekre vonatkozó sűrűségfüggvény csak K értékén keresztül függ a paramétertől, ezért a Neyman-féle faktorizációs tételt is felhasználva kapjuk, hogy K (vagy \frac{K}{N}) valóban elégséges statisztika.

Harmadik lépésként azt igazoljuk, hogy a \frac{K}{N} statisztika teljes is a c értékére nézve. Ehhez azt kellene belátni, hogy ha egy olyan y függvényre, amelyre

E_{c}(y(\frac{K}{N})) = 0

minden c lehetséges értékre, akkor y azonosan nulla azokon a helyeken, amelyek a c lehetséges értékei által adott eloszlások tartójában vannak. Ezek minden esetben az egyes \frac{K}{N} értékek, az egyes c paraméterek által adott eloszlások pedig a különböző P (Y > z) paraméterű binomiális eloszlások. Tegyük fel, hogy a statisztika mégsem teljes, mert létezik olyan y valós értékű függvény, amely értelmezve van a \frac{K}{N} (K = 0 … N) helyeken, ott nem azonosan nulla, viszont minden lehetséges c értékre teljesül a

E_{c}(y(\frac{K}{N})) = 0

egyenlőség. Jelöljük egy rögzített c-re és j-re a c értékéből előállított binomiális eloszlás j-edik tagját b_{c,j}-vel, vagyis

b_{c,j} := \binom{N}{j}p^{j}(1 − p)^{N−j},

ahol

p = \frac{1}{2} − \frac{z − c}{2t}.

Ezzel a jelöléssel élve az előző feltétel a következőképpen írható:

E_{c}(y(\frac{K}{N})) = y(0)b_{c,0} + y(\frac{1}{N})b_{c,1} + … + y(1)b_{c,N} = 0

teljesül minden megengedett c paraméterre. Vegyük észre, hogy ez a binomiális eloszlás paraméterére, mint változóra nézve egy N-edfokú polinom. Ezért, ha ez c bármely értékére nulla, akkor biztosan az azonosan nulla polinomról van szó, hiszen c különböző értékeihez különböző paraméterértékek tartoznak, és így a polinom legalább N + 1 helyen (sőt, végtelen sok helyen is) nulla lesz. Tehát a polinom minden együtthatójának nullának kell lennie. A konstans tag együtthatója viszont csakis a fenti összeg első tagjából adódik, ami azt jelenti, hogy az y(0) = 0 egyenlőségnek fenn kell állnia. Ekkor viszont az elsőfokú tag csakis y(\frac{1}{N})b_{c,N−1}-ből adódik, ami azt jelenti, hogy y(\frac{1}{N}) is nulla. Hasonlóan látható, hogy az y függvény összes további \frac{j}{N} (j = 2 … N) helyen felvett értéke is nulla, ami ellentmond eredeti feltevésünknek. Ezzel igazoltuk a szóbanforgó statisztika teljességét is.

Végül arra a tényre hivatkozunk, hogy ha egy torzítatlan becslés függvénye egy teljes elégséges (így persze minimális elégséges) statisztikának, akkor az az egyetlen hatásos becslés. A fentiekben felhasznált fogalmak, állítások megtalálhatók például a [2] könyvben.

A normális eloszlású zaj esete

Vizsgáljuk meg most azt az esetet, amikor a zaj normális eloszlású. Módosítsuk az előző részben bevezetett jelöléseket, amelyekkel a következőképpen írható le a mostani probléma: Adott egy z − c várható értékű normális eloszlás, amelynek szórásnégyzetét optimális módon akarjuk meghatározni úgy, hogy ennek az Y ∼ N(z − c, σ^{2}) eloszlásnak a várható értékét az adott típusú megfigyelések alapján a lehető legjobban tudjuk becsülni. Ismét z értékét tekintjük méréshatárnak, azaz z pozitív, csakúgy, mint z−c. Ezt minden rögzített c-re meg akarjuk határozni, és az optimális paraméterérték c-től való függésére vagyunk kíváncsiak. A mintavétel ismét azt jelenti, hogy csak azt tudjuk, Y értéke meghaladt-e (egy ismeretlen) c értéket. Természetesen több erre vonatkozó mérés is rendelkezésre áll. Legyen ismét N a mérések száma, K pedig az észlelt átlépések száma. Tekinthetjük tehát a mérést egy X indikátor valószínűségi változónak, ahol

P (X = 1) = Φ_{σ^{2}}(−c) P (X = 0) = 1 − Φ_{σ^{2}}(−c),

ahol Φ_{σ^{2}} jelöli a σ^{2} szórásnégyzetű, nulla várható értékű normális eloszlás eloszlásfüggvényét. Használjuk a standard normális esetre a szokásos Φ jelölést is. Az előző részben követett okoskodás mintájára azt írhatjuk tehát, hogy

E(X) = P (X = 1) = Φ_{σ^{2}}(−c),

azaz az E(X)-re vonatkozó tapasztalati értékből és az eloszlásfüggvény inverzére vonatkozó

Φ_{σ^{2}}^{−1}(x) = σΦ^{−1}(x)

képletből azt kapjuk, hogy

−c = σΦ^{−1}(\frac{K}{N}),

vagyis −σΦ^{−1}(\frac{K}{N}) := f(X) lesz c becslése. Azonban felmerül a probléma, hogy az így adódó becslés „nem mond semmit”, azaz amennyiben nullát kapunk a mérés során, akkor a c értékére −∞ adódik, míg ha a megfigyelt érték 1, akkor pedig végtelen. Sőt, már \frac{K}{N} > \frac{1}{2} esetén is c negatív lesz, amit nem tudunk értelmezni. Ennek a becslésnek ilyen formában még várható értéke sincsen, azaz eleve nem is lehet konzisztens. Ezért a becslést a következőkben módosítani fogjuk. Olyan c értékekket szeretnénk kapni tehát, amelyekre z ≥ c ≥ 0 és a megfigyelt minta valószínűsége a legnagyobb. A fenti feltételből a

0 ≥ σΦ^{−1}(\frac{K}{N}) ≥ −z,

azaz a

0 ≥ Φ^{−1}(\frac{K}{N}) ≥ \frac{−z}{σ} ,

vagyis a

\frac{1}{2} ≥ \frac{K}{N} ≥ Φ(−\frac{z}{σ})

reláció teljesülése esetén közvetlenül értelmezni tudjuk az eredményt, azaz a c értékére adott becslésünk továbbra is −σΦ^{−1}(\frac{K}{N}), míg ha ez nem igaz, akkor \frac{1}{2} ≤ \frac{K}{N} esetében legyen c becslése értelemszerűen c = 0, illetve ha \frac{K}{N} ≤ Φ(−\frac{z}{σ}) teljesül, akkor pedig c = z. Itt tulajdonképpen egy momentumbecslésről van szó, amelyet a „szélén” módosítani kellett, ott a momentummódszer helyett maximum likelihood becslést használtunk. Az így kapott becslést jelölje f(X). Írjuk most fel a veszteségfüggvényt, amit minimalizálni akarunk. Ez sem lesz olyan egyszerűen átalakítható, mint az előző esetben, ugyanis most nem tudjuk, hogy torzítatlan becslésről volna szó. A veszteségfüggvény ebben az esetben is

E(f(X) − c))^{2} = E(f^{2}(X)) − 2cE(f(X)) + c^{2},

ahol ezt nem lesz az előző módszerrel könnyen kiszámítani, mert itt már nem független megfigyelések összegének szórásnégyzete szerepel. A fenti veszteségfüggvény értékének kiszámításához meg kell határozni f(X) eloszlását. Mivel ez a leírt becslési módszerben csak \frac{K}{N} -en keresztül függ a mintától, azért a továbbiakban f(X) helyett egyszerűen csak f(\frac{K}{N})-et írunk. Tudjuk, hogy K eloszlása N-edrendű Φ(−\frac{c}{σ}) várható értékű binomiális eloszlás, így a \frac{K}{N} valószínűségi változó eloszlása ebből könnyen meghatározható. Nagy mintaelemszám esetén alkalmazhatjuk a centrális határeloszlástételt, melyből azt kapjuk, hogy

\frac{K − NΦ(−\frac{c}{σ})}{\sqrt{NΦ(−\frac{c}{σ})(1 − Φ(−\frac{c}{σ}))}} ∼ Z

, ahol Z ∼ N(0, 1), azaz kis átalakítással azt kapjuk, hogy

\frac{K}{N} = \begin{matrix}√ & \overline{Φ(−\frac{c}{σ})(1 − Φ(−\frac{c}{σ}))} \\ √ & \overline{N} \end{matrix}Z + Φ(−\frac{c}{σ})

A becslés, ill. négyzete várható értékének kiszámításához először az f(\frac{K}{N}) = 0 esemény valószínűségét számoljuk ki. Tudjuk, hogy

P (ĉ = 0) = P (\frac{1}{2} ≤ \frac{K}{N})

vagyis azt írhatjuk, hogy

P (ĉ = 0) = P (\frac{1}{2} ≤ \begin{matrix}√ & \overline{Φ(−\frac{c}{σ})(1 − Φ(−\frac{c}{σ}))} \\ √ & \overline{N} \end{matrix}Z + Φ(−\frac{c}{σ})) =

= P (\frac{\frac{1}{2} − Φ(−\frac{c}{σ})}{\sqrt{Φ(−\frac{c}{σ})(1 − Φ(−\frac{c}{σ}))}}\sqrt{N} ≤ Z) = Φ(\frac{Φ(−\frac{c}{σ}) − \frac{1}{2}}{\sqrt{Φ(−\frac{c}{σ})(1 − Φ(−\frac{c}{σ}))}}\sqrt{N}).

Ha \frac{K}{N} ≤ Φ(−\frac{z}{σ}) áll fenn, akkor tudjuk, hogy a becslés z, ennek valószínűsége az előzőek szerint

P (Φ(−\frac{z}{σ}) ≥ \begin{matrix}√ & \overline{Φ(−\frac{c}{σ})(1 − Φ(−\frac{c}{σ}))} \\ √ & \overline{N} \end{matrix}Z + Φ(−\frac{c}{σ})),

azaz

P (\frac{Φ(−\frac{z}{σ}) − Φ(−\frac{c}{σ})}{\sqrt{Φ(−\frac{c}{σ})(1 − Φ(−\frac{c}{σ}))}}\sqrt{N} ≥ Z) = Φ(\frac{Φ(−\frac{z}{σ}) − Φ(−\frac{c}{σ})}{\sqrt{Φ(−\frac{c}{σ})(1 − Φ(−\frac{c}{σ}))}}\sqrt{N}).

A formulák áttekinthetősége érdekében vezessük be a következő jelöléseket.

a := Φ(−\frac{z}{σ}) d := Φ(−\frac{c}{σ})

P (ĉ = z) = P_{z} := Φ(\frac{a − d}{\sqrt{d(1 − d)}}\sqrt{N}) P (ĉ = 0) = P_{0} := Φ(\frac{d − \frac{1}{2}}{\sqrt{d(1 − d)}}\sqrt{N})

Az előzőek figyelembevételével

E(f(\frac{K}{N})) = −\frac{σ}{\sqrt{2π}} ∫a^{\frac{1}{2}} exp(−\frac{(x − d)^{2}N}{2d(d − 1)} )Φ^{−1}(x)dx(1 − P_{0} − P_{z}) + P_{z}z.

A veszteségfüggvényben szereplő másik nem konstans tag kiszámításához vegyük figyelembe, hogy a becslés négyzete, amelyet most meg kell határozni

[σΦ^{−1}(\frac{K}{N})]^{2} ha \frac{1}{2} ≥ \frac{K}{N} ≥ Φ(−\frac{z}{σ})

z^{2} ha \frac{K}{N} ≤ Φ(−\frac{z}{σ})

0 egyébként.

Ezek alapján felírhatjuk a keresett második tagot is, ahol használjuk az előzőekben bevezetett egyszerűsítő jelöléseket.

E(f^{2}(\frac{K}{N})) = \frac{σ^{2}}{\sqrt{2π}} ∫a^{\frac{1}{2}} exp(−\frac{(x − d)^{2}N}{2d(d − 1)} )(Φ^{−1}(x))^{2}dx(1 − P_{z} − P_{0}) + P_{z}z^{2}.

Mindkét tagban ugyanis feltételes várható értékről van szó, mégpedig az a feltétel, hogy olyan értéket kapunk-e, hogy az eloszlásfüggvény segítségével sztenderd módon tudunk becsülni. Ha ez az eset áll fenn, akkor ennek az eseménynek a valószínűségét szorozzuk az ilyen módon kapott becslés várható értékével, ha pedig nem teljesül a feltétel, akkor még két eset lehetséges. Ha a becslés nulla, azt a számolás során nem kell figyelembe venni, ha pedig z, akkor ennek az esemánynek a valószínűségét kell szorozni z-vel. Ezt azért tudjuk ilyen egyszerűen megtenni, mert pozitív a feltételek valószínűsége. A minimalizálandó kifejezés tehát összefoglalva a következő

2d\frac{σ}{\sqrt{2π}} ∫a^{\frac{1}{2}} exp(−\frac{(x − d)^{2}N}{2d(d − 1)} )Φ^{−1}(x)dx(1 − P_{0} − P_{z}) − P_{z}z+

+\frac{σ^{2}}{\sqrt{2π}} ∫a^{\frac{1}{2}} exp(−\frac{(x − d)^{2}N}{2d(d − 1)} )(Φ^{−1}(x))^{2}dx(1 − P_{z} − P_{0}) + P_{z}z^{2}

Analitikus módon nehéz volna ezt minimalizálni, ezért a feladatot ettől a ponttól numerikusan kezeljük. A MATLAB programot hívjuk segítségül, és a becsles2.m m-fájlban megírt függvény minimumát keressük úgy, hogy c értékét rögzítjük és σ-ban minimalizálunk. A néhány rögzített c értékhez tartozó optimális σ értékekből álló grafikont a mellékletben találhatjuk.

Megjegyzés: Az eredményt látva joggal merül fel a kérdés, hogy nagyon kis értékek esetén nem lesz-e kisebb a veszteségfüggvény értéke, ha egyáltalán nincsen zaj, tehát nem volna-e jobb, ha nem adnánk hozzá semmit. A veszteségfüggvény értéke zajmentes esetben valóban kisebb, csakhogy mivel a rendszerben egy alapzaj mindig jelen van, akkor azt érdemes az optimális értékre beállítani, így a kis tényleges értékek esetében is hasznos a zaj hozzáadása. A másik lehetséges észrevétel, hogy az ábra nem teljesen „szabályos”, azaz a grafikon a szélén megtörik. Ez a numerikus számítások miatt lehet, mert ekkor a minimalizáló eljárásban már igen kis számok jelennek meg.

2.1.2. Egy másik sztochasztikus rezonanciával kapcsolatos modell vizsgálata

Tekintsünk egy determinisztikusan, adott alakú differenciálegyenlettel leírható jelenséget, amelyben az ismeretlen megoldásfüggvény (ami speciálisan koncentráció, belső energia vagy akár egyedszám is lehet) monoton növő. Bizonyos (nyilván ismert) időpontokban méréseket végzünk, ahonnan csak annyit tudunk meg, hogy az ismeretlen függvény meghaladt-e már egy kritikus értéket. Természetesen azt tudjuk, hogy mennyi ez a kritikus érték. Ilyen jelenség fordulhat elő például ingerületátvitelnél, mikor mindössze az inerület létrejöttét tudjuk meg ill. bizonyos kémiai reakcióknál, ahol egy újabb reakció beindításához egy kritikus tömegre van szükség. Olyan rendszerek leírásakor is szükség lehet az alábbiakra, amikor a modell csak egy bizonyos telítődési határig írja le jól a folyamatot, mint például a Lotka-Volterra modellek egyes változatai, melyek a terület eltartóképességét nem veszik figyelembe. A fenti feltételeket matematikailag a következőképpen formalizáhatjuk. Adottak a t_{1}, t_{2}, …, t_{k} időpontok, amelyekben a méréseket végezzük. Adott továbbá egy M kritikus érték, és a folyamatot determinisztikus módon leíró ẋ(t) = f(x(t)) autonóm differenciálegyenlet. Tudjuk továbbá a mérések alapján, hogy az x(t_{1}), x(t_{2}), … x(t_{k−1}) értékek mind kisebbek, mint a kritikus érték, de a k-adik mérés során már azt meghaladó mennyiséget mértünk, tehát x(t_{k}) ≥ M. A fentiek alapján keressük az x(0) kezdeti értéket. A modell idáig teljesen determinisztikus, a kezdeti értékre vonatkozóan a differenciálegyenlet analitikus vagy numerikus megoldása alapján bizonyos határokat tudunk mondani, de amint azt egy konkrét példán keresztül látni fogjuk, az így kapott becslésnek túl nagy lesz a hibája, vagy legalábbis csökkenthető lesz. Az előző részhez hasonlóan azt reméljük, hogy ha a rendszerhez megfelelő zajt adunk, akkor a kezdeti értékre adott becslés jobb lesz. Természetesen ez csak úgy lehetséges, ha a rögzített időpontokban többször, egymástól függetlenül is mérhetünk. Jelöljük a lehetséges mérések számát n-el. A mérések nyilván úgy is értelmezhetők, hogy csak azt a t_{k} időpontot tudjuk meg, amikor a függvény a kritikus értéket biztosan elérte. Ezek az időpontok a keresett kezdeti feltételre nézve elégséges statisztikát adnak (mivel a mérések függetlensége miatt a mért értékek sorrendje „nem számít”). Megjegyezzük, hogy ez a modell több szempontból is különbözik az előzőtől. Egyrészt a mostani esetben nem egy konkrét, rögzített értéket kell mérni, hanem egy folyamatot, másrészt viszont most nem csak egy indikátorokból álló mintánk van, hanem a mérés során több lehetséges (persze előre meghatározott és diszkrét) időpontot is kaphatunk eredményként. Zajon most egy eleve egy minden időpillanatban független N(0, σ^{2}) eloszlású valószínűségi változókbó álló fehérzajfolyamatot értünk. A zaj hozzáadásával módosított modell a következő: A kezdeti érték egy X(0) = X ∼ N(τ, σ^{2}) valószínűségi változó. Ismerjük azt a k értéket, illetve több független mérésből többet is, amelyre

X(t_{k−1}) < M ≤ X(t_{k}).

A feladat τ becslése oly módon, hogy az előző részhez hasonlóan E(τ − τ̂)^{2} a lehető legkisebb legyen, ill. olyan σ érték meghatározása, amely esetén ilyen értelemben a lehető legjobb becslés adható. Mindezt egy konkrét példán szemléltetjük, ebben az esetben fogjuk a keresett σ értéket meghatározni, de a leírt módszer alkalmas arra, hogy más konkrét esetekben is kiszámoljuk a keresett értéket. Legyen most tehát a differenciálegyenlet a legegyszerűbb ẋ(t) = x(t), amely exponenciális növekedést fejez ki. Ha a tényleges kezdeti feltétel e^{τ_{0}}, akkor a megoldásfüggvény x(t) = e^{τ_{0}+t}, vagyis a mérések alapján azt tudjuk, hogy

e^{τ_{0}+t_{k−1}} < M ≤ e^{τ_{0}+t_{k}},

azaz

log M − t_{k} ≤ τ_{0} < log M − t_{k−1}.

De mivel a sztochasztikus modell értelmében τ_{0} helyett igazából kezdeti feltételként a fent bevezetett X valószínűségi változó által felvett érték szolgál, azért azt írhatjuk, hogy

log M − t_{k} ≤ X < log M − t_{k−1},

ahol X ∼ N(τ_{0}, σ^{2}). A kezdeti feltétel meghatározása a mostani speciális esetben tehát e^{τ_{0}} becslését jelenti. Azaz a következő kifejezést kell minimalizálni (ez lesz ugyanis a veszteségfüggvény):

E(e^{τ_{0}} − ê^{τ})^{2}.

Most vizsgáljuk meg a determinisztikus esetet, tehát nézzük meg, mennyire pontatlan a becslés akkor, ha determinisztikus modellt használunk. Ekkor az előzőek szerint

e^{log M−t_{k}} ≤ e^{τ} < e^{log M−t_{k−1}},

ahol csak ezt a két határt ismerjük, köztük a keresett kezdeti érték mindenhol azonos valószínűséggel lehet, vagyis e^{τ} ebben az intervallumban egyenletes eloszlást követ. Ennek szórásnégyzete pedig a következő lesz

\frac{(−e^{log M−t_{k}} + e^{log M−t_{k−1}})^{2}}{12} .

Megfelelő szórásnégyzetű zaj hozzáadásától azt reméljük, hogy ez az érték csökkenni fog.

A sztochasztikus modell esetén is próbáljunk meg hasonló módon becsülni, ahogyan azt az előbb tettük. Amennyiben a mérés során egy t_{i} értéket figyeltünk meg, akkor legyen a keresett kezdeti érték becslése most is a két határ számtani közepe, amiből a felírt veszteségfüggvényre vonatkozó képletet folytatva kapjuk, hogy

E(e^{τ_{0}} − ê^{τ})^{2} = E[(e^{τ_{0}} − \frac{1}{2} ∑_{i=1}^{n} \frac{e^{log M−t_{i}} + e^{log M−t_{i−1}}}{n})^{2}].

Itt jogosan merül fel a kérdés, hogy τ értékére tesszük fel, hogy egyenletes eloszlású a megengedett intervallumban, vagy pedig ugyanezt e^{τ} -ról feltételezzük. A két eset ugyanis különbözik, mi az első feltevéssel fogunk élni. E szerint az egyes mérésekből adódó becslés maximum likelihood becslés, azonban a több mintaelemre vonatkozó fenti becslésről még az sem látszik, hogy megengedhető volna. Ha azonban már erre a becslésre is igaz, hogy az előző, determinisztikus modellből kapott becslésnél kisebb veszteségű, akkor mindenképp érdemes zajt adni a kezdeti feltételhez. A formulák egyszerűbbé tétele érdekében vezessük be a következő jelölést:

X_{i} = e^{τ_{0}} − \frac{1}{2}(e^{log M−t_{i}} + e^{log M−t_{i−1}}).

Mivel méréseink egymástól függetlenek voltak, azért az X_{i} valószínűségi változók is függetlenek lesznek. A fenti kifejezést pedig egyszerűen így írhatjuk:

E(e^{τ_{0}} − ê^{τ})^{2} = E[(\frac{1}{n} ∑_{i=1}^{n} X_{i})^{2}] = E \frac{1}{n^{2}} (∑_{i=1}^{n} X_{i}^{2} + ∑_{i≠j}^{n} X_{i}X_{j}) =

= \frac{1}{n}E(X_{i}^{2}) + (1 − \frac{1}{n})(E(X_{i}))^{2}.

Ennek kiszámításához meg kell határoznunk X_{i} eloszlását. Annak valószínűségét kell tehát kiszámítani, hogy t_{k} lesz a kritikus idő a mérés során. Ez a fentiek alapján azt jelenti, hogy

log M − t_{k} − τ_{0} < X < log M − t_{k−1} − τ_{0},

ahol X ∼ N(0, σ^{2}). Ebből a keresett valószínűségre a következő adódik:

P (t_{k} a kritikus idő) = P (log M − t_{k} − τ_{0} < X < log M − t_{k−1} − τ_{0}) =

= \frac{1}{σ\sqrt{2π}} ∫_{log M−t_{k}−τ_{0}}^{log M−t_{k−1}−τ_{0}} exp(−x^{2}/2σ^{2})dx.

Mindezt a veszteségfüggvényt megadó képletbe helyettesítve azt kapjuk, hogy a veszteségfüggvény a következő módon számítható ki:

\frac{1}{n} ∑_{k=1}^{∞}[e^{τ_{0}} − \frac{1}{2}(e^{log M−t_{k}} + e^{log M−t_{k−1}})]^{2} \frac{1}{σ\sqrt{2π}} ∫_{log M−t_{k}−τ_{0}}^{log M−t_{k−1}−τ_{0}} e^{−x^{2}/2σ^{2}}dx+

+(1− \frac{1}{n})[∑_{k=1}^{∞} e^{τ_{0}} − \frac{1}{2}(e^{log M−t_{k}} +e^{log M−t_{k−1}})\frac{1}{σ\sqrt{2π}} ∫_{log M−t_{k}−τ_{0}}^{log M−t_{k−1}−τ_{0}} e^{−x^{2}/2σ^{2}}dx]^{2}.

A veszteségfüggvény minimális értékét megkaphatjuk tehát úgy, hogy a fenti kifejezést σ-ban minimalizáljuk. Így akkor egyúttal azt is megkapjuk, hogy melyik az a σ érték, amire az teljesül, hogy ha kezdetben olyan szórásnégyzetű fehérzajt adunk a jelenséghez, akkor az előzőekben előállított becslés a lehető legkisebb szórásnégyzetű lesz. Vegyük azonban észre, hogy még mindig nem egy konkrét, jól meghatározott képlet áll előttünk, hiszen előre nem ismerjük τ_{0} értékét, mert éppen azt keressük. Másrészt lehetséges (sőt, igaz is), hogy τ_{0} különböző értékeire a becslés veszteségfüggvénye különböző lesz. Így egy átlagos veszteségről kell beszélni úgy, hogy τ_{0}-ról a megfelelő intervallumon valamilyen eloszlást feltételezünk. A determinisztikus eset vizsgálatakor is ez történt, konkrétan egyéb τ_{0}-ra vonatkozó információ hiányában azt mondtuk, hogy e^{τ_{0}} eloszlása a megfelelő intervallumon egyenletes legyen. Ugyanezt fogjuk feltételezni a most tárgyalt sztochasztikus modell esetében is. A minimalizálást természetesen itt is numerikusan végezzük, ismét a MATLAB programcsomag segítségével. Az e^{τ_{0}} eloszlására vonatkozó feltevést úgy építjük be a számításba, hogy a megfelelő intervallumon egyenletes eloszlás szerint felvett pontokra számítjuk ki az integrált, majd a végén a kapott értékek átlagát vesszük. A mellékelt ábrán látható vízszintes vonal a determinisztikus esetben kapott veszteségfüggvény értéke, a másik grafikon pedig a különböző szórásnégyzetü zajok hozzáadása után (a leírt módszerrel kiszámítva) adódó veszteséget mutatja. A program futása során kapott értékek azonos σ esetén is különböznek, melynek oka az e^{τ_{0}} kezdeti feltételek egyenletes eloszlásának szimulációja. Ezért néhány eredmény átlagával dolgoztunk. További részletek a programhoz írt kommentárban találhatók.

Összefoglalás

A dolgozatban leírt eredmények, eljárások sok új kérdést vetnek fel. A legfontosabb, gyakorlati alkalmazást is érintő probléma talán az, hogy miként lehet a második fejezetben leírtakat egy rendszer optimális szabályozására, ill. pontosabb mérések megvalósítására felhasználni. Foglalkozni kellene a konstans érték becslésén kívül tetszőleges folyamat becslésével is. Érdemes volna megvizsgálni annak lehetőségét, hogy mesterségesen valamilyen periodikus függvény vagy egy idősor egyes értékeit adjuk a folyamathoz. A becslések statisztikai tulajdonságainak vizsgálatát is pontosítani kellene, a lehetséges torzítás csökkentésére, ill. az optimális becslés megtalálására kellene törekedni. A zaj hozzáadását nem csak a determinisztikus, hanem a sztochasztikus modellhez is illeszteni kellene. Az első rész utolsó fejezetében említett programok működését is hatékonnyá kellene tenni, vagy rámutatni arra, hogy elvileg sem lehet sokkal jobbakat írni, mert az eloszlások valamilyen értelemben túl közel vannak egymáshoz, azaz például a megkülönböztethetőségüket mérő információs divergenciájuk nagyon kicsi.

Irodalomjegyzék

[1] Arnold, L., „On the consistency of the mathematical models of chemical reactios”, in: Dynamics of Synergetic Systems Ed. H. Haken (Springer Verlag, Berlin-Heidelberg-New York, 1980), 107–118.

[2] Borovkov, A.A.: Matematikai statisztika (TYPOTEX, Budapest, 1999.), 153–155.

[3] Gardiner, C.W.: Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences (Springer-Verlag, Berlin-Heidelberg-New York-Tokyo, 1983), 237–244.

[4] Jetschke, G.: Mathematik der Selbstorganisation (Vieweg, Braunschweig-Wiesbaden, 1983.), 246–249.

[5] Karlin, S. -Taylor, H.M.: Sztochasztikus folyamatok (Műszaki Könyvkiadó, Budapest, 1984.), 144–145.

[6] Lowe, G. - Gold, G.H.: „Olfactory transduction is intrinsically noisy” Proc. Natl. Acad. Sci. USA Vol.92,pp. 7864–7868, August 1995.

[7] Tóth J. - Érdi P.: A sztochasztikus kinetikai modell nélkülözhetetlensége, In: Nemlineáris dinamika és egzotikus kinetikai jelenségek kémiai rendszerekben (Bazsa Gy. szerk.), Debrecen-Budapest-Gödöllő, 1992. pp. 117–143.