Näytetään tekstit, joissa on tunniste mallinnus. Näytä kaikki tekstit
Näytetään tekstit, joissa on tunniste mallinnus. Näytä kaikki tekstit

maanantai 12. helmikuuta 2024

Ekologinen kaavakokoelma

Tähän postaukseen olen pyrkinyt kokoamaan lyhyesti ekologisen sovellusalan kaavoja aiemmista postauksistani. Pidän tarkoituksella turhat selittelyt mahdollisimman lyhyenä, ja pyrin keskittymään vain olennaiseen. Lisätietoa malleista ja niitä koskevista ominaisuuksista löytyy aihetta tarkastelevista blogipostauksista. Poikkeuksena kaksi viimeistä mallia, joista en ainakaan kirjoitushetkellä ole julkaissut erillisiä kirjoitelmia.

Kasvumallit

Kahdesta vakiona käytetystä kasvumallista eksponentiaalinen kasvaa kiihtyvää tahtia ja sopiikin siksi ekologisia sovelluksia paremmin esimerkiksi koronkorkoa tai taloudellista kasvua kuvaaviin malleihin. Logistinen kasvu taas olettaa kasvun tasaantuvan, ja mallin on siksi havaittu sopivan edellistä paremmin myös ekologisiin ja väestötieteellisiin sovelluksiin.

Eksponentiaalinen kasvu

\begin{equation}\frac{dN}{dt}=(b-d)N(t)  \end{equation}

Logistinen kasvu

\begin{equation}\frac{dN}{dt}=rN(1-\frac{N}{K}). \end{equation}

Peto-saalismallit

Peto-saalis-mallit kuvaavat petoeläinten sekä niiden saaliiden välistä dynamiikkaa. Lotka-Volterra on alan klassikko ja yksi vanhimmista lajissaan, mutta mallin epämääräiset oletukset mm. populaatioiden kasvun ja petojen loputtoman ruokahalun suhteen johtivat nopeasti vaihtoehtoisten mallien kehittämiseen. Kahdesta suositumpi, Arditi-Ginzburg, painottaa petojen ja saaliseläinten välistä suhdetta dynamiikassaan, perusajatuksena siis, että liian vähä saalistiheys johtaa petojen harvenemiseen, mutta myös liian suuri petojen tiheys saaliiden vähenemiseen. Rosenzweig-MacArthur ottaa huomioon pelkän saalistuksen ohella myös eläinten luonnollisen poistuman. Molemmat mallit ovat erittäin herkkiä lähtöoletuksille sekä käytetyille funktioille, joten näiden kanssa on oltava tarkkoja.

 Lotka-volterra

\begin{equation} \frac{dx}{dt} = \alpha x - \beta xy \end{equation}
\begin{equation} \frac{dy}{dt} = \delta xy - \gamma y, \end{equation}

Arditi-Ginzburg

\begin{equation} \frac{dx}{dt}=f(x)x-g(\frac{x}{y})y \end{equation}
\begin{equation} \frac{dy}{dt}=\alpha g(\frac{x}{y})y -\beta y \end{equation}

Rosenzweig-MacArthur

\begin{equation} \frac{dx}{dt}=rx(1-\frac{x}{K})-yh(x) \end{equation}
\begin{equation} \frac{dy}{dt}=y(-\delta+\gamma h(x)) \end{equation}

Lajien välinen kilpailu

Lajien välisessä kilpailussa yksi tai useampi samasta elintilasta ja ekologisesta lokerosta kilpaileva laji on mallinnettu logistisen kasvumallin mukaan. Tässä tapauksessa oletuksena siis on, että elintilan resurssit kykenevät elättämään vain tietyn määrän eläimiä. Kahden lajin kesken voi syntyä myös tasapaino, jolloin molemmat kykenevät ylläpitämään elinkelpoista populaatiota.

Logistinen kilpailumalli

\begin{equation} \frac{dN_i}{dt} = r_i N_i (\frac{K_i-N_i}{K_i}) \end{equation}

Evoluutiomallit

Evoluutiomalleissa mallinnetaan valittujen alleelien siirtymistä sukupolvelta toiselle. Wright-Fisher mallissa oletetaan, että siirtymä sukupolvien välillä tapahtuu kerralla ilman hallitsemattomasti sekaisin lisääntyviä yksilöitä. Wright-Fisherin binomilaskentaan perustuva malli perustuukin eläimiä paremmin kausittain lisääntyvän ja lakastuvan kasvikunnan ilmiöiden kuvaamiseen. Eläinkunnan mallintamiseen soveltuvassa Moran-mallissa Markov-säännön mukaisesti yksilöt voivat vapaasti siirtyä tilasta toiseen, eli lisääntyä tai kuolla, ja eri sukupolvet voivat elää ja lisääntyä sekaisin keskenään.

Wright-Fisher -malli

\begin{equation}\binom{2N}{k}p^q q^{2N-k}=\frac{(2N)!}{k!(2N-k)!}p^q q^{2N-k}\end{equation}

Moran-malli

\begin{equation}P_{i,i-1}=\frac{N-1}{N}\frac{i}{N} \end{equation}
\begin{equation}P_{i,i}=1-P_{i,i-1}-P_{i,i+1}\end{equation}
\begin{equation}P_{i,i+1} = \frac{1}{N} \frac{N-1}{N} \end{equation}

lauantai 10. helmikuuta 2024

Peliteoria ja ekologia, osa 2: tasapainotiloja ja ESS

Sekastrategia-Nash

Tähän sopii lähteeksi tämä. Palataan aiempaan vangin dilemmaan.

Oletetaan, että puhtaan valinnan sijaan puhutaankin tunnetuista todennäköisyyksistä. Oletetaan esimerkiksi, että molempien pelaajien strategioiden todennäköisyydet ovat p,1-p,q ja 1-q, jolloin strategioiden mukaisiksi odotusarvoiksi saadaan
\begin{equation}2*p+3*(1-p)-3*p-(1-p)=0,\end{equation}

jolloin siis p = 2/3, sekä
\begin{equation}q*2+(1-q)*3-q*3-(1-q)=0,\end{equation}

jolloin q = 2/3
Näillä todennäköisyyksillä molempien odotusarvot ovat n. 2.3. Piirretään näiden perusteella kuvaaja

josta nähdään molemmat puhtaat sekä sekastrategia-Nashit. Toisin sanoen, kun q=p=0 ja q=p=1, on kyseessä Nashin tasapaino, sekä tässä tapauksessa sekastrategia-Nash kun q=p=2/3 ja q=p=1/3. Kummassakin tapauksessa tasapainopisteitä ovat siis (2,2) ja (1,1).

Korreloiva sekä bayesilainen Nash

Pätevä artikkeli löytyy Wikipediasta. Tähän kaveriksi Rogersin Game Theory: Analysis of Conflict. Bayesiläiseen peliin erikseen juttua täällä.   

Korreloivassa Nash-tasapainossa oletetaan pelaajien lisäksi jokin kolmas taho, kuten todennäköisyydet laskeva tietokone tai asiantuntija, joka antaa pelaajille suosituksen ja luvut valintojensa suhteen. Jos pelaajat päättävät pitäytyä annetuissa suosituksissa poikkeamatta niistä, on kyseessä korreloiva Nashin tasapaino. Kyseessä on vähän kuin edelleen jatkokehitelty versio aiemmasta sekastrategiamallista.

Bayesilaisen pelin Nash-tasapaino jatkaa edelleen saman todennäköisyysajattelun kehittelyä. Tässä versiossa pelaajat eivät tiedä toistensa voitoista. Tästä voisi toimia esimerkkinä vaikkapa suljettu huutokauppa, jossa ainoastaan tarjoajat tietävät oman hintansa. Tässä tapauksessa pelaajat pyrkivät jälleen maksimoimaan odotettavissa olevan hyödyn valinnoillaan, ottaen samalla huomioon hyötyyn valintoja koskevat todennäköisyydet. Valittujen maksimoitujen hyötyjen mukainen strategia vastaa bayesilaista Nashin tasapainoa.

Koska kaikki todennäköisyyspainotteiset tasapainotilat ovat niin lähellä toisiaan, en viitsi edes keksiä niille ylimääräisiä esimerkkejä.

Stackelberg-tasapaino

Aiemmissa tasapainotiloissa pelaajat tekivät päätöksensä samanaikaisesti tietämättä toistensa siirroista. Stackelberg-tasapaino perustuu kahden hengen pelattavaan Markov-ketjuna pelattavaan vuoropohjaiseen peliin, jossa pelaajat tietävät toistensa edellisen siirron päätöstä tehdessään. Tällön tasapainopiste on kohdassa, jossa molemmat pelaajat tekevät optimisiirtonsa. Tyypillisesti tällöin ensimmäisenä siirtonsa tekevä pelaaja saa seuraajaa paremman tuloksen.

Seuraava esimerkki on ensimmäisestä osasta:

Stackelbergin kohdalla pelaajat eivät siis tee samanaikaista umpimähkäistä päätöstä vaan ottavat vuoroja. Jos ensimmäinen pelaaja maksimoi hyötynsä, tämä valitsisi vaikenemisen. Havaittuaan ensimmäisen pelaajan valinnan toinen pelaaja maksimoisi oman hyötynsä vaikenemalla. Näin päädytään Stackelbergin tasapainoon, joka tällä kertaa on myös ESS ja Nashin tasapaino.

Käytännössä Stackeklberg-tasapainon idea on helposti sovellettavissa markkinoiden duopoli-tilanteeseen: tällöin samoista apajista kilpailee kaksi yritystä, jotka tekevät siirtoja esim. hinnoittelullaan.

keskiviikko 7. helmikuuta 2024

Peliteoria ja ekologia, osa 1: Nashin tasapaino, evoluutio ja ESS

ESS-lähde otettu täältä, minkä ohella ohella käytössä on ollut myös Smithin ja Pricen kirja. En päässyt käsiksi Thomasin vaihtoehtoista ESS-määritelmää käsittelevään aineistoon, joten tyydyin sen kohdalla tylsästi vain Wikipedian artikkeliin.

Vangin dilemma. Nashin tasapaino merkitty punaisella. Tässä tilanteessa kumpikaan osapuoli ei voi yksipuolisesti parantaa asemaansa muuttamalla strategiaansa, joten Nashin tasapaino on ruudussa (2,2).




Formaalimmin muotoiltuna: Olkoon kuvatunlainen nelikenttä sekä strategiat S ja T, jolloin siis diagonaali vasemmalta oikealle on (S,S) ja (T,T). Tällöin Nashin tasapaino saavutetaan tilanteessa, jossa
\begin{equation}E(S,S) \geq E(S,T),\end{equation}
eli kumpikaan osapuoli ei voi yksipuolisesti parantaa asemaansa muuttamalla strategiaansa. Nashin tasapaino ei välttämättä kuitenkaan ole pareto-optimaalinen, jolloin siis saavutettu tasapaino olisi kirjaimellinen optimitilanne kummallekin pelaajalle. Vangin dilemmassa ylempänä vangit pääsisivät yhteensä vähemmällä vaikenemalla.

Evolutionary Stable Strategy eli ESS voidaan kirjoittaa edellisen määritelmän tavoinS != T
\begin{equation}E(S,S) > E(S,T) \ \textbf{tai} \end{equation}
\begin{equation}E(S,S)=E(S,T) \ \textbf{ja} \ E(S,T)>E(T,T) \end{equation}

Määritelmän ylempi osa tarkoittaa "vahvaa Nashin ehtoa", eli tilannetta, jossa strategiaa muuttamalla ei voi yksipuolisesti parantaa tilannettaan, mutta jossa ei ole edes yhtä hyvää lopputulosta tarjolla, vaan pelkästään heikompia. Alemmassa oletetaan kolme lopputulemaltaan identtistä strategiaa sekä yksi muita huonompi. Jo määritelmästä siis huomataan, että Nashin tasapaino ei välttämättä aina ole ESS. Yllä kuvatussa vangin dilemmassa Nashin tasapaino on ESS, sillä kilpailevat strategiat (S,T) ja (T,S) ovat strategiaa (S,S) heikompia lopputulokseltaan. Tässä tapauksessa myös E(T,T) on ESS. E(T,T) on myös pareto-optimaalinen.

Tavallisen ESS-mallin ohella on myös olemassa "vahva" määritemä jonka mukaan kun S != T, pätee
\begin{equation}E(S,S) \geq E(S,T) \ \textbf{ja} \end{equation}
\begin{equation} E(S,T)>E(T,T). \end{equation}

Vahvan määritelmän mukaan ESS on Nashin tasapaino, jolle pätee lisäksi aiempi ESS-ehto. Vahva ESS ei välttämättä ole tavanomainen tai seuraavassa esitettävä heikko ESS. Tässä tapauksessa Nashin tasapaina jo tavallinen ESS eivät ole vahvasti ESS, sillä lopputulokseltaan paras strategia ei ole vaihtoehtoja heikompi. Vahvan määritelmän mukainen ESS voisikin siis löytyä tässä tapauksessa strategiasta (T,T), sillä tällöin se on kaikkia vaihtoehtoja parempi, mutta strategia (S,S) ei voi täyttää määritelmän alemman osan vaatimusta. Tällöin strategia (T,T) on myös pareto-optimaalinen.

Weakly Evolutionary Stable Strategy eli WESS voidaan kirjoittaa edellisen määritelmän tavoin, S != T
\begin{equation}E(S,S) > E(S,T) \ \textbf{tai} \end{equation}
\begin{equation}E(S,S)=E(S,T) \ \textbf{ja} \ E(S,T) \geq E(T,T) \end{equation}

Tällöin siis avautuu mahdollisuus, että kaikki strategiayhdistelmät ovatkin lopputulemaltaan yhtä hyviä. Vangin dilemman Nashin tasapaino ja ESS on määritelmän ylemmän osan perusteella myös heikko. Strategia (T,T) on paitsi heikko ESS, myös pareto-optimaalinen.

Jatketaan aiempia määritelmiä uusilla tuloksilla.

Jos strategialle S != T pätee
\begin{equation}E(S,T)>E(T,T), \end{equation}
on kyseessä on lokaalisti ylivertainen strategia. Tällöin kahden strategia S tuotta lähtökohtaisesti aina paremman lopputuloksen strategiaa T vastaan pelattuna.
Jos taas on olemassa sellainen ũ>0 ja u välillä (0,ũ), joilla pätee
\begin{equation}E(S, uT + (1-u)S)>E(T,uT+(1-u)S),\end{equation}
on strategialla S tasainen invaasiomuuri. Tällöin jos populaation alueella on pienempi kilpaileva populaati strategialla T, strategiaa S käyttävä populaatio hiljalleen syrjäyttää kilpailijansa luonnollisen kilpailuetunsa kautta. Oletuksena kaikki ryhmän jäsenet pitäytyvät samassa yksittäisessä strategiassa.

Seuraavat kolme määritelmää ovat yhtenevät:
1. Strategia on vahva ESS.
2. Strategialla on tasainen invaasiomuuri.
3. Strategia on lokaalisti ylivertainen.

Viimeiset kaksi määritelmää ja teoreema ovat Hofbauerin artikkelista, johon en pääse toistaiseksi käsiksi. Hofbauerin artikkelissa on tähän kuitenkin olemassa myös todistus, joten jos joskus satun sen löytämään, saatan sen käydä täällä tarkemmin läpi. Paremman puutteessa käytetään lähteenä tässä Stanfordin yliopiston filosofian laitoksen tieteellisiin julkaisuihin erikoistuneelta mainiolta SEP-sivustolta/tietokannasta löytyvää artikkelia. Strategian ylivertaisuus voidaan ainakin johtaa suoraan vahvan ESS:n määritelmästä.

Seuraavaan lähde täältä. En pääse lukemaan Sigmundin ja Nowakin alkuperäisartikkelia, joten käytän tässä vain sekundaarista lähdettä. Seuraavassa on sovellettu ESS-teoriaa populaatiokokoihin. Todellisessa maailmassa populaatioiden välinen dynamiikka nimittäin riippuu pelkän strategian ohella myös elintilasta ja resursseista kilpailevien yksilöiden määristä. Tätä tarkoitusta varten määritellään  myös erilliset kriittiset rajat populaatiokoille sekä näiden vaikutukselle ryhmien välisessä kilpailussa.

Olkoon evolutiivisia strategioita kuvaava matriisi muotoa.



Matriisin strategioiden E_A ja E_B kelpoisuus voidaan ilmaista muodossa
\begin{equation}E_A = ap_A +bp_A\end{equation}
\begin{equation}E_b = cp_B +dp_B.\end{equation}

Tällöin ehto mutaation A invaasiolle on

\begin{equation}E_A > E_B, \ eli \ p_A > \frac{d-b}{a-c}p_B. \end{equation}

Edellisestä huomataan, että strategioiden välinen siirtymäpiste on ehto

\begin{equation}p_A = \frac{d-b}{a-c}p_B\end{equation}

Oletetaan seuraavaksi, että a>c, tällöin voidaan sama kirjoittaa
\begin{equation}\lambda = \frac{p_A}{p_B}\end{equation}sekä invaasiokynnys  muodossa
\begin{equation}\lambda_0 = \frac{d-b}{a-c}.\end{equation}

Kun pätee
\begin{equation}\lambda_0 = \lambda, \end{equation}
kaksi strategiaa ovat tasapainossa, mutta kun siirrytään epäyhtälöihin, valloittaa toinen strategia elintilaa toiselta. Tällöin mutaatiot ja kilpailevat populaatiot syrjäyttävät vähitellen alkuperäispopulaatiota - tai toisin päin. Invaasiokynnyksen ollessa pieni vähäisemmänkin ulkopuolisen ryhmän on helppoa tunkeuta vallitsevan populaation elinalueille, ja homogeenisen populaation ylläpitäminen käy vaikeammaksi. Yllä annettu esitys siis perustuu sellaisenaan viimeksi annettuun lähteeseen.

Aiemmasta peto-saalismallista ja sen yhteydestä ESS-malliin täällä. En nyt ala turhia referoimaan erinomaisen tekstin sisältöä tänne, mutta sanotaan nyt vaikka, että dynaamisessa tapauksessa aiemmat stabiliteettitulokset sekä erityisesti asymptoottista stabiiliutta koskeva teoria osuvat näppärästi yhteen ESS-strategian kanssa. Järjestelmän ollessa stabiili populaatio kykenee siis pitämään puolensa mutaatioita vastaan.

Peliteorian dynaamisempa puolta tullaan käsittelemään erikseen ESS-sarjan myöhemmissä osissa.

lauantai 3. helmikuuta 2024

Stokastisen analyysin ja martingaalien lyhyt oppimäärä, osa1

Tämä blogipostaus perustuu pitkälti Gregory F. Fawlerin erinomaiseen luentomonisteeseen. Tässä postauksessa on lyhyesti referoitu keskeisimpiä osia ja jätetty monia tärkeitä yksityiskohtia erikseen mainitsematta ja todistamatta. Suosittelenkin lämpimästi lähdeaineiston lukemista, jos aihe vain vähääkään kiinnostaa. 

Ehdollinen todennäköisyys ja odotusarvo

Ehdollinen todennäköisyys kirjoitetaan
\begin{equation}P(Y|X),\end{equation}
eli tapahtuman Y kohdalla tiedetään X sekä kahden olevan keskenään toisistaan riippuvaisia. Tällainen voisi olla esimerkiksi perinnöllinen sairaus perheessä: tällöin riippumattomuusoletus johtaisi yhden lapsen sairastumisesta muiden perheenjäsenten sairauden todennäköisyyden vähenemiseen, kun riippuvuusoletuksella tulos on päinvastainen, ja yhden sairastuminen lisää muidenkin riskiä sairastua, sillä nämä jakavat keskenään vanhempiensa geneettisen perimän.

Ehdollinen odotusarvo voidaan kirjoittaa muodossa
\begin{equation}E(Y|F_n),\end{equation}
missä F_n eli filtraatio viittaa tapahtumaa Y edeltävään informaatioon.

Kävin tämän aiheen perusteellisemmin läpi gradussani, joten en mene tässä aivan perustasoa pidemmälle.

Martingaalit

Martingaalien teoria perustuu ehdollisen odotusarvon perustalle. Tällöin oletuksella n < t, martingaalille pätee
\begin{equation}E(M_t|F_n)=M_n.\end{equation}
Toisin sanoen martingaalin M_t odotusarvo on sen ehtoa vastaavan informaation hetken martingaali M_n. Martingaalia voisi käytännössä havainnollistaa reilulla kolikonheitolla, jossa kummallakin pelaajalla on tasaväkiset mahdollisuudet voittaa.

Martingaali voidaan kirjoittaa myös muodossa

\begin{equation}E(M_n-M_m | F_m)=0,\end{equation}
eli riippumatta hetkestä tai historiasta, tämän hetken odotetun voiton suuruus on 0.

Kyseessä on supermartingaali jos

\begin{equation}E(M_t|F_n) \leq M_n\end{equation}
ja submartingaali jos

\begin{equation}E(M_t|F_n) \geq M_n.\end{equation}
Jos martingaalille pätee jompi kumpi annetuista määritelmistä, aiempaa kolikonheittoesimerkkiä soveltaen, peli suosii jompaa kumpaa pelaajista, eikä mahdollisuus voittaa siksi ole tasaväkinen.

Brownin liike, stokastinen integrointi

Brownin liike tai Wiener-prosessi on jatkuva satunnaisuutta kuvaava malli. Olkoon B_t=B(t) satunnainen arvo hetkellä t. Näiden satunnaisten arvojen välille oletetaan korrelaatiota sekä t -> B(t) oletetaan satunnaisfunktioksi. Satunnaisprosessia B(t) sanotaan yksiulotteiseksi Brownin liikkeeksi satunnaiskululla m ja varianssilla sigma^2, jos seuraavat ehdot pätevät:

1. B_0=0
2. Kun s<t, S_t-S_s on normaalijakauma varianssilla sigma² ja satunnaiskululla m.
3. Jos s<t, satunnaismuuttuja B_t-B_s on riippumaton muuttujasta B_r kun r <= s.
4. Funktio t -> B_t on muuttujan t jatkuva funktio todennäköisyydellä 1.

Jos lisäksi sigma² = 0 ja m = 1, on kyseessä standardi Brownin liike. Oletetaan seuraavaksi B_t standardiksi Brownin liikkeeksi. Tällöin Brownin liike voidaan kirjoittaa muodossa

\begin{equation}Y_t=\sigma B_t + mt,\end{equation}
missä Y_t on normaalisti jakautunut Brownin liike, jolla variaatio on sigma² ja satunnaiskulku m. Muuttujan B_t standardiudesta siis seuraa muuttujan Y_t yleinen normaalijakautuneisuus.

Stokastinen differentiaaliyhtälö voidaan kirjoittaa muodossa

\begin{equation}dX_t = m(t,X_t)+\sigma (t,X_t)B_t, \end{equation}
missä B_t on standardi Brownin liike. Stokastinen differentiaaliyhtälö kirjoitetaan siis samoin kuin aiemmin mainitti normaali Brownin liike Y_t, mutta muuttujan B_t tilalle voidaan valita jokin useasta eri tavoin jakautuneesta satunnaisprosessista. Annetun differentiaaliyhtälön ratkaisu on muotoa

\begin{equation}X_t = X_0 + \int_0^t m(s,X_s)ds + \int_0^t \sigma (s,X_s)dB_s.\end{equation}
Iton kaava, kun oletukset ovat yksinkertaisimmillaan:

\begin{equation}df(B_t)=f'(B_t)dB_t+\frac{1}{2}f''(B_t)dt.\end{equation}
Iton kaava, kun otetaan huomioon aika ja paikka:
\begin{equation}df(B_t)=\partial_x f(B_t)dB_t+(\partial_{t} f(B_t)+ \frac{1}{2}\partial_{xx}f''(B_t))dt.\end{equation}
Prosessia X_t sanotaan geometriseksi Brownin liikkeeksi, jos standardilla Brownin liikkeellä B_t, varianssilla sigma² sekä satunnaiskululla m pätee
\begin{equation}X_t=mX_t dt+ \sigma X_t dB_t = X_t(mdt+\sigma dB_t).\end{equation}
Samainen prosessi on keskeisellä sijalla rahoitusalan sovelluksissa ja Black-Scholesissa, joten tämän kanssa kannattaa olla tarkkana.

Radon-Nikodym, Girsanovin teoreema

Olkoon u,v sigma-äärellisiä avaruudessa (omega,F) ja v << u. Tällöin on olemassa funktio f, jolle pätee kaikilla E
\begin{equation}v(E)=\int_E du.\int\end{equation}
Funktiota f sanotaan mitan v Radon-Nikodym-derivaataksi suhteessa mittaan u. Funktiota merkitään

\begin{equation}f=\frac{dv}{du}.\end{equation}
Oletetaan seuraavaksi todennäköisyysavaruus (omega,F,P) ja todennäköisyysmitta Q, jolla pätee Q << P. Tällöin Radon-Nikodym-derivaatta

\begin{equation}X=\frac{dQ}{dP}\end{equation}
on ei-negatiivinen muuttuja, jonka odotusarvo on E(X)=1. Jos A on tapahtuma, voidaan kirjoittaa

\begin{equation}Q(A)=E_P(X1_A).\end{equation}
Todennäköisyysmittojen P ja Q odotusarvojen suhde voidaan vielä erikseen kirjoittaa muodossa
\begin{equation}E_Q(A)=E_P(X \frac{dQ}{dP})\end{equation}

Sovelletaan aiempia tuloksia pohjustamalla Girsanovin teoreemaa muutamilla aputuloksilla. Olkoon ei-negatiivinen martingaali M_t, jolle on voimassa

\begin{equation}dM_t = M_t A_t dB_t,\end{equation}
ja jonka ratkaisu on muotoa

\begin{equation}M_t = e^{Y_t}, Y_t = \int_0^t A_s dB_s - \frac{1}{2} \int_0^t (A_s)^2 ds.\end{equation}

Oletetaan seuraavaksi, että ratkaisu on martingaali, vaikka näin ei aina välttämättä olekaan. Määritellään sen avulla mitta P*: jos V on F_t-mitallinen tapahtuma, niin pätee

\begin{equation}P*=E(1_V M_t).\end{equation}
Jos oletamme, että mitat P ja P* on määritelty F_t-mitallisissa tapahtumissa, voidaan kirjoittaa

\begin{equation}\frac{dP*}{dP}=M_t,\end{equation}
mikä vastaa siis aiemmin määriteltyä Radon-Nikodym-derivaattaa. Seuraavaksi huomataan, että jos s < t ja V on F_s-mitallinen tapahtuma, niin V on myös F_t-mitallinen tapahtuma, jolloin siis pätee

\begin{equation}E(1_V M_s)=E(1_V M_t)\end{equation}
sekä F_t-mitallisella muuttujalla X

\begin{equation}E*(X)=E(X M_t).\end{equation}

Olkoon mitta P* kuten äskeisessä. Tällöin Girsanovin teoreeman mukaan yhtälössä

\begin{equation} dBt = A_t dt + dW_t\end{equation}
W_t on P*-mitallinen standardi Brownin liike, jolla on voimassa

\begin{equation}W_t = B_t - \int_0^t A_s ds.\end{equation}

Seuraava lyhyt esimerkki on peräisin täältä.

Finanssimatematiikassa Girsanovia käytetään riskineutraalin stokastisen satunnaiskulun, rahoitusintrumentin hinnan, määrittelemiseen. Oletetaan seuraavassa todennäköisyysavaruus riskineutraalilla todennäköisyysmitalla Q.


Kirjoitetaan ensin stokastinen differentiaaliyhtälö

\begin{equation}dY_t = (\mu-r)Y_t dt + \sigma Y_t d B_t, \end{equation}
ja muunnetaan se muotoon
\begin{equation}dY_t = \sigma \frac{(\mu-r)}{\sigma}Y_t dt + \sigma Y_t d B_t, \end{equation}
ja edelleen

\begin{equation}dY_t = \sigma Y_t d(\frac{(\mu-r)}{\sigma}t + B_t)=\sigma Y dW_t,\end{equation}
eli toisin sanoen pätee
\begin{equation}dW_t = B_t - \frac{(r-\mu)}{\sigma}t = B_t - \int_0^t A_s ds.\end{equation}
aivan Girsanovin teoreeman mukaisesti.

Euler-Maruyama ja stokastisten funktioiden visualisointi

Stokastisten differentiayhtälöiden ja satunnaiskulun visualisoinniksi on olemassa perinteisestä Eulerin rekursiomenetelmästä johdettu hyödyllinen stokastinen sovellus, joka kantaa muotoa
\begin{equation}Y_{n+1}=Y_n + a(Y_N,\tau_n)\triangle t + b(Y_n,\tau_n) \triangle W_n, \end{equation}
missä kolmio-t kuvaa aikaintervallia sekä kolmio-W normaalijakautuneen satunnaisluvun hetken t erotusta sitä edeltäneestä hetken t-1 satunnaisluvusta. Seuraavassa esimerkkikuvassa on simuloitu koron satunnaiskulkua kuvaavan stokastisen Vasicek-mallin eri realisaatioita:


Euler-Maruyama on välttämätön työkalu stokastisten yhtälöiden parissa työskennellessä, sillä menetelmä yksinkertaistaa ja selkeyttää käytännön sovellusten suhteen niin tehokkaasti muuten varsin vaikeaselkoista stokastista analyysiä. Minulla on tästä pitkät pätkät määritelmineen päivineen eräässä vanhassa harjoitustehtävässäni, joten en pistä täydellistä kuvausta aiheesta tänne.

Feynman-Kac, rahoituksen stokastiset menetelmät ym. lisäteoriaa jatketaan myöhemmin toisessa, aiheeseen syvällisemmin keskittyvässä blogipostauksessa.


Käytetty R-koodi:

T0 <- 1
r0 <- 0.03
kappa <- 0.35
thetaQ <- 0.0035/kappa
sigma <- 0.05

hedge.points <- 500
simulations.total <- 5
maturity <- T0

dt <- maturity/(hedge.points+1)
timeline <- seq(0,maturity, dt)

f <- matrix(r0,(hedge.points+2),simulations.total)

vasicek_rate <- function(r,kappa,theta,sigma,dt){
    expkappadt <- exp(-kappa*dt)
    vasi_vola <- (sigma^2)*(1-expkappadt^2)/(2*kappa)
    result <- r*expkappadt+theta*(1-expkappadt)+sqrt(vasi_vola)*rnorm(1)
    return(result)
}

for(i in 2:(hedge.points+2)){
    for(j in 1:simulations.total){
        f[i,j]<- vasicek_rate(f[i-1,j],kappa,thetaQ,sigma,dt)
    }
}

plot(timeline,f[,1], ylim=range(f,thetaQ), type="l", col=1,ylab="Korko (r)", xlab = "Hetki (t)")
for(j in 2:simulations.total){
 lines(timeline,f[,j], col=j )
}

title(main="Koron satunnaiskulun simulaatio Vasicek-mallin mukaan")

sunnuntai 28. tammikuuta 2024

Lyhyt johdatus matemaattiseen mallinnukseen osa 2: bifurkaatiot ja kaaosteoria

Isokliinit, virtaus- ja vaihekartat

Lotka-Volterran tasapainopisteitä ovat alfa/beta ja gamma/sigma. Piirretään näiden pohjalta isokliinit (alla), ja huomataan, että kun tasapainopiste 1:n yläpuolella populaatio A pienenee ja alapuolella kasvaa, populaatio B taas vähenee tasapainopiste 2:n vasemmalla ja kasvaa oikealla puolella. Kahden populaation dynamiikan yhteisvaikutusta on kuvattu mustalla nuolella:


Luciano da F. Costa
on piirtänyt Lotka-Volterran dynamiikasta virtaus- ja vaihekartan, josta huomataan virtaussuunnan seuraavan isokliinien avulla saatua vaikutelmaa. Costan mallissa kaikki parametrit ovat 1, jolloin myös kahden isokliinin leikkauspiste on (1,1) ja kyseessä on epästabiili tasapainopiste. Koska pisteessä ominaisarvot olivat imaginäärisiä ilman reaalisia osia, on järjestelmä vaimentamaton oskillaattori.


Bifurkaatiot

Bifurkaatioilla tarkoitetaan väljimmillään yllättävää muutosta funktion laadullisessa käytöksessä. Satula-solmu-bifurkaatiossa kaksi tasapainopistettä kohtaavat ja eliminoivat toisensa, transkriittisessä bifurkaatiossa tasapainopisteiden stabiilisuusominaisuudet vaihtuvat, ja kahdentumisbifurkaatiossa tasapainopisteiden määrä vaihtelee arvojen muutoksen mukaan. Näiden lisäksi on myös Hopfin bifurkaatiossa, jossa tasapainopisteiden stabiilisuus- sekä jaksollisuusominaisuudet vaihtuvat. Aiheesta löytyy lisää informaatiota sekä muodolliset todistukset Parvisen luentomonisteesta (kevät 2022).

Seuraavassa nähtävän logistisen kartan ja bifurkaation tarkemmat laskennalliset perusteet löytyvät täältä. Aiemmin käsiteltiin logistista kasvumallia

\begin{equation}\frac{dN}{dt}=F(N)=rN(1-\frac{N}{K}). \end{equation}

Kirjoitetaan tässä logistinen funktio rekursiomuodossa

\begin{equation}x_{t+1}=rx_t(1-x_t),\end{equation}

jolla on haarukkabifurkaatiot kun r= 3,3.449, 3.544... Hyvin yksinkertaisena ideana on laskea funktion rekursiokaavan avulla uusia ilmentymiä ajassa t, laskea kiintopisteet ja arvioida niiden stabiilius sekä käyttäytyminen eri kasvuvakio r:n arvoilla. Kuten aiemmassa linkissä huomataan, alkaa rekursiofunktio monimutkaistuessaan saada yhä useampia juuria, joiden vuoksi koko yhtälö alkaa hajota toistuvien bifurkaatioiden sekä kiintopistedynamiikan vuoksi yhä kaoottisemmaksi ja hallitsemattomammaksi sotkuksi.Seuraavassa nähtävän logistisen funktion dynamiikan bifurkaatioiden Feigenbaumin kuvaajan R-koodi löytyy täältä


Kiehtovinta logistisessa funktiossa on, että vaikka sillä onkin kaoottisia vaiheita ja näennäisen loputtomasti bifurkaatioita, löytyy r=3 ja r<4 väliltä myös vakaita, ei-kaoottisia ja ei-bifurkaatioita aiheuttavia poikkeamia. Kun r >= 4, populaatio romahtaa. Logistisen bifrukaatiokasauman kiehtovista poikkeamista, kaaoksen teoriasta ja säännönmukaisuuksista on myös erinomainen Wikipedia-artikkeli myöhempää käyttöä varten.

Kaaos

Arkikielessä kaaoksella tarkoitetaan löyhästi jonkinlaista hallitsematonta epäjärjestystä. Matematiikassa kaaosta koskeva teoria on huomattavasti "tylsempää" ja määritelmät suurpiirteisempiä. Määritelmistä puhutaankin tässä kohtaa monikossa, sillä niitä riittää aiheelleen uskollisesti pökerryttävän suuri määrä. Seuraavassa muutama esimerkki:

(Kaaos Lin ja Yorken mukaan): Olkoon J reaaliakselin väli. Kuvaus f : J → J
on kaoottinen Lin ja Yorken mielessä, jos se on jatkuva ja sillä on kolmejaksoinen rata. Toisinaan vaaditaan lisäksi herkkä riippuvuus alkuehdoista.

(Kaaos Devaneyn mukaan). Olkoon X metrinen avaruus. Jatkuva kuvaus
f : X → X on kaoottinen Devaneyn mielessä joukossa X, jos seuraavat kolme ehtoa ovat voimassa:

1. f on topologisesti transitiivinen
2. kuvauksen f jaksolliset pisteet ovat tiheässä joukossa X
3. kuvauksella f on herkkä riippuvuus alkutilasta

Yllä mainitut määritelmät löytyvät tarkempine lisätietoineen Parvisen luentomonisteesta (kevät 2022).

Näiden lisäksi Feldman antaa kirjassaan Chaos and Dynamical Systems vielä seuraavan löyhän määritelmän: dynaaminen systeemi on kaoottinen, jos

1. sen aikariippuvainen kulku on deterministinen funktio
2. sen radat ovat rajoitettuja
3. sen radat ovat herkkiä lähtöoletuksille
4. sen radat ovat epäsäännöllisiä

Erikseen on mainittava, että vaikka kolme eri määritelmää tietyssä mielessä risteävätkin, ei yhden määritelmän mukainen kaaos välttämättä ole sitä kahden muun mukaan.

Ehkä kuuluisin kaikista kaoottisista epälineaarisita dynaamisista järjestelmistä on Lorenzin systeemi
\begin{equation}\frac{dx}{dt}\sigma (y-x)\end{equation}
\begin{equation}\frac{dy}{dt}x(\rho-z)-y\end{equation}
\begin{equation}\frac{dz}{dt}xy-\beta z,\end{equation}

joka kuvaa tasaisen, alta lämpenevän ja yltä viilenevän nestekerroksen dynamiikkaa. Kyseessä on epälineaarinen, epäjaksollinen, kolmiulotteinen deterministinen systeemi, jolla on kokonaiset kolme tasapainopistettä. Kun p=1, on funktiossa haarukkabifurkaatio. Kun taas p=28, sigma=10 ja beta = 8/3, on systeemillä kaoottisia ratkaisuja. Tällöin systeemi jää kiertämään äärettömästi Lorentzin attraktorin, kirjallisuudesta riippuen myös "oudon attraktorin" ympärille, ja systeemi hakeutuu kaikkiin kolmeen tasapainopisteeseen. Tällöin myös yllä annetut määritelmät pätevät, ja ratojen kuvaajaksi saadaan Wikipedian koodilla


Kaaosteorian attraktorit ovat oma kiehtova ryhmänsä, josta voisi kirjoittaa vaikka kuinka paljon. Tässä kuitenkin vain koskettiin aiheeseen näin yleisellä tasolla, joten ei siitä sen enempää. Ehkä joskus myöhemmin kirjoittanen aiheesta lisää.

lauantai 20. tammikuuta 2024

SIS ja SIRVI-malli

SIS-malli

Lähde löytyy täältä.

SIS-malli kuvaa tilannetta, jossa taudin sairastaminen ei tuotakaan pysyvää immuniteettia, vaan henkilö palaakin takaisin lähtötilanteeseen eli taudille alttiiksi.
\begin{equation}I(t+\triangle t) = I(t)+\beta \triangle t S(t) I(t) - \gamma \triangle t I(t), \end{equation}
kun kolmio-t lähestyy nollaa, saadaan differentiaaliyhtälö

\begin{equation}\frac{dI}{dt} =\frac{\beta SI}{N} - \frac{\gamma I}{N}. \end{equation}

sekä vastineeksi

\begin{equation} \frac{dS}{dt}=-\frac{\beta SI}{N}+\frac{\gamma I}{N}.\end{equation}

Oletus I+S=N saadaan käännettyä muotoon S=N-I, sekä sijoittamalla saatu muoto edelliseen differentiaalityhtälöön saadaan

\begin{equation}\frac{dI}{dt} \beta (N-I)I - \gamma I. \end{equation}

sekä edelleen muotoilemalla

\begin{equation}\frac{dI}{dt} = (\beta N - \gamma)I(1-\frac{\beta I}{N-\gamma}), \end{equation}

joka on siis logistinen kasvumalli.

Edellinen differentiaaliyhtälö voidaan edelleen ratkaista. Ratkaisun pituuden ja monimutkaisuuden vuoksi en käy ratkaisumenetelmää yksityiskohtaisesti tässä läpi, täydellinen ratkaisu löytyy tästä Khanin, Sabbirin ja Sadiqin artikkelista. Differentiaaliyhtälön ratkaisu lähtee ylempänä muotoillun yhtälön integroinnista. Kun muuttujan I differentiaaliyhtälö on saatu ratkaistua, saadaan muuttujan S ratkaisu yksinkertaisesti käyttämällä identiteettiä S = N-I, jolloin ratkaisuiksi saadaan

\begin{equation}S=N - \frac{\beta N - \gamma}{\beta - (\beta N - \gamma)\frac{\beta N - \gamma - i_0 \beta}{(\beta N - \gamma)i_0}e^{-(\beta N - \gamma )t}}. \end{equation}

sekä muuttujan I yhtälöksi

\begin{equation}I=\frac{\beta N - \gamma}{\beta - (\beta N - \gamma)\frac{\beta N - \gamma - i_0 \beta}{(\beta N - \gamma)i_0}e^{-(\beta N - \gamma )t}}. \end{equation}


Skipataan tässä yksinkertaisemmassa mallissa tasapainopisteiden haku sekä niiden stabiiliusominaisuuksien analyysi.  Aiheesta löytyy erinomainen esitelmä täältä.
 
Tartuntoja kuvaavien differentiaaliyhtälöiden kuvaaja on


SIS-mallissa perustartuntaluku on muotoa
\begin{equation}R_0 = \frac{\beta N}{\gamma}, \end{equation}

eli jos R_0 > 1, epidemia käynnistyy. Tällöin sairastuneella on vähintään yksi kontakti, ja tauti voi siten lähteä leviämään.

SIRVI-malli 

SIRVI-malli lähtee pitkälti samoista oletuksista kuin kaksi edellistäkin, sillä erotuksella, ettei rokotus tällä kertaa annakaan pysyvää immuniteettiä. SIRVIn asetelma muistuttaa siis vahvasti koronarokotteen vaikutuksia, ja on siksi hyvin ajankohtainen. Koronarokote ei siis antanut kestävää immuniteettia tautia vastaan, estänyt ketään sairastamasta tai levittämästä tautia, mutta vahvisti hetkellisesti immuniteettia tautia vastaan. SIRVIin löytyy erinomaisen perusteellinen lähde täältä, ja käytän sitä myös tämän blogipostauksen pohjana. Kaikki olennainen tieto on peräisin Turkyilmazoglun artikkelista, enkä ota tästä kunniaa itselleni.

SIRVI-malli muodostuu yhtälöistä:

\begin{equation} \frac{dS}{dt}=-\beta SI -vS\end{equation}

\begin{equation} \frac{dI}{dt}=\beta SI + \alpha VI -\gamma I\end{equation}

\begin{equation} \frac{dR}{dt}=\gamma I \end{equation}

\begin{equation} \frac{dV}{dt}=v S - \alpha VI \end{equation}

missä alfa, beta, gamma ja v ovat rokotettujen tartuntamäärä, epäiltyjen tartuntamäärä, parantumismäärä ja epäiltyjen rokotusmäärä.

Hetkellä t=0 annettujen funktioiden muuttujille pätee:

\begin{equation} S_0 = 1-\tau \end{equation}

\begin{equation} I_0 = \tau \end{equation}

\begin{equation} R_0 = 0 \end{equation}

\begin{equation} V_0 = 0 \end{equation}.

Tällöin on voimassa myös yhtälö

\begin{equation}S+I+R+V = 1 \end{equation}.
Mallin perustartuntaluku on

\begin{equation}R_0=\frac{S_0 \beta -V_0 \alpha}{\gamma}\end{equation}
sekä tasapainopisteet

\begin{equation}(S,I)=(0,0)\end{equation}
\begin{equation}(V,R)=(0 \leq V, R \leq 1)\end{equation}

SIRVI-mallin Jacobin matiiris on muotoa

\begin{equation}J=\begin{bmatrix}-v & -SB & 0 & 0 \\ 0 & V \alpha + S \beta - \gamma & 0 & 0 \\ v & -V \alpha & 0 & 0 \\ 0 & \gamma & 0 & 0  \end{bmatrix}\end{equation}

SIRVI-mallin stabiilisuusominaisuuksista voidaan ainakin sanoa, että nollasta eroavat ominaisarvot

\begin{equation}\lambda_1 = -v \end{equation}
sekä

\begin{equation}\lambda_2 = V \alpha + S \beta - \gamma, \end{equation}
mistä huomataan, että stabiilisuus pätee rokotetulla ja taudista vapaalla populaatiolla kun

\begin{equation}V \leq \frac{\gamma}{\alpha}\end{equation}

Kuvaaja on siis muotoa



Vaikka tässä käsiteltiinkin vain kolmea erityistä epidemiologista mallia, riittää vastaavia sovellutuksia maailmassa vaikka kuinka paljon. Epidemiologia on yksi sovelletun matematiikan merkittävimmistä haaroista.


Käytetyt koodit, perusten aiempaan SIR-mallissa käytettyyn Python-koodiin:

# SIS-malli
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Total population, N.
N = 1000
# Initial number of infected and recovered individuals, I0 and R0.
I0 = 1
# Everyone else, S0, is susceptible to infection initially.
S0 = N - I0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma = 0.09, 2.0
# A grid of time points (in days)
t = np.linspace(0, 160, 160)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma):
    S, I = y
    dSdt = -beta * S * I / N + gamma * I/N
    dIdt = beta * S * I / N - gamma * I/N
    return dSdt, dIdt

# Initial conditions vector
y0 = S0, I0,
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I = ret.T

# Plot the data on three separate curves for S(t), I(t) and R(t)
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()

#################################################################

# SIRVI-malli
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Total population, N.
N = 1000
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0, V0 = 1, 0, 0
# Everyone else, S0, is susceptible to infection initially.
S0= N - I0 - R0 -V0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma, v ,alpha= 0.2, 2.10, 0.03, 0.03
# A grid of time points (in days)
t = np.linspace(0, 160, 160)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma, v, alpha):
    S, I, R, V = y
    dSdt = -beta * S * I/N - v * S
    dIdt = beta * S * I/N + alpha*V*I/N - gamma * I/N
    dRdt = gamma * I/N
    dVdt = v * S-alpha*V*I/N
    return dSdt, dIdt, dRdt, dVdt

# Initial conditions vector
y0 = S0, I0, R0, V0
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma,v,alpha))
S, I, R, V = ret.T

# Plot the data on three separate curves for S(t), I(t) and R(t)
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
ax.plot(t,V/1000, 'y', alpha=0.5, lw=2, label = 'vaccinated' )
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()

torstai 18. tammikuuta 2024

Vaihtoehtoja Lotka-Volterralle

Diclaimer: Bloggerin matematiikkatyökalu on järjettömän buginen ja epäluotettava. Siksi olenkin tekstissäni pyrkinyt mahdollisuuksien mukaan välttämään ylimääräisten symbolien käyttöä.

Edellisessä, Lotka-Volterran peto-saalis -mallia käsittelevässä postauksessa käytiin läpi mallin historiaa sekä yleisiä ominaispiirteitä ja rajoituksia. Huomattiin, että kyseinen malli perustuu hyvin epärealistisille oletuksille peto- ja saaliseläinten kantojen ominaisuuksista. Käytännössä onkin huomattu, ettei Lotka-Volterra aina pädekään erityisen tarkasti kaikissa tilanteissa, minä vuoksi mallille onkin vuosien varrella kehitetty vaihtoehtoisia, oletuksiltaan joustavampia peto-saalis -malleja. Tässä postauksessa käydään läpi kaksi keskeisintä ja ehkäpä yleisimmin käytettyä: Arditi Ginzburg sekä Rosenzweig-MacArthur.

Mistä Lotka-Volterraa siis tarkalleen ottaen kritisoitiin? Ensinnäkin: Lotka-Volterra lähti oletuksesta, että petoeläimillä on rajaton ruokahalu, mikä ei reaalimaailmassa voi pitää eikä pidäkään paikkanasa. Samoin Lotka-Volterrassa petoeläinten puuttumisen seurauksena saaliseläimen kanta kasvaa Malthusin lain mukaisesti rajattomasti - tämäkään oletus ei välttämättä kestä päivänvaloa. Tässä  erinomaisesa luennossa on käsitelty ekologisten mallien ominaisuuksia syvällisemmin ja kattavammin kuin tässä blogipostauksessa.

Arditi-Ginzburg

Arditi-Ginzburgin peto-saalis -malliksi kutsutaan yhtälöparia

\begin{equation} \frac{dx}{dt}=f(x)x-g(\frac{x}{y})y \end{equation}

\begin{equation} \frac{dy}{dt}=\alpha g(\frac{x}{y})y -\beta y \end{equation}
jossa funktio f() kuvaa saaliseläimen kannan luonnollisia muutoksia, funktio g() liittyy saaliseläinten sekä petojen kantojen suhteeseen, alfa on vakio petojen saaliiksi saamasta lisääntymisedusta sekä beta petojen kuolemaa kuvaava suhdeluku. Arditi-Ginzburg ottaa siis huomioon monipuolisemmin populaatioiden muutoksen muutenkin kuin näiden keskinäisen riippuvuuden kautta. Malli onkin vahvasti suhdelukupainotteinen. Arditin ja Ginzburgin päätyminen mainittuihin suhdelukuihin on pitkä tarina, josta voi lukea tarkemmin täällä, referoin keskeisiä kohtia tänne muutamin täydennyksin ja huomioin. 

Tässä kohtaa on syytä huomata, ettei funktioita f() ja g() oleteta yleisesti annetuiksi, vaan niiden paikalle voidaan sovittaa tarpeen mukaan erilaisia suhdeluvuilla toimivia troofisia funktioita eriävine ominaisuuksineen. Tässä blogipostauksessa on tarkoituksella käytetty mahdollisimman yksinkertaisia troofisia funktioita

\begin{equation}f(x)=\gamma x \end{equation}

sekä

\begin{equation}g(\frac{x}{y})= \frac {\gamma x}{y} \end{equation}

jolloin yhtälöpari voidaan kirjoittaa muotoon

\begin{equation} \frac{dx}{dt}=\gamma x² - \gamma x \end{equation} 

\begin{equation} \frac{dy}{dt}=\alpha \gamma x -\beta y \end{equation}

ja asettamalla derivaatat nolliksi, saadaan tasapainopisteiksi (0,0) , (1,gamma*alfa/beta), joiden muodostama Jacobin matriisi on muotoa

\begin{bmatrix} 2 \gamma x - \gamma & 0 \\ \alpha \gamma & -\beta \end{bmatrix}

mistä huomataan, että mikäli beta ja gamma oletetaan aina nollaa suuremmaksi, ovat matriisin ominaisarvot origossa aina samanmerkkiset, eli järjestelmä on annetuilla arvoilla ja oletuksilla lähtökohtaisesti epävakaa. Toisen tasapainopisteen kohdalla itseisarvoista on toinen positiivinen ja toinen negatiivinen, jolloin kyseessä on epävakaa satulapiste. Ominaisarvoista ja stabiiliudesta lisää tämän linkin takaa.

Demonstroisin mieluusti mallia soveltuvalla kuvaajalla, mutta annetuilla oletuksilla ja valinnoilla kannat nousevat aluksi vain pudotakseen nopeasti nollaan, joten vaikka esimerkkinä käytetyt funktiot malliin sopivatkin, eivät ne yhdessä näytä tuottavan juurikaan kelvollista jälkeä. Tässä yksi pätevä opetus siitä, miksi lähtöoletusten kanssa on oltava tarkkana!

Rosenzweig-MacArthur

Rosenzweig-MacArthur -malliksi kutsutaan yhtälöparia

\begin{equation} \frac{dx}{dt}=rx(1-\frac{x}{K})-yh(x) \end{equation}

\begin{equation} \frac{dy}{dt}=y(-\delta+\gamma h(x)) \end{equation}

 Rosenzweig-MacArthur -mallista voi lukea lisää täällä. Yhtenä keskeisenä tavoitteena Rosenzweigilla ja MacArthurilla oli saaliseläimen populaation malthuslaisen kasvun rajoittaminen logistisella kasvufunktiolla. Yksinkertaisuudessaan ylempi differentiaaliyhtälö kuvaa saaliseläimen lisääntymistä ja kuolemaa muista kuin petoeläimestä johtuvista syistä vähennettynä petoeläimen saalistamisesta johtuva väheneminen. Alempi differentiaaliyhtälö kuvaa petoeläinkannan kasvua ja vähenemistä. Tässä mallissa funktio h(x) kuvaa petoeläimen saalistamien saaliseläinten lukumäärää, K populaatioiden koko yhteensä, delta petoläinten poistumaa ja gamma lisääntymistä.

Yllä annetussa artikkelissa on erinomaisen syvällistä analyysia aiheesta, ja tässäkin blogissa tullaan ehdottomasti myöhemmin käsittelemään mainittua teoriaa, mutta tässä vaiheessa tyydymme käymään lyhyesti mallin pääkohdat läpi. Smith yksinkertaistaa yllä annetun kaavan Holling II -tyyppiseksi peto-saalis -malliksi, jolloin sen muodoksi saadaan
\begin{equation} \frac{dx}{dt} = x(1-\frac{x}{K})-\frac{mxy}{1+x} \end{equation}
\begin{equation} \frac{dy}{dt} = y(-c + \frac{mx}{1+x})\end{equation}


Tutkitaan seuraavaksi lyhyesti mallin stabiilisuusominaisuuksia. Tasapainopisteet saadaan yhtälöparista
\begin{equation} 0 = x(1-\frac{x}{K}-\frac{my}{1+x})\end{equation}
\begin{equation} 0 = y(-c + \frac{mx}{1+x})\end{equation}

Tässä tapauksessa tasapainopisteet ovat (0,0), (k,0). Smithin \ määritellään myös eläinten välisen yhteiselon tasapainopiste, mutta emme ota sitä tällä kertaa huomioon. Jacobin matriisiksi saamme tasapainopisteellä (k,0)
\begin{equation}J=\begin{bmatrix} -1 & -\frac{mk}{1+k} \\ 0 & \frac{mk}{1+k}-c \end{bmatrix}\end{equation}
mistä huomataan, että piste on epästabiili, sillä toinen ominaisarvoista on aina positiivinen, ja alempikin on negatiivinen vain jos pätee
\begin{equation}\frac{mk}{1+k}<c.\end{equation}
Tasapainopisteen (0,0) Jacobin matriisi on
\begin{equation}J=\begin{bmatrix} 1 & -m \\ 0 & -c \end{bmatrix},\end{equation}
mistä taas huomataan, että pisteet ovat epästabiileja, sillä toinen ominaisarvo on välttämättä positiivinen, sillä oletuksena c > 0.

Ensimmäisessä kuvaajassa on voimassa
\begin{equation}\frac{mk}{1+k}<c,\end{equation} jolloin piste (k,0) ja (0,0) ovat epästabiileja

Alemmassa epäyhtälö on käännetty, jolloin pisteet (k,0) ja (0,0) edelleen epästabiileja, mutta käyrien muoto on täysin erilainen.

Käytetyt koodit, muokattu Lotka-Volterran mallista:

import numpy as np
import matplotlib
matplotlib.use('tkagg')
import matplotlib.pyplot as plt
from scipy.integrate import odeint

y0 = [10,1] # fish, bears, units in hundreds

t = np.linspace(0,50,num=1000)

m = 1.1
c = 2
k = 100

params = [m,c,k]

def sim(variables, t, params):
    # fish population
    x = variables[0]
    # bear population
    y = variables[1]
    
    m = params[0]
    c = params[1]
    k = params[2]
    
    dxdt = x*(1-x/k) - m*x*y/(1+x)
    dydt = -c*y + m*x*y/(1+x)
    
    return([dxdt,dydt])

y = odeint(sim, y0, t, args=(params,))

f,(ax1,ax2) = plt.subplots(2)

line1, = ax1.plot(t,y[:,0], color="b")

line2, = ax2.plot(t,y[:,1], color="r")

ax1.set_ylabel("Fish (hundreds)")
ax2.set_ylabel("Bears (hundreds)")
ax2.set_xlabel("Time")

plt.show()

keskiviikko 17. tammikuuta 2024

Hyvin lyhyt johdatus matemaattiseen mallinnukseen

Mallien johtaminen

Oletetaan, että olemme tehneet joukon havaintoja ilmiöstä A suhteessa B. Huomataan, että havaintojen perusteella voidaan kirjoittaa kaksiulotteinen suora. Havaintojemme perusteella voidaan siis kirjoittaa suoralle yhtälö
\begin{equation}y-y_0 = k(x-x_0) \end{equation}
missä k on suoran kulmakerroin, eli

\begin{equation}k = \frac{\triangle y}{\triangle x}=\frac{y_2-y_1}{x_2-x_1}.\end{equation}
Näin siis yksinkertaisten havaintojen avulla voidaan kirjoittaa ilmiötä yleisesti tai monimutkaisempaa likimain hetkellisesti kuvaava suoran yhtälö, johon havaintoarvoja x asettamalla voidaan arvioida y suuruutta sekä toisin päin.

Otetaan seuraavaksi tarkasteluun logistista mallia käsittelevästä postauksesta tuttu esimerkki eksponentiaalisen kasvun kaavan johtamisesta. Olkoon siis x(t) populaation koko hetkellä t, b syntyneiden määrä ja d kuolleiden määrä hetkellä t. Saadaan

\begin{equation}x(t+\triangle t)=x(t)+b \triangle t x(t) -d \triangle t x(t)  \end{equation}\begin{equation}\frac{x(t+\triangle t)-x(t)}{\triangle t}=(b-d)x(t)  \end{equation}
kun t lähestyy nollaa, saadaan differentiaaliyhtälöksi
\begin{equation}\frac{dx}{dt}=(b-d)x(t)  \end{equation}
sekä differentiaaliyhtälön ratkaisuksi
\begin{equation}x(t)=x_0 e^{rt},  \end{equation}
missä r on kasvuvauhti. Differentiaaliyhtälön ratkaisua kutsutaan eksponentiaaliseksi yhtälöksi. 

Otetaan vielä yksinkertainen diskreetti esimerkki yleisen funktion muotoilemisesta, kun ilmiön tila tunnetaan hetkellä t. Kirjoitetaan
\begin{equation}x(t+1)=x(t)+\alpha, \end{equation}
menemällä ajassa taaksepäin ja sijoittamalla menneitä muuttujia x(t) sisäkkäin, saadaan ensin

\begin{equation}x(t+1)=(x(t-1)+\alpha)+\alpha\end{equation}
\begin{equation}x(t+1)=(x(t-2)+\alpha)+2\alpha,\end{equation}
ja niin edespäin. Menemällä tarpeeksi kauas taaksepäin saadaan yhtälöiden perusteella muotoiltua yleinen muoto:
\begin{equation}x(t+1)=x_0 + (t+1) \alpha \end{equation}

Järjestelmän stabiilisuus

Järjestelmän stabiiliudesta lähde täällä sekä täällä.

Aloitetaan Ljaputinovin stabiiliudesta. Olkoon autonominen dynaaminen systeemi
\begin{equation}x' = f(x(t)), x(0)=x_0, \end{equation}
jolla on olemassa tasapainopiste 
\begin{equation}x_{\epsilon}, f(x_{\epsilon})=0. \end{equation}
Tasapainopistettä sanotaan Ljapunov-stabiiliksi, jos jokaista nollaa suurempaa epsilonia kohden on olemassa myös nolla suurempi delta, jolla pätee kaikilla t >0
\begin{equation} ||x(0)-x_{\epsilon}|| < \delta \Longrightarrow ||x(t)-x_{\epsilon}|| < \epsilon \end{equation}
kaikilla t >0.

Jos tasapainopiste toteuttaa ehdon
\begin{equation} ||x(0)-x_{\epsilon}|| < \delta \Longrightarrow \lim_{t \rightarrow \infty} ||x(t)-x_{\epsilon}|| < \epsilon, \end{equation}
on kyseessä asymptoottisesti stabiili piste.

Asymptoottista tasapainopistettä sanotaan eksponentiaalisesti stabiiliksi, jos mielivaltaisilla muuttujilla alfa,beta ja delta >0 pätee
\begin{equation} ||x(0)-x_{\epsilon}|| < \delta \Longrightarrow ||x(t)-x_{\epsilon}|| < \epsilon \leq \alpha ||x(t)-x_{\epsilon}||e^{- \beta t} \end{equation}
Ehdoista ensimmäinen tarkoittaa, että mikäli yhtälöt lähtevät tasapainopisteen läheisyydestä, ne jäävät kiertämään sen välitöntä ympäristöä kykenemättä erkaantumaan siitä. Asymptoottinen stabiilius tarkoittaa, että yhtälöt eivät vain kierrä tasapainopistettä, vaan myös suppenevat sitä kohti. Eksponentiaalisessa tapauksessa suppenemisnopeus on annettua ylärajaa nopeampaa.

Esimerkkejä Ljapunovin stabiiliuden selvittämisestä löytyy enemmän annetuista lähteistä.

Järjestelmän stabiilisuus voidaan myös selvittää seuraavasti. Olkoon mielivaltainen yhtälöpari
\begin{equation} \frac{dy}{dt}=ax-bz. \end{equation}
\begin{equation} \frac{dz}{dt}=-ax+bz+d. \end{equation}
Ensin selvitetään differentiaaliyhtälöiden tasapainopisteet, eli nollakohdat yllä mainituille yhtälölle:
\begin{equation} 0=ax+b. \end{equation}
\begin{equation} 0 =-ax+bz+d. \end{equation}
Seuraavaksi kirjoitetaan samaisen differentiaaliyhtälökaksikon Jacobin matriisi

\begin{equation}J=\begin{bmatrix} a & 0 \\ a & b \end{bmatrix}\end{equation}
josta laskettujen ominaisarvojen

\begin{equation} det(J-\lambda I) = 0 \end{equation}
avulla saadaan selvitettyä funktion keskeiset stabiilisuusominaisuudet. Eli siis jos jatkuvan dynaamisen järjestelmän ominaisarvot ovat reaalisia ja

\begin{equation}\lambda_{1,2} < 0,\end{equation}
on järjestelmä stabiili, kun taas jos on voimassa
\begin{equation}\lambda_{1,2} > 0\end{equation}
\begin{equation}\lambda_{1} < 0, \lambda_{2}>0,\end{equation}
on järjestelmä epästabiili. Lisäksi jos reaaliosa on nolla, on järjestelmä silloinkin epästabiili. Jos ominaisarvot ovat imaginäärisiä, ja reaalinen osa on negatiivinen, on järjestelmä stabiili. Muutoin kyseessä on epästabiili järjestelmä. Diskreetissä tapauksessa järjestelmä on asymptoottisesti stabiili jos ja vain jos ominaisarvon reaaliselle osalle pätee

\begin{equation}|\lambda_{1,2}|<1.\end{equation}

Jatkuvassa tapauksessa järjestelmä on asymptoottisesti stabiili jos ja vain jos ominaisarvojen reaaliosat ovat negatiivisia.

Nämä tiedot riittävät tässä vaiheessa mallinnusta koskevien blogipostausten ymmärtämiseen. Tilastollista mallintamista, optimointia ja todennäköisyyslaskentaa koskevia osuuksia käsitellään erikseen myöhemmin.

sunnuntai 14. tammikuuta 2024

SIR- ja SIRV-mallit sekä rokotuksen vaikutuksen mallinnus

SI-mallin johtaminen

Aiheen johdannoksi esitellään Charnovin monisteiden esittelemä tapa tartuntatautimallin määrittelemiseksi. Aloitetaan seuraavasti tunnetuksi oletetusta kaavasta
\begin{equation}I(t+\triangle t) = I(t)+ \beta \triangle t S(t) I(t)\end{equation}
missä t kuvaa aikaa, I(t) tartuntojen määrää ja S(t) alttiiden määrää hetkellä t sekä beta kontaktien määrää aikayksikössä, eli leviämisintensiteettiä. Viemällä I(t) funktion vasemmalle puolelle, jakamalla molemmat puolet kolmio-t:llä, kun kolmio-t lähestyy nollaa, saadaan tartuntojen määrän muutosta kuvaava differentiaaliyhtälö
\begin{equation} \frac{dI}{dt}=\beta SI, \end{equation}
joka on tuttu myöhemmistä aiheen variaatioista. Myöhempiä malleja ei johdeta yhtä perusteellisesti, vaan käytännössä aloitetaan valmiista ja annetuiksi oletetuista differentiaaliyhtälöistä.

SIR-malli

SIR-malli kuvaa epidemian dynamiikkaa olettamalla olevaksi kolme tilaa: alttiiden S, todetun tautitapauksen I sekä parantuneen R. Tässä mallissa ei oteta huomioon syntymiä eikä kuolemia. Kerran saatu immuniteetti ei häviä, eikä samaan tautiin voi sairastua toistamiseen. Näistä oletuksista on erilaisia vapaampia muotoiluja, kuten SIRI-malli, jossa henkilölle ei muodostu kestävää immuniteettiä, SIRD-malli, jossa henkilön kuolema otetaan huomioon, sekä joukko erilaisia rokotusmalleja, joista tässä käsiteltävässä immuniteetti on pysyvä. Vaihtoehtoisia muotoiluja oletuksineen käsitellään erikseen myöhemmin.

Yksinkertainen SIR-malli muodostuu differentiaaliyhtälöistä

\begin{equation} \frac{dS}{dt}=\frac{-\beta SI}{N} \end{equation}

\begin{equation} \frac{dI}{dt}=\frac{\beta SI}{N} -\gamma I\end{equation}

\begin{equation} \frac{dR}{dt}=\gamma I. \end{equation}
SIR-mallin parametreista beta kuvaa tartuntaintensiteettiä ja gamma poistumis- tai parantumisintensiteettiä riippuen miten mallia tulkitaan. Kahden parametrin avulla voidaan edelleen määritellä perustarttuvuusluku

\begin{equation} R_0 = \frac{\beta}{\gamma}, \end{equation}
joka karkeasti tarkoittaa odotettujen tartuntojen määrän. Seuraavaksi huomataan, että sijoittamalla perustarttuvuusluku tartuntojen muutosta kuvaavaan differentiaaliyhtälöön, saadaan

\begin{equation}\frac{dI}{dt}=(R_0 \frac{S}{N}-1)\gamma I, \end{equation}
mistä edelleen huomataan, että kun pätee

\begin{equation}R_0 S > N, \end{equation}
on infektioiden muutos samalla kasvava, eli tartuntojen määrä lisääntyy ja epidemia on kiihtymisvaiheessa. Käänteisessä tapauksessa tartuntojen määrä vähenee ja epidemia on asettumassa. Sama voitaisiin todeta toisin, eli olettamalla likimain jokaisen populaatioon kuuluvan yksilön olevan altis eli S = N, seuraisi siitä, että epidemian käynnistymiseen riittää ehto
\begin{equation}R_0 > 1, \end{equation}
sekä talttumiseen käänteinen tulos.

SIR-mallin kuvaajassa näkyvät eri ryhmien muutosten suunnat sekä kasvunopeudet. Seuraavan SIR-mallin kuvaajaan löytyy koodi täältä.


Kuvaajasta huomataan, että alttiiden määrä vähenee taudin leviämisen sekä siitä parantuneiden kasvamisen mukana. Kaikki muutosnopeudet tasaantuvat ajan mittaan. Samaista koodia sovelletaan mahdollisuuksien mukaan myös saman mallin eri sovelluksissa, kuten seuraavassa SIRV-mallissa.

SIRV-malli

Määritellään johdannoksi rokotusmalliin seuraava tulos. Olkoon S_0 alttiiden prosentuaalinen osuus hetkellä 0, olkoon lisäksi p rokotettujen osuus väestöstä, ja olkoon lisäksi voimassa S_0=(1-p). Jos oletetaan rokotuksen tuoma immuniteetti pysyväksi, voidaan kirjoittaa edellisen kohdan perustarttuvuuslukua hyödyntäen seuraava ehto epidemian käynnistymiselle:
\begin{equation} R_0 (1-p) > 1. \end{equation}
Annetun ehdon avulla voidaan määritellä rokotettavien minimäärä p* epidemian estämiseksi, jolloin voidaan kirjoittaa
\begin{equation} p^* = 1-\frac{1}{R_0}.\end{equation}
SIRV-malli on perusteiltaan kuin aiemmin kuvailtu yksinkertainen SIR-malli, mutta siinä on otettu huomioon rokotuksen vaikutus. Rokotus antaa tässä mallissa pysyvän immuniteetin, eikä tautiin kuolleita oteta tässäkään huomioon.

\begin{equation} \frac{dS}{dt}=\frac{-\beta SI}{N} -vS\end{equation}

\begin{equation} \frac{dI}{dt}=\frac{\beta SI}{N} -\gamma I\end{equation}

\begin{equation} \frac{dR}{dt}=\gamma I \end{equation}

\begin{equation} \frac{dV}{dt}=v S \end{equation}
Tässä mallissa parametrit tarkoittavat samaa kuin edellisessä mallissa, mutta uutena SIRV-malliin on lisätty rokotusintensiteettiä kuvaava parametri v. Kun muut parametrit ovat samat kuin ensimmäisessä ja rokotusintensiteetti v = 0.005, saadaan kuvaaja


Kuvasta siis näemme, että rokotuksilla saadaan taudin leviäminen hallintaan ja immuniteetti puolelle samasta määrästä ihmisiä, jotka saivat saman sairastamalla taudin. Infektoiden määrän kasvuvauhtia saatiin näin loivennettua ja huippu madallettua.


SIRV-mallissa käytetty Python-koodi, muokattu aiemmin linkatusta esimerkistä:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

N = 1000
I0, R0, V0 = 1, 0, 0
S0 = N - I0 - R0 -V0
beta, gamma, v = 0.2, 1./10, 0.005
t = np.linspace(0, 160, 160)

def deriv(y, t, N, beta, gamma, v):
    S, I, R, V = y
    dSdt = -beta * S * I / N - v * S
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    dVdt = v * S
    return dSdt, dIdt, dRdt, dVdt

y0 = S0, I0, R0, V0

ret = odeint(deriv, y0, t, args=(N, beta, gamma,v))
S, I, R, V = ret.T

fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
ax.plot(t,V/1000, 'y', alpha=0.5, lw=2, label = 'vaccinated' )
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()

Logistiset kasvumallit

Diclaimer: Bloggerin matematiikkatyökalu on järjettömän buginen ja epäluotettava. Siksi olenkin tekstissäni pyrkinyt mahdollisuuksien mukaan välttämään ylimääräisten symbolien käyttöä.

Lähdetään oletuksesta, että N(t) on populaation koko hetkellä t, b on syntyneiden määrä ja d kuolleiden määrä. Saadaan

\begin{equation}N(t+\triangle t)=N(t)+b \triangle t N(t) -d \triangle t N(t)  \end{equation}

\begin{equation}\frac{N(t+\triangle t)-N(t)}{\triangle t}=(b-d)N(t)  \end{equation}

kun t lähestyy nollaa, saadaan differentiaaliyhtälöksi

\begin{equation}\frac{dN}{dt}=(b-d)N(t)  \end{equation}

sekä differentiaaliyhtälön ratkaisuksi

\begin{equation}N(t)=N_0 e^{rt},  \end{equation}
missä r on kasvuvauhti. Kyseessä on toisin sanoen eksponentiaalisesti kasvava yhtälö. Malthusin lähde, Chasnovin Mathemathical Biology -moniste löytyy ilmaiseksi täältä. Seuraavassa yksinkertainen diskreetti esimerkki eksponentiaalisen kasvun mallin kuvaajasta:

Perusominaisuuksiltaan eksponentiaalinen malli sopii erinomaisesti taloudellisen kasvun ja kassaantumisen kuvaamiseen, mutta resurssien luonnossa asettamien rajoitusten vuoksi erityisen huonosti populaatioiden kasvun mallintamiseen. Tästä syystä Malthusin kasvumallin pohjalta kehitettiin Verholstin malli eli logistinen kasvumalli, jossa otetaan huomioon kasvun yläraja K:

\begin{equation}\frac{dN}{dt}=F(N)=rN(1-\frac{N}{K}). \end{equation}
Verhulstin mallissa N kuvaa populaation kokoa, r kasvuvakiota ja K populaation koon ylärajaa. Logistisen kasvumallin kasvunopeus hidastuu ajan ja populaation koon kasvun mukana. Logistisen kasvumallin differentiaaliyhtälölle voidaan muotoilla seuraavanlainen kuvaaja populaation muutosnopeuden ja koon mukaan:


Logistinen kasvumalli ratkaisemalla differentiaaliyhtälöiden ratkaisumenetelmiä hyödyntäen päädytään kaavaan

\begin{equation}N(t)=\frac{P_0*K*e^{r*t}}{(K-P_0)+P_0*e^{rt}} \end{equation}

jonka pohjalta saadaan muotoiltua populaation koon kuvaaja:

 

Logistiselle kasvumallille saadaan tasapainopisteiksi N=0 ja N=K. Tasapainopisteiden stabiilius saadaan tarkistettua derivoimalla logistinen malli, jolloin saadaan
\begin{equation}F'(N)=r(1-\frac{2N}{K}). \end{equation}

Sovittamalla tasapainopisteet derivoituun yhtälöön, huomataan, että kun N=0, yhtälö saa muodon F'(0)=r ja kun F'(K)=-r. Kyseessä on toisin sanoen satulapiste, jossa nollapiste on epästabiili sekä K stabiili. Populaatio voi siis lähestyä nollapistettä mielivaltaisen lähelle kyeten kuitenkin kasvamaan siitä pois, mutta lähestyminen kapasiteetin ylärajaa (voi) johtaa ajautumiseen kohti ääripistettä.

Logistisella mallilla on olemassa myös Allee-efektillä täydennetty muoto

\begin{equation}\frac{dN}{dt}=-rN(1-\frac{N}{A})(1-\frac{N}{K}), \end{equation}
missä populaation kasvuvauhti on negatiivinen kun 0<A<N sekä positiivinen kun N<A<K, oletuksella että on aina voimassa A<K. Ensimmäisessä tapaksessa Allee-efektiä sanotaan "vahvaksi", jälkimmäisessä "heikoksi".  Tällöin populaation koolla on siis vaikutusta kasvunopeuteen, mutta vaikutuksen suuruus ja suunta riippuvat muuttujan A arvoista. Allee-efektillä varustetun logistisen kasvumallin kuvaajat ovat vahvalla Allee-efektillä kun K=102 ja A=40:

sekä heikolla kun A = 10 ja N = 11-100:


Allee-efektillä täydennetyn mallin stabiilisuusominaisuudet seuraavat pitkälti tavallisen logistisen mallin vastaavia, joten niitä ei tässä erikseen käydä läpi. Aiheesta on syvällisemmin kirjoittanut Lv, Chen ja Chen artikkelissaan.

Wallace ja Tsoularis muotoilivat esseessään seuraavan yleisen muodon logistisen kasvun funktiolle:

\begin{equation} \frac{dN}{dt}=rN^{\alpha}(1-(\frac{N}{K})^{\beta})^{\gamma} .\end{equation}

Huomataan, että kun alfa = 1 ja gamma = 0, saadaan eksponentiaalinen kasvumalli. Kun taas alfa = 1, beta=1 ja gamma=1, saadaan tavallinen logistinen kasvumalli. Muuttujien eri arvoilla voidaan edelleen muotoilla suoraan ainakin Gompertzin malli, Blumbergin malli, Hyper-Gompertz... käytännössä kaikki kuvaavat samanlaista käyrää, mutta pienin erilaisin painotuksin. Mainittuja malleja ei tässä käsitellä enempää.

Käytetyt R-koodit:


eksponentiaalifunktio <- function(r,t){
    funktio <- r*exp(t)
    return(funktio)
}

plot(x=1:7,y=eksponentiaalifunktio(0.1,1:7),type = "l", col="red",xlab="aika",ylab = "populaation koko")
 
logistinen_funktio <- function(r,N,K){
funktio <- r*N*(1-N/K)
return(funktio)
}


plot(x=0:100,y=logistinen_funktio(0.5,0:100,100),type = "l", col="red",xlab="N",ylab = "F(N)")
 
logistinen_funktio1 <- function(p_0,K,t,r){
funktio <- p_0*K*exp(r*t)/((k-p_0)+p_0*exp(r*t))
return(funktio)
} 
 
plot(x=0:30,y=logistinen_funktio1(1,100,0:30,0.5),type = "l", col="red",xlab="aika",ylab = "populaation koko")

allee_funktio <- function(r,N,K,A){
    funktio <- -r*N*(1-N/A)*(1-N/K)
    return(funktio)
}

plot(x=0:100,y=allee_funktio(0.5,0:100,102,40),type = "l", col="red",xlab="N",ylab = "F(N)")
lines(x, y2, type = "l", col = "black")   
plot(x=0:89,y=allee_funktio(0.5,11:100,102,10),type = "l", col="red",xlab="N",ylab = "F(N)")
lines(x, y2, type = "l", col = "black")