Stochastyczny i deterministyczny model logistycznego wzrostu

Łańcuchy Markowa z czasem ciągłym

Rozważamy stochastyczny model logistycznego wzrostu. Podobnie jak dla modelu SIR wykorzystamy wygenerowane trajektorie, aby określić własności procesu. Rezultaty zestawimy z (deterministycznym) rozwiązaniem odpowiedniego równania różniczkowego zwyczajnego.

W wersji deterministycznej modelem logistycznego wzrostu nazywamy rozwiązanie

\(n = n(t)\) równania różniczkowego zwyczajnego \[ \frac{dn}{dt}=rn(1-\frac{n}{K}), \] gdzie \(n(0)=n_0>0\) oraz \(r,K\) są parametrami.

Interpretacja jest następująca: \(n(t)\) opisuje wielkość populacji w chwili \(t\), \(r\) oznacza intensywność wzrostu populacji („przyrost naturalny”), a \(K\) pojemność środowiska.

W modelu stochastycznym modelujemy wprost „każdego osobnika” w populacji; wielkość populacji \(n(t)\) określa liczbę osobników, z których każdy może się rozmnożyć \((n \rightarrow n + 1)\) lub umrzeć \((n \rightarrow n - 1)\). Załóżmy, że w chwili \(t\) mamy \(n(t) = i\geq 0\) osobników oraz, że z intensywnością \(\lambda_i\) w populacji pojawia się nowy osobnik \((i \rightarrow i + 1)\), a z intensywnością \(\mu_i\) jeden z osobników umiera \((i \rightarrow i - 1)\).

Oczywiście \(\lambda_0=\mu_0=0\).

Stochastyczny model logistyczny powinien spełniać następujące warunki:

\[\lambda_n-\mu_n = rn-\frac{r}{K}n^2,\] \[\lambda_K = \mu_K.\] Rozważmy dwa różne przypadki spełniające powyższe założenia: \[ (i)= \left\{ \begin{array}{ll} i-\frac{i^2}{100} & i=0,1,\ldots ,100,\\ 0 & i>100,\\ \end{array} \right. \] \[\mu_i=\frac{i^2}{100} , i=1,2,\ldots\] \[ (ii)\lambda_i=i, \ \mu_i=\frac{i^2}{50}, \ i=0,1,2,\ldots \]

Przyjmijmy, że \(n_0=5.\)

Dla modelu deterministycznego, przedstaw rozwiązanie numeryczne \(n(t)\) (pakiet deSolve), dla różnych wartości parametrów \(r\), \(K\) oraz różnych warunków początkowych \(n_0\).
Następnie dla modelu stochastycznego: wyznacz wartość oczekiwaną \(E[n(t)]\) oraz wariancję \(var[n(t)]\) dla ustalonej chwili np. \(t = 10\);
Wyznacz funkcję wartości oczekiwanej \(E[n(t)]\) oraz wariancji \(var[n(t)]\), \(t \geq 0\).
Porównaj \(E[n(t)]\) z rozwiązaniem deterministycznym.
Opisz różnice pomiędzy modelem deterministycznym oraz modelem stochastycznym (oraz między wersjami (i) oraz (ii)).

Teoria

Zanim przejdziemy do właściwego rozwiązania zadania, zbadajmy nasze równanie, w celu wyznaczenia stanów stacjonarnych, które pozwolą nam w pełni zrozumieć analizę modelu. Zacznijmy zatem od znalezienia stanów stacjonarnych.

\[rn(1-\frac nK)=0 \Leftrightarrow n=0 \lor n=K\] Zatem równanie posiada dwa stany stacjonarne: \(n_1=0\) i \(n_2=K.\) Badając funkcje znaku, możemy wyczytać, że rozwiązanie \(n(t)\) rośnie w przedziale \([0,K]\) oraz maleje dla \(t\geq K\).

portret fazowy

Z powyższego portretu fazowego odczytujemy, że punkt \(n_1=0\) jest niestabilny, natomiast punkt \(n_2=K\) jest asymptotycznie stabilny.

Następnie obliczamy drugą pochodną w celu zbadania wypukłości funkcji. \[\frac{d^2n}{dt^2}=r\frac{dn}{dt}\Big(1-\frac{2n}{K}\Big).\] Po przyrównaniu drugiej pochodnej do zera, otrzymujemy, że funkcja \(n(t)\) ma punkt przegięcia o wartości \(\frac K2\). Zatem możemy sformułować następujące wnioski:

rozwiązanie naszego równania zależy od początkowej liczebności \(n_0\);
dla \(n_0>K\) rozwiązanie maleje i dąży do punktu równowagi \(K\);
gdy \(n_0=K\) populacja pozostaje na stałym poziomie \(K\);
dla \(n_0<K\) funkcja rośnie i dąży do punktu równowagi \(K\), przy czym w przedziale \((\frac K2, K)\) mamy powolny wzrost, a dla \(n_0<\frac K2\) mamy najpierw szybki wykładniczy wzrost, a następnie tempo tego wzrostu spada.

Model deterministyczny

W celu przedstawienia rozwiązania numerycznego \(n(t)\) dla modelu deterministycznego zainstalowano pakiet “deSolve” i zbadano rozwiązania w zależności od różnych wartości parametrów r,K oraz różnych warunków początkowych \(n_0\).

\(n_0 = 5\)-stałe, \(r = 2\)-stałe, K-zmienne

#a) 
#Instalujemy pakiet "deSolve" w celu przedstawienia rozwiązania numerycznego N(t),dla różnych wartości
#parametrów r, K oraz różnych warunków początkowych N (n0).

#install.packages("deSolve")
library("deSolve")

## Warning: pakiet 'deSolve' został zbudowany w wersji R 4.1.3

T=50

Markov_chain<- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dN<- r*N*(1-N/K)
    return(list(c(dN)))
  })
}


pars <- c(r=2,K=10)

yini <- c(N=5)
times <- seq(0, 50, by = 0.01)
out <- ode(yini, times, Markov_chain, pars)
#head(out)
plot(out[,1],out[,2],type="l",xlim=c(0,T),xlab="Czas (t)",ylab="Wielkosc populacji",
     main="Wielkosc populacji w chwili t - Model deterministyczny-n0=5,r=2,K=10")

pars <- c(r=2,K=5)

yini <- c(N=5)
times <- seq(0, 50, by = 0.01)
out <- ode(yini, times, Markov_chain, pars)
#head(out)
plot(out[,1],out[,2],type="l",xlim=c(0,T),xlab="Czas (t)",ylab="Wielkosc populacji",
     main="Wielkosc populacji w chwili t - Model deterministyczny-n0=5,r=2,K=5")

pars <- c(r=2,K=3)

yini <- c(N=5)
times <- seq(0, 50, by = 0.01)
out <- ode(yini, times, Markov_chain, pars)
#head(out)
plot(out[,1],out[,2],type="l",xlim=c(0,T),xlab="Czas (t)",ylab="Wielkosc populacji",
     main="Wielkosc populacji w chwili t - Model deterministyczny-n0=5,r=2,K=3")

\(n_0 = 5\)-stałe, \(K = 20\)-stałe, r-zmienne

pars <- c(r=20,K=20)

yini <- c(N=5)
times <- seq(0, 50, by = 0.01)
out <- ode(yini, times, Markov_chain, pars)
#head(out)
plot(out[,1],out[,2],type="l",xlim=c(0,T),xlab="Czas (t)",ylab="Wielkosc populacji",
     main="Wielkosc populacji w chwili t - Model deterministyczny-n0=5,r=20,K=20")

pars <- c(r=30,K=20)

yini <- c(N=5)
times <- seq(0, 50, by = 0.01)
out <- ode(yini, times, Markov_chain, pars)
#head(out)
plot(out[,1],out[,2],type="l",xlim=c(0,T),xlab="Czas (t)",ylab="Wielkosc populacji",
     main="Wielkosc populacji w chwili t - Model deterministyczny-n0=5,r=30,K=20")

\(n_0\) -zmienne, \(K = 50\)-stałe, \(r=2\)-stałe

pars <- c(r=2,K=50)

yini <- c(N=15)
times <- seq(0, 50, by = 0.01)
out <- ode(yini, times, Markov_chain, pars)
#head(out)
plot(out[,1],out[,2],type="l",xlim=c(0,T),xlab="Czas (t)",ylab="Wielkosc populacji",
     main="Wielkosc populacji w chwili t - Model deterministyczny-n0=15,r=2,K=50")

pars <- c(r=2,K=50)

yini <- c(N=35)
times <- seq(0, 50, by = 0.01)
out <- ode(yini, times, Markov_chain, pars)
#head(out)
plot(out[,1],out[,2],type="l",xlim=c(0,T),xlab="Czas (t)",ylab="Wielkosc populacji",
     main="Wielkosc populacji w chwili t - Model deterministyczny-n0=35,r=2,K=50")

pars <- c(r=2,K=50)

yini <- c(N=75)
times <- seq(0, 50, by = 0.01)
out <- ode(yini, times, Markov_chain, pars)
#head(out)
plot(out[,1],out[,2],type="l",xlim=c(0,T),xlab="Czas (t)",ylab="Wielkosc populacji",
     main="Wielkosc populacji w chwili t - Model deterministyczny-n0=75,r=2,K=50")

Uwzględniając różne możliwości otrzymujemy następujące wnioski:

Przy zadanych parametrach \(n_0\) i \(r\), zauważamy, że zmieniając parametr \(K\) (pojemność środowiska) uzyskujemy wartość wielkości populacji na różnych poziomach. Ponieważ \(K\) jest punktem stałym rozwiązania równania różniczkowego, rozwiązanie zawsze stabilizuje się na poziomie \(K\), to znaczy, że populacja osiąga stan równowagi w zadanym przez nas punkcie. Zatem zmieniając pojemność środowiska ustalamy końcową wielkość populacji.
Przy zadanych parametrach \(n_0\) i \(K\), wraz ze wzrostem parametru \(r\) (intensywność wzrostu populacji) populacja szybciej dąży do zadanego punktu równowagi. Interpretując \(r\) jako przyrost naturalny, zwiększanie tego parametru oznacza większą intensywność rozmnażania się nowych osobników. Na wykresie objawia się to krzywą, która co raz bardziej przybliża się do osi \(Y\).
Przy zadanych parametrach \(r\) i \(K\) zauważamy, że dobierając \(n_0\) (wielkość populacji) w przedziale \([0,\frac{K}{2}]\) funkcja początkowo szybko rośnie, stabilizując się w punkcie równowagi. Dla \(n_0\) w przedziale \((\frac{K}{2},K]\) funkcja również rośnie do momentu osiągnięcia punktu równowagi. Natomiast dla \(n_0>K\) funkcja maleje do momentu stabilizacji. Dla \(n_0=K\) wielkość populacji nie ulega zmianie.

Model Stochastyczny

Wersja (i)

Wyznaczamy funkcję wartości oczekiwanej i wariancji dla zadanego procesu.

Poniżej zaprezentowano wykresy dla \(n_0=5\).

#WERSJA (i)
T=50
#Funkcja generujaca trajektorie zadanego procesu
trajektorie1 <- function(N){
  t=c(0)
  i=1
while(t[i]<T & N[i]>0){
  if(N[i]<100){
  lambda=N[i]-N[i]^2/100
  }else{
    lambda=0}
  mi=N[i]^2/100
  t[i+1]=t[i]+rexp(1, lambda+mi)
  
  if(runif(1)< lambda/(lambda+mi)){
    N[i+1]=N[i]+1
    
  }else{
    N[i+1]=N[i]-1
    
  }
  
  i=i+1
}
 return(cbind(t,N))
}

#c)
#Wyznaczamy wartosc oczekiwana i wariancje
#Narysujemy średnią wartość m trajektorii procesu -> w tym celu każdą wygenereowaną trajektorię
#Zapiszemy w dyskretnym podziale czasu t_i [0,T] co dt, te trajektorie będziemy dodawać do siebie.

  sumy<-function(m,N,T,dt){
  dysk=seq(0,T,by=dt)
  trajsuma=numeric(length(dysk))
  trajsuma2=numeric(length(dysk))
  for(i in 1:m){
    t=c(0)
    i=1
    while(t[i]<T & N[i]>0){
      if(N[i]<100){
        lambda=N[i]-N[i]^2/100
      }else{
        lambda=0}
      mi=N[i]^2/100
      t[i+1]=t[i]+rexp(1, lambda+mi)
      
      if(runif(1)< lambda/(lambda+mi)){
        N[i+1]=N[i]+1
        
      }else{
        N[i+1]=N[i]-1
        
      }
      
      i=i+1
    }
    trajdysk<-sapply(dysk, function(x) N[sum(t<=x)])
    trajsuma=trajsuma+trajdysk
    trajsuma2=trajsuma2+trajdysk^2
  }
  return(cbind(dysk,srednia=trajsuma/m,srednia2=trajsuma2/m))
}
#Otrzymujemy dyskretny podział czasu [0,T], średnią N(t) i średnią kwadratów (N^2(t))
T=50
N=5
dt=0.1
m=1000
wynik=sumy(m,N,T,dt)
plot(wynik[,1],wynik[,2],type='l',lwd=2,xlab='Czas (t)',ylab='E(N(T))',xlim=c(0,T),
     main="Wartosc oczekiwana")

#Wariancja = (EN)^2 - (E(N^2))
plot(wynik[,1],wynik[,3]-wynik[,2]^2,type='l',lwd=2,xlab='t',ylab='Var(N(T))',main="Wariancja")

Dla porównania zaprezentowano również wykres pojedynczej trajektorii.

#Rysujemy wykres trajektorii startujac od N=5
traj1=trajektorie1(5)
plot(traj1[,1],traj1[,2],type='s',xlab="Czas (t)",ylab="Wielkosc populacji (N(t))",
     main="Wielkosc populacji w chwili t")

Na wykresie pojedynczej trajektorii widać, że wielkość populacji po pewnym czasie stabilizuje się, ale nie jest stała. Wahania liczebności osobników oscylują wokół rozwiązania równania różniczkowego, czyli do punktu równowagi.

W celu wyznaczenia funkcji wartości oczekiwanej i wariancji dodajemy i bierzemy średnią z wielu trajektorii, dyskretyzując czas w przedziale [0,T].

Na pierwszym wykresie widzimy, że wartość oczekiwana dąży do punktu równowagi. W porównaniu do pojedynczej trajektorii, wahania liczebności są znacznie mniejsze. Na wykresie wariancji widać, że po pewnym czasie funkcja stabilizuje się.

W celu wyliczenia wartości oczekiwanej i wariancji w konkretnej chwili t, np.t=10 podstawiamy do funkcji zadaną wartość t otrzymując następujące wyniki:

#b)
wynik[wynik[,1]==10,2]                             #wartosc oczekiwana dla czasu t=10

## srednia 
##  49.301

wynik[wynik[,1]==10,3]-wynik[wynik[,1]==10,2]^2    #wariancja dla czasu t=10

## srednia2 
##  22.9664

Zatem spodziewamy się, że populacja, która w chwili t=0 ma 5 osobników po czasie t=10 będzie miała ich około 49.

Wersja (ii)

Wyznaczamy funkcję wartości oczekiwanej i wariancji dla zadanego procesu.

Poniżej zaprezentowano wykresy dla \(n_0=5\).

T=50
#WERSJA (ii)
trajektorie2<- function(N){
  t=c(0)
  i=1
  while(t[i]<T & N[i]>0){
    lambda=N[i]
    mi=N[i]^2/50
    t[i+1]=t[i]+rexp(1, lambda+mi)
    
    if(runif(1)< lambda/(lambda+mi)){
      N[i+1]=N[i]+1
      
    }else{
      N[i+1]=N[i]-1
      
    }
    
    i=i+1
  }
  return(cbind(t,N))
}
traj2=trajektorie2(5)

#c)
sumy<-function(m,N,T,dt){
  dysk=seq(0,T,by=dt)
  trajsuma=numeric(length(dysk))
  trajsuma2=numeric(length(dysk))
  for(i in 1:m){
    t=c(0)
    i=1
    while(t[i]<T & N[i]>0){
      lambda=N[i]
      mi=N[i]^2/50
      t[i+1]=t[i]+rexp(1, lambda+mi)
      
      if(runif(1)< lambda/(lambda+mi)){
        N[i+1]=N[i]+1
        
      }else{
        N[i+1]=N[i]-1
        
      }
      
      i=i+1
    }
    trajdysk<-sapply(dysk, function(x) N[sum(t<=x)])
    trajsuma=trajsuma+trajdysk
    trajsuma2=trajsuma2+trajdysk^2
  }
  return(cbind(dysk,srednia=trajsuma/m,srednia2=trajsuma2/m))
}
#Wartosc oczekiwana
T=50
N=5
dt=0.1
m=1000
wynik2=sumy(m,N,T,dt)
plot(wynik2[,1],wynik2[,2],type='l',lwd=2,xlab='Czas (t)',ylab='E(N(T))',xlim=c(0,T),
     main="Wartosc oczekiwana")

#Wariancja = (EN)^2 - (E(N^2))
plot(wynik2[,1],wynik2[,3]-wynik2[,2]^2,type='l',lwd=2,xlab='t',ylab='Var(N(T))',main = 'Wariancja')

Dla porównania zaprezentowano również wykres pojedynczej trajektorii.

#Rysujemy wykres trajektorii startujac od N=5
traj2=trajektorie2(5)
plot(traj2[,1],traj2[,2],type='s',xlab="Czas (t)",ylab="Wielkosc populacji (N(t))",
     main="Wielkosc populacji w chwili t")

Podobnie jak w przypadku modelu (i) pojedyncza trajektoria prezentuje wielkość populacji, która oscyluje wokół punktu równowagi, natomiast funkcja wartości oczekiwanej charakteryzuje się mniejszymi wahaniami i widocznie zbiega rozwiązania równania różniczkowego. W przypadku funkcji wariancji, po pewnym czasie ustala się ona na stałym poziomie, wyższym niż w modelu (i).

W celu wyliczenia wartości oczekiwanej i wariancji w konkretnej chwili t, np.t=10 podstawiamy do funkcji zadaną wartość t otrzymując następujące wyniki:

#b)
wynik2[wynik2[,1]==10,2]                             #wartosc oczekiwana dla czasu t=10

## srednia 
##  49.503

wynik2[wynik2[,1]==10,3]-wynik2[wynik2[,1]==10,2]^2    #wariancja dla czasu t=10

## srednia2 
## 51.66799

Zatem spodziewamy się, że populacja, która w chwili t=0 ma 5 osobników po czasie t=10 będzie miała ich około 49.

wersja(i)

W podpunkcie d) porównano wartość oczekiwaną z rozwiązaniem deterministycznym i zaprezentowano wynik na wykresie. Kolorem czerwonym zaznaczono wartość deterministyczną.

Parametry rozwiązania deterministycznego wyliczono w następujący sposób:

K wyliczono przyrównując \(\mu_i\) i \(\lambda_i\), wyznaczając wielkość \(i\),
r wyliczono podstawiając do równości \(\lambda_n-\mu_n = rn-\frac{r}{K}n^2\), biorąc n=5 i K=50.

pars <- c(r=1,K=50)

yini <- c(N=5)
times <- seq(0, 50, by = 0.01)
out <- ode(yini, times, Markov_chain, pars)

plot(wynik[,1],wynik[,2],type='l',lwd=2,xlab='Czas (t)',ylab='E(N(T))',xlim=c(0,T),
     main="Wartosc oczekiwana")
lines(out[,1],out[,2],type="l",col="red")

Widzimy, że wartość oczekiwana modelu stochastycznego dąży do rozwiąznia modelu deterministycznego.

wersja(ii)

W podpunkcie d) porównano wartość oczekiwaną z rozwiązaniem deterministycznym i zaprezentowano wynik na wykresie. Kolorem niebieskim zaznaczono wartość deterministyczą.

plot(wynik2[,1],wynik2[,2],type='l',lwd=2,xlab='Czas (t)',ylab='E(N(T))',xlim=c(0,T),
     main="Wartosc oczekiwana")
lines(out[,1],out[,2],type="l",col="blue")

Podobnie jak w modelu (i) widzimy, że wartość oczekiwana modelu stochastycznego dąży do rozwiąznia modelu deterministycznego.

Podsumowanie

Różnice pomiędzy modelem deterministycznym a stochastycznym.

Podstawową różnicą pomiędzy tymi modelami jest to, że model stochastyczny za każdym uruchomieniem programu wyznacza różne wartości badanych wielkości. Natomiast model deterministyczny wyznacza jedno, teoretyczne rozwiązanie. Ze względu na konstrukcję parametrów \(\lambda\) oraz \(\mu\) wykresy pojedynczych trajektorii, w przypadku modelu stochastycznego, charakteryzują się pewnymi wahaniami, natomiast rozwiązanie równania różniczkowego prezentuje “gładki” wykres i dodatkowo wielkości \(n(t)\) przyjmują wartości w liczbach rzeczywistych, a nie tylko w liczbach naturalnych.

Różnice pomiędzy modelem (i) a (ii).

Zważając na różnicę w definiowaniu parametrów zauważamy, że w modelu (ii) uzyskujemy wyższą wariancję niż w modelu (i), w modelu (i) po pewnym czasie oscyluje ona w przedziale 20-30, natomiast w modelu (ii) w przedziale 40-50. Rożnice te wynikają z tego, że w okolicach punktu równowagi ułamek \(\frac{\lambda}{\lambda+\mu}\) (który to odpowiada za zmianę stanu) w modelu (ii) wolniej zmienia wartość. To oznacza, że częściej będą występować sytuację, że zanim umrze jeden osobnik, musi narodzić się więcej osobników (lub na odwrót). W modelu (i) po narodzinach jednego osobnika ułamek \(\frac{\lambda}{\lambda+\mu}\) zmniejsza się, co sprawia, że w następnej chwili prawdopodobieństwo śmierci rośnie (tym samym prawdopodobnieństwo narodzin spada) i wahania będą mniejsze niż w modelu (ii).