2018

Dynamical System

  • Recall the system of ordinary differential equations (ODEs) \[\frac{d\mathbf{x}(t)}{dt} = F(\mathbf{x}(t)).\]

  • Interaction amongs the variables \(\mathbf{x}=(x_1,x_2,x_3)\) and \(N=x_1+x_2 + x_3\)

  • Consider the following system \[\frac{dx_1}{dt} = - \beta \frac{x_1 x_2}{N} \\ \frac{dx_2}{dt} = \beta \frac{x_1 x_2}{N} - \gamma x_2 \\ \frac{dx_3}{dt} = \gamma x_2 \mbox{ and } x_2\geq0\]

Solving the system

library("deSolve")
ODEfunc=function(t, x, vparameters){
   X_1 = x[1]; X_2 = x[2]; X_3 = x[3]  
   if (X_2<0) X_2=0 
   with(as.list(vparameters),{
      npop = X_1+X_2+X_3  
      dX_1 = -beta*X_1*X_2/npop 
      dX_2 = +beta*X_1*X_2/npop - gamma*X_2 
      dX_3 = +gamma*X_2                 
      out = c(dX_1,dX_2,dX_3)
      list(out)
   })
}

Solving the system

gamma = 1/3; beta = 0.5
N = 1000; X2_0 = 10; X1_0 = N-X2_0; X3_0 = 0
vt = seq(0,1000,1)

vparameters=c(gamma=gamma,beta=beta)
inits=c(X_1=X1_0,X_2=X2_0,X_3=X3_0)
odemodel = as.data.frame(lsoda(inits, vt, ODEfunc, vparameters))
tail(odemodel)
##      time      X_1          X_2      X_3
## 996   995 406.3657 8.556423e-51 593.6343
## 997   996 406.3657 7.571048e-51 593.6343
## 998   997 406.3657 6.699150e-51 593.6343
## 999   998 406.3657 5.927662e-51 593.6343
## 1000  999 406.3657 5.245020e-51 593.6343
## 1001 1000 406.3657 4.640993e-51 593.6343

Solving the system

Markov Model

  • Each (random) component \(X_i\) in \(\mathbf{x}\) contains small changes in an infinitesimal time interval.

  • Recall Poisson approximation: the smallest change that can occur is one unit.

  • For a network: \[X_1 \longrightarrow X_2 \longrightarrow X_3\] The smallest change for \(X_i\) is plus or minus one.

  1. \(X_1\) can change by \(-1\), and \(X_2\) changes by \(+1\).
  2. \(X_2\) can change by \(-1\), and \(X_3\) changes by \(+1\).

State Variations

  • The number of possible change patterns is \(J=2\). \[\begin{array}{c} \Delta X_{1}\:\Delta X_{2}\:\Delta X_{3}\\ \lambda=\left(\begin{array}{ccc} -1 & 1 & 0\\ 0 & -1 & 1 \end{array}\right) \end{array}\]

  • Based on the current state of the system, determine the time step, \(\delta_t\), needed for just one unit of change \[\delta_t = \frac{1}{\beta X_1 X_2/N + \gamma X_2}\] which is one unit divided by the total state variation

Transition Probability

  • Based on the current state of the system, that each of the possible transitions will take place in time \(\delta_t\) such that \[\begin{array}{ccc} J & \Lambda_{j} & p_{j}/\delta_{t}\\ 1 & [-1,1,0] & \frac{\beta X_{1}X_{2}}{N}\\ 2 & [0,-1,1] & \gamma X_{2} \end{array}\]

  • If \(Y_1 \sim \mbox{Exp}(\theta_1)\) and \(Y_2 \sim \mbox{Exp}(\theta_2)\), then \(\min (Y_1, Y_2)\) is distributed as \(\mbox{Exp}(\theta_1+\theta_2)\). For a model with two possible state transitions, with expected rate \(\theta_1\) and \(\theta_2\), respectively, we expect the rate for at least one of the transitions to take place to be \(\theta_1+\theta_2\).

Poisson Approximation

  • Sample Poisson distributed random numbers based on these probabilities. \[ [-\Delta, \Delta, 0] \sim \mbox{Poi}(p_1 \times \delta_t) \\ [0, -\Delta, \Delta] \sim \mbox{Poi}(p_2 \times \delta_t)\] Note \(\mbox{Poi}=\mbox{Exp}\) when there is only \(1\) jump. Then \[ \mathbb{E} [ \Delta \mathbf{x}]= \mathbb{E}[-\Delta, \Delta, 0] + \mathbb{E}[0, -\Delta, \Delta] =(p_1 +p_2)\times \delta_t = 1.\]

  • Simulate the increment \(\Delta \mathbf{x}\) accordingly by the time step \(\delta_t\) such that \[\mathbf{x}(t) = \mathbf{x}(t-1) + \Delta \mathbf{x}.\]

Simulation

set.seed(20180420)
time = 0
vstate = c(X1_0,X2_0,X3_0)
K = length(vstate); J = 2              
lambda = matrix(0,nrow=J,ncol=length(vstate))
lambda[1,] = c(-1,1,0)
lambda[2,] = c(0,-1,1)
mcstate = list(); i = 1

Simulation

while(vstate[2]>0&vstate[1]>0){
  mcstate[[i]] = c(vstate,time)

  X_1 = vstate[1]; X_2 = vstate[2]; X_3 = vstate[3]
  vec_p = c(beta*X_1*X_2/N, gamma*X_2)
  delta_t = 1/sum(vec_p); vec_l = rpois(length(vec_p),vec_p*delta_t)
  
  vstate = vstate + vec_l%*%lambda
  vstate[vstate<0] = 0; i = i+1; time = time + delta_t
  if (i%%100==1) cat(i,time,vstate,"\n")}
## 101 8.126566 937 24 39 
## 201 13.10404 892 28 80 
## 301 17.0431 841 41 118 
## 401 20.35251 787 41 172 
## 501 23.89327 744 39 217 
## 601 27.16426 691 39 270 
## 701 31.13568 643 35 322 
## 801 34.93098 598 34 368 
## 901 42.43374 549 16 435 
## 1001 51.7809 499 14 487

Simulation

head(mcstate)
## [[1]]
## [1] 990  10   0   0
## 
## [[2]]
## [1] 990.0000000  10.0000000   0.0000000   0.1207243
## 
## [[3]]
## [1] 990.0000000  10.0000000   0.0000000   0.2414487
## 
## [[4]]
## [1] 990.000000  10.000000   0.000000   0.362173
## 
## [[5]]
## [1] 990.0000000  10.0000000   0.0000000   0.4828974
## 
## [[6]]
## [1] 990.0000000  10.0000000   0.0000000   0.6036217

Simulation

MCMC

  • The previous simulation is a low-cost computation.

  • Note that one realisation of the stochastic model is useless to infer anything.

  • When it comes to stochastic modelling, we need to run more simulations to make the model realistic.

MCMC

niter = 25  
mcstate = list()
i = 1
for (iter in 1:niter){
  ~COPY THE PREVIOUS SIMULATION CODE HERE~
  cat("Doing realisation:",iter,niter," ",time,vstate,"\n")
}
## Doing realisation: 1 25   60.01304 416 0 584 
## Doing realisation: 2 25   43.54139 425 0 576 
## Doing realisation: 3 25   58.62811 394 0 606 
## Doing realisation: 4 25   49.96172 503 0 498 
## Doing realisation: 5 25   49.52646 278 0 722 
## Doing realisation: 6 25   45.02408 338 0 662 
## Doing realisation: 7 25   45.14087 534 0 466 
## Doing realisation: 8 25   53.79226 462 0 538 
## Doing realisation: 9 25   76.38771 407 0 594 
## Doing realisation: 10 25   55.84398 398 0 602 
## Doing realisation: 11 25   59.78341 547 0 453 
## Doing realisation: 12 25   65.82722 396 0 604 
## Doing realisation: 13 25   60.54757 355 0 645 
## Doing realisation: 14 25   41.65695 477 0 523 
## Doing realisation: 15 25   59.01596 391 0 610 
## Doing realisation: 16 25   53.19685 543 0 457 
## Doing realisation: 17 25   49.07036 335 0 665 
## Doing realisation: 18 25   64.23467 399 0 601 
## Doing realisation: 19 25   54.79813 400 0 600 
## Doing realisation: 20 25   44.58744 401 0 599 
## Doing realisation: 21 25   75.84119 422 0 578 
## Doing realisation: 22 25   54.83471 504 0 496 
## Doing realisation: 23 25   81.24339 439 0 561 
## Doing realisation: 24 25   83.29638 464 0 536 
## Doing realisation: 25 25   88.35492 396 0 606

MCMC

Estimation

  • The parameter \(\beta\) and \(\gamma\) could be unknown. Unknown parameter vector: \(\theta= (\beta, \gamma)\)

  • Let \(Y(t)=X_2(t)\).

  • Observe \(y_t = (y_2, y_7, y_{12}, \dots y_{97})\), \(X_1(0)=990\) and \(X_3(0)=0\).

  • \(y_t\) can be thought of being Poisson distributed.

Estimation

  • Poisson maximum likelihood : the likelihood function is \[L(\alpha, y_1,\dots,y_n) = \Pi_{j=1}^{n} \exp(-\alpha)\frac{1}{y_j !} \alpha^{y_j}.\]

  • The log-likelihood function is \[\ell(\alpha, y_1,\dots,y_n) = -n \alpha - \sum_{j=1}{n} \ln (y_j !) + \ln (\alpha) \sum_{j=1}^{n}y_j\]

  • Negative log-likelihood for time varying parameter: \[-\ell(Y_1(\theta),\cdots, Y_n(\theta), y_1,\dots,y_n) = \sum_{j=1}^{n} Y_j (\theta) + \sum_{j=1}^{n} \ln (Y_j(\theta)) y_j.\]

Simulated Likelihood

  • The parameterized system of ODEs \[\frac{d\mathbf{x}(\theta, t)}{dt} = F(\mathbf{x}(\theta, t)).\]

  • Set up an optimization problem \[\hat{\theta} = \arg \min_{\theta} \mathrm{d}( (y_2(\theta), y_7(\theta),\dots y_{97}(\theta)), (y_2, y_7, y_{12}, \dots y_{97}))\\ \mbox{where } (y_2(\theta), y_7(\theta),\dots y_{97}(\theta)) \mbox{ from } \frac{dy(\theta, t)}{dt} = F(y(\theta, t))\] where \(\mathrm{d}(\cdot,\cdot)\) is the distant function.

Maximum Likelihood

times = odemodel[seq(2,100,by=5),1]
data = odemodel[seq(2,100,by=5),3]

params = c('beta'=0.3,'gamma'=0.2)

yinit = ode(c(X1_0, X2_0, X3_0), times, ODEfunc, params, method='ode45')

par(mfrow=c(1,2))
plot(times, yinit[,3], type='l')
points(times, data)

Maximum Likelihood

Maximum Likelihood

ML=function(params,times,data){
  params = abs(params)
  ycurr = ode(c(X1_0, X2_0, X3_0), times, ODEfunc, params, method='ode45')
  y = ycurr[,3]
  NLL =  sum(y) - sum(data*log(y))
  # NLL = -sum(log(dpois(round(data),round(y))))  
  # ML using normally distributed measurement error
  # NLL = -sum(log(dnorm(data,y,0.1*mean(data)))) or
  # NLL = sum((y - data)^2)  
}
res = optim(params,fn=ML,times=times,data=data)
estimators = res$par
yest = ode(c(X1_0, X2_0, X3_0), times, ODEfunc, estimators, method='ode45')
plot(times, yest[,3], type='l')
points(times, data)

Maximum Likelihood

Stochastic Differential Equations

  • Stochastic differential equations (SDE) are based on the deterministic model.

  • Add stochastic terms onto the model ODE’s to simulate random noise in the system. \[\frac{d\mathbf{X}(t)}{dt} = F(\mathbf{X}(t)) + \mathbf{G} \mathbf{W}\] where \(\mathbf{G}\) is \(K\times J\) matrix, \(\mathbf{W}\) is a \(J\)-vector white-noise.

  • The possible state transitions of one-units is described by a matrix, \(\Lambda=(\Lambda_1, \Lambda_2)\).

Stochastic Terms

  • It turns out that the \(\mathbf{G}\) matrix is intimately related to \(\Lambda\) and the \(p\) vector \[\mathbf{G} = (\sum_{j=1}^{2}\Lambda_{j}\times\sqrt{p_j/\delta_t})^{T}.\]

  • Recall \(J=2\) \[\begin{array}{ccc} J & \Lambda_{j} & p_{j}/\delta_{t}\\ 1 & [-1,1,0] & \frac{\beta X_{1}X_{2}}{N}\\ 2 & [0,-1,1] & \gamma X_{2}. \end{array}\]

Stochastic Terms

The matrix \(\mathbf{G}\) is \[ \mathbf{G}^{T}= \left(\begin{array}{ccc} -\sqrt{\frac{\beta X_{1}X_{2}}{N}} & \sqrt{\frac{\beta X_{1}X_{2}}{N}} & 0\\ 0 & -\sqrt{\gamma X_{2}} & \sqrt{\gamma X_{2}} \end{array}\right).\]

SDE

The SDE is \[\frac{dX_1}{dt} = - \beta \frac{X_1 X_2}{N} - \sqrt{\frac{\beta X_{1}X_{2}}{N}} W_1 \\ \frac{dX_2}{dt} = \beta \frac{X_1 X_2}{N} + \sqrt{\frac{\beta X_{1}X_{2}}{N}} W_1 - \gamma X_2 - \sqrt{\gamma X_{2}} W_2 \\ \frac{dX_3}{dt} = \gamma X_2 + \sqrt{\gamma X_{2}}W_2.\]

Time Series

  • The equations for implementation are \[ \Delta X_1 = - \Delta t \beta \frac{X_1 X_2}{N} - \sqrt{\Delta t} \sqrt{\frac{\beta X_{1}X_{2}}{N}} W_1 \\ \Delta X_2 = \Delta t \beta \frac{X_1 X_2}{N} + \sqrt{\Delta t} \sqrt{\frac{\beta X_{1}X_{2}}{N}} W_1 \\- \Delta t \gamma X_2 - \sqrt{\Delta t} \sqrt{\gamma X_{2}} W_2 \\ \Delta X_3 = \Delta t \gamma X_2 + \sqrt{\Delta t} \sqrt{\gamma X_{2}}W_2.\]

SDE

niter = 50; ldynamic_time_step = 0; Delta_t = 0.1
mcstate = list(); i = 1
for (iter in 1:niter){
  ~ COPY THE SIMULATION SETUP ~
  while(vstate[2]>0&vstate[1]>0){
    X_1 = vstate[1]; X_2 = vstate[2]; X_3 = vstate[3]
    vec_p = c(beta*X_1*X_2/N,gamma*X_2)
    scale = 10
    if (ldynamic_time_step) Delta_t = scale/sum(vec_p)
    G = lambda
    for (irow in 1:nrow(G)){
      G[irow,] = G[irow,]*sqrt(vec_p[irow])
    }
    G = t(G)
    W = rnorm(ncol(G)) # continued

SDE

    #... from last slide
    delta_X1 = Delta_t*(-beta*X_1*X_2/N)+
      sqrt(Delta_t)*sum(G[1,]*W)
    delta_X2 = Delta_t*(beta*X_1*X_2/N-gamma*X_2)+
      sqrt(Delta_t)*sum(G[2,]*W)
    delta_X3 = Delta_t*(gamma*X_2)+
      sqrt(Delta_t)*sum(G[3,]*W)

    vstate[1] = vstate[1] + delta_X1 
    vstate[2] = vstate[2] + delta_X2 
    vstate[3] = vstate[3] + delta_X3
    vstate[vstate<0] = 0
    i = i+1
    time = time + Delta_t
  }
  cat("Doing realisation:",iter,niter," ",time,vstate,"\n")
}

SDE

## Doing realisation: 1 50   55.9 354.6329 0 645.9328 
## Doing realisation: 2 50   43.9 439.9795 0 560.0653 
## Doing realisation: 3 50   2.5 992.4058 0 7.608571 
## Doing realisation: 4 50   129.9 438.0895 0 561.9113 
## Doing realisation: 5 50   63.7 427.0659 0 574.1617 
## Doing realisation: 6 50   59.8 387.1143 0 612.8917 
## Doing realisation: 7 50   65.2 374.9314 0 625.7895 
## Doing realisation: 8 50   46.1 398.1762 0 603.5565 
## Doing realisation: 9 50   70.6 454.1321 0 545.8773 
## Doing realisation: 10 50   46.7 385.647 0 614.3637 
## Doing realisation: 11 50   89 400.8219 0 599.9955 
## Doing realisation: 12 50   49.7 440.2219 0 559.8129 
## Doing realisation: 13 50   70.2 422.3841 0 577.6187 
## Doing realisation: 14 50   54.4 332.6087 0 667.8057 
## Doing realisation: 15 50   60.5 365.3872 0 634.6638 
## Doing realisation: 16 50   56.7 395.6705 0 604.34 
## Doing realisation: 17 50   53.6 380.6143 0 619.3871 
## Doing realisation: 18 50   62.1 408.9576 0 591.0432 
## Doing realisation: 19 50   62 382.5405 0 617.7588 
## Doing realisation: 20 50   55.3 330.7728 0 669.5614 
## Doing realisation: 21 50   55.3 523.3571 0 476.654 
## Doing realisation: 22 50   56.8 359.8165 0 640.1912 
## Doing realisation: 23 50   54 375.7596 0 625.764 
## Doing realisation: 24 50   64.1 402.4108 0 597.6158 
## Doing realisation: 25 50   50.4 537.1013 0 463.8532 
## Doing realisation: 26 50   66.6 367.0908 0 632.9147 
## Doing realisation: 27 50   58.4 364.9871 0 635.0362 
## Doing realisation: 28 50   69.1 406.9527 0 593.0497 
## Doing realisation: 29 50   76.7 370.4262 0 629.5811 
## Doing realisation: 30 50   68.9 390.4173 0 609.583 
## Doing realisation: 31 50   50.5 395.7199 0 604.3316 
## Doing realisation: 32 50   55 322.0314 0 678.405 
## Doing realisation: 33 50   57.1 428.6268 0 571.6581 
## Doing realisation: 34 50   66.3 507.8484 0 492.1822 
## Doing realisation: 35 50   67.3 416.5946 0 583.4577 
## Doing realisation: 36 50   50.6 470.9972 0 529.0682 
## Doing realisation: 37 50   77.9 360.1881 0 639.8548 
## Doing realisation: 38 50   57.4 398.7668 0 601.2376 
## Doing realisation: 39 50   63.8 401.117 0 600.3609 
## Doing realisation: 40 50   84.1 400.7222 0 599.6636 
## Doing realisation: 41 50   59.3 418.8213 0 581.193 
## Doing realisation: 42 50   50.6 397.7451 0 602.3121 
## Doing realisation: 43 50   69.3 524.6143 0 475.3943 
## Doing realisation: 44 50   55 403.319 0 596.7641 
## Doing realisation: 45 50   61.3 441.0079 0 558.9935 
## Doing realisation: 46 50   44.6 343.3439 0 657.5541 
## Doing realisation: 47 50   61.4 482.938 0 517.0656 
## Doing realisation: 48 50   63.3 360.2089 0 639.8183 
## Doing realisation: 49 50   64.3 394.2126 0 605.8488 
## Doing realisation: 50 50   44.3 352.4943 0 648.2694

SDE