ODE Event/Forcing Examples

knitr/ode-examples/data/two_comp.data.R

@@ -0,0 +1,13 @@
+K <- 2
+S <- 2
+T <- 91
+t0 <- 0.9
+theta <- 
+c(6, 0.6)
+ts <- 
+c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4,
+3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1,
+6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7,
+8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10)
+y0 <- 
+c(3, 6)

knitr/ode-examples/data/two_comp_e.data.R

@@ -0,0 +1,20 @@
+E <- 9
+end <- 
+c(11, 21, 31, 41, 51, 61, 71, 81, 91)
+K <- 2
+S <- 2
+start <- 
+c(1, 12, 22, 32, 42, 52, 62, 72, 82)
+T <- 91
+t0 <- 
+c(0.9, 2, 3, 4, 5, 6, 7, 8, 9)
+theta <- 
+c(6, 0.6)
+ts <- 
+c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4,
+3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1,
+6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7,
+8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10)
+y0 <- 
+structure(c(0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0),
+.Dim = c(9, 2))

knitr/ode-examples/data/two_comp_ef.data.R

@@ -0,0 +1,20 @@
+E <- 9
+end <- 
+c(11, 21, 31, 41, 51, 61, 71, 81, 91)
+K <- 2
+S <- 2
+start <- 
+c(1, 12, 22, 32, 42, 52, 62, 72, 82)
+T <- 91
+t0 <- 
+c(0.9, 2, 3, 4, 5, 6, 7, 8, 9)
+theta <- 
+c(6, 0.6)
+ts <- 
+c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4,
+3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1,
+6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7,
+8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10)
+y0 <- 
+structure(c(0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0),
+.Dim = c(9, 2))

knitr/ode-examples/data/two_comp_ef_fit.data.R

@@ -0,0 +1,50 @@
+E <- 9
+end <- 
+c(11, 21, 31, 41, 51, 61, 71, 81, 91)
+K <- 2
+S <- 2
+start <- 
+c(1, 12, 22, 32, 42, 52, 62, 72, 82)
+T <- 91
+t0 <- 
+c(0.9, 2, 3, 4, 5, 6, 7, 8, 9)
+theta <- 
+c(6, 0.6)
+ts <- 
+c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2,
+3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7,
+5.8, 5.9, 6, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2,
+8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10)
+y <- 
+structure(c(0.0501581870130331, 0.171390081674708, 0.154594268549219, 0.062065259109517, 0.0456829424379411, 0,
+0.208754305578589, 0, 0.0789583981577013, 0, 0, 1.13180612537173, 0.488747295719004, 0.175938234736721,
+0.267220370371945, 0.0596273820158646, 0.114899295074526, 0.272005811828668, 0.00321200420427842, 0.0389646648821001, 0,
+1.20316038685346, 0.904118171078214, 0.98223423254061, 0.565295590488372, 0.522927186207073, 0.450513348921578,
+0.620365646294086, 0.4678699568707, 0.37057962028136, 0.555095261235546, 1.74251402715741, 1.14647265763887,
+0.908613731853189, 0.784366309313526, 0.679373169916111, 0.417099037784891, 0.653246649036365, 0.675125974398131,
+0.716075219103793, 0.278019100183284, 1.35021531291007, 0.759361286445196, 0.465848106411765, 0.303854355878931, 0,
+0.044578906469163, 0, 0, 0, 0, 1.08801243790365, 0.559402702889078, 0.365565008029735, 0.161785136396892,
+0.168557908358781, 0.0127705580341626, 0.04979082649047, 0.0881241905230884, 0, 0, 1.49806956362223, 0.884173967451606,
+0.86524922688302, 0.502676374667655, 0.562566837439883, 0.476620165921046, 0.519539509440848, 0.218662655660154,
+0.507152267513044, 0.513108897994652, 1.49996791821863, 1.17971436361413, 0.772735263706271, 0.630423972389203,
+0.601900948517316, 0.582514140732464, 0.377620333755696, 0.635375982811837, 0.490251481506544, 0.569144292781708,
+1.46665189789528, 0.651212389778649, 0.314880821297867, 0.133348004125639, 0.117135004199271, 0.0739408633489209, 0,
+0.0789415707251588, 0.0717642520375456, 0, 0, 0, 0, 0, 0.0138889574457143, 0.0699245986925465, 0, 0.0317266826234932,
+0.0563205469874709, 0.051618034181338, 0, 0.841733488377438, 1.25733816758221, 1.29762986984564, 1.43653036429298,
+1.49264448556212, 1.47832947165913, 1.39076145170894, 1.44881690940123, 1.29443317998131, 1.1183347031059,
+2.05961490452341, 2.565206781713, 2.95672006361366, 3.17636228719576, 3.25413343388872, 3.5490595350775,
+3.51759017530972, 3.6994981311587, 3.57849480197177, 3.86912766768025, 4.8210520239222, 5.31095419533127,
+5.47286924999716, 5.63455729766036, 5.6922340334613, 5.66667261852888, 5.89035701971403, 5.69128413842573,
+5.5653802812262, 5.73401806024937, 6.25614629211007, 6.54052967860694, 6.40081961001448, 6.17014388043333,
+6.16265043518996, 5.46368006613957, 5.56019481866796, 5.20963902503991, 4.8578719579305, 4.650638197639,
+5.05538464031898, 5.1948775633154, 5.3505287015876, 5.1999262717203, 4.96108694353612, 4.7424579450011, 4.25759754692586,
+4.23333121655404, 4.02914711154008, 3.59490771430027, 4.50205709132745, 4.92517194331467, 5.12484869975415,
+5.12002491226594, 5.2107712706115, 5.06519828158445, 5.30513662266828, 5.21038342568496, 5.17183165177851,
+5.29201340365431, 6.1766443939395, 6.63712849299992, 6.49368187291553, 6.84264585701854, 6.71851118465657,
+6.71572385763204, 6.6640121268618, 6.43922293809072, 6.40336702821227, 6.10826396009771, 7.10335309546733,
+7.36477065392934, 7.0881024310677, 6.8154777231462, 6.61436329940234, 6.11434070682372, 5.93348503591031,
+5.46484119601089, 5.33824209235697, 5.20892459607444),
+.Dim = c(91, 2))
+y0 <- 
+structure(c(0, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0),
+.Dim = c(9, 2))

knitr/ode-examples/data/two_comp_f.data.R

@@ -0,0 +1,13 @@
+K <- 2
+S <- 2
+T <- 91
+t0 <- 0.9
+theta <- 
+c(6, 0.6)
+ts <- 
+c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4,
+3.5, 3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6, 6.1,
+6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7,
+8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10)
+y0 <- 
+c(3, 6)

knitr/ode-examples/models/one_comp.stan

@@ -0,0 +1,39 @@
+functions {
+  real[] one_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
+    real dydt[1];
+    dydt[1] = -theta[1] * y[1];
+    return dydt;
+  }
+}
+data {
+  int<lower=1> T;
+  int<lower=1> E;
+  int<lower=1> S;
+  real y0[E,S];
+  real t0[E];
+  real ts[T];
+  real theta[1];
+  int start[E];
+  int end[E];
+}
+transformed data {
+  real x[0];
+  int x_int[0];
+}
+model {
+}
+generated quantities {
+  real y_hat[T,S];
+  // 1. stop the integrator at t
+  // 2. take the last value of y_hat (value at t-1) as the new y0 (initial state)
+  // 3. combine y0 with the event at t
+  // 4. restart the integrator
+  for (i in 1:E) {
+    real y_state[S];
+    if (i == 1)
+      y_state = y0[1,];
+    else
+      y_state = to_array_1d(to_vector(y0[i,]) + to_vector(y_hat[end[i-1],]));
+    y_hat[start[i]:end[i],] = integrate_ode_bdf(one_comp, y_state, t0[i], ts[start[i]:end[i]], theta, x, x_int);
+  }
+}

knitr/ode-examples/models/two_comp.stan

@@ -0,0 +1,30 @@
+// two compartment (sim)
+functions {
+  real[] two_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
+    real dydt[2];
+    real a = theta[1];
+    real b = theta[2];
+    dydt[1] = -a*y[1];
+    dydt[2] = a*y[1] - b*y[2];
+    return dydt;
+  }
+}
+data {
+  int<lower=1> T;
+  int<lower=1> S;
+  int<lower=1> K;
+  real y0[S];
+  real t0;
+  real ts[T];
+  real theta[K];
+}
+transformed data {
+  real x[0];
+  int x_int[0];
+}
+model {
+}
+generated quantities {
+  real y_hat[T,S];
+  y_hat = integrate_ode_bdf(two_comp, y0, t0, ts, theta, x, x_int);
+}

knitr/ode-examples/models/two_comp_e.stan

@@ -0,0 +1,42 @@
+// two compartment event (sim)
+functions {
+  real[] two_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
+    real dydt[2];
+    dydt[1] = -theta[1] * y[1];
+    dydt[2] = theta[1] * y[1] - theta[2] * y[2];
+    return dydt;
+  }
+}
+data {
+  int<lower=1> T;
+  int<lower=1> E;
+  int<lower=1> S;
+  int<lower=1> K;
+  real y0[E,S];
+  real t0[E];
+  real ts[T];
+  real theta[K];
+  int start[E];
+  int end[E];
+}
+transformed data {
+  real x[0];
+  int x_int[0];
+}
+model {
+}
+generated quantities {
+  real y_hat[T,S];
+  // 1. stop the integrator at t
+  // 2. take the last value of y_hat (value at t-1) as the new y0 (initial state)
+  // 3. combine y0 with the event at t
+  // 4. restart the integrator
+  for (i in 1:E) {
+    real y_state[S];
+    if (i == 1)
+      y_state = y0[1,];
+    else
+      y_state = to_array_1d(to_vector(y0[i,]) + to_vector(y_hat[end[i-1],]));
+    y_hat[start[i]:end[i],] = integrate_ode_bdf(two_comp, y_state, t0[i], ts[start[i]:end[i]], theta, x, x_int);
+  }
+}

knitr/ode-examples/models/two_comp_ef.stan

@@ -0,0 +1,64 @@
+// two compartment event and forcing (sim)
+functions {
+    int max_indx (real t, real[] x) {
+    int N = num_elements(x);
+    int indx = 1;
+    while (indx <= N && x[indx] <= t) {
+      indx = indx + 1;
+    }
+    indx = indx - 1;
+    return indx;
+  }
+
+  real force_constant (real t, real[] x_step, real[] y_step) {
+    int indx;
+    real value;
+    indx = max_indx(t, x_step);
+    if (indx == 0)
+      value = 0;
+    else
+      value = y_step[indx];
+    return value;
+  }
+
+  real[] two_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
+    real dydt[2];
+    real c = force_constant(t, {0.0,3,5,7,9}, {0.0,3,0,3,0});
+    dydt[1] = -theta[1] * y[1] + c;
+    dydt[2] = theta[1] * y[1] - theta[2] * y[2];
+    return dydt;
+  }
+}
+data {
+  int<lower=1> T;
+  int<lower=1> E;
+  int<lower=1> S;
+  int<lower=1> K;
+  real y0[E,S];
+  real t0[E];
+  real ts[T];
+  real theta[K];
+  int start[E];
+  int end[E];
+}
+transformed data {
+  real x[0];
+  int x_int[0];
+}
+model {
+}
+generated quantities {
+  real y_hat[T,S];
+  // 1. stop the integrator at t
+  // 2. take the last value of y_hat (value at t-1) as the new y0 (initial state)
+  // 3. combine y0 with the event at t
+  // 4. restart the integrator
+  for (i in 1:E) {
+    real y_state[S];
+    if (i == 1)
+      y_state = y0[1,];
+    else
+      y_state = to_array_1d(to_vector(y0[i,]) + to_vector(y_hat[end[i-1],]));
+    y_hat[start[i]:end[i],] = integrate_ode_bdf(two_comp, y_state, t0[i], ts[start[i]:end[i]], theta, x, x_int);
+  }
+}

knitr/ode-examples/models/two_comp_ef_fit.stan

@@ -0,0 +1,71 @@
+// two compartment event and forcing (fit)
+functions {
+    int max_indx (real t, real[] x) {
+    int N = num_elements(x);
+    int indx = 1;
+    while (indx <= N && x[indx] <= t) {
+      indx = indx + 1;
+    }
+    indx = indx - 1;
+    return indx;
+  }
+
+  real force_constant (real t, real[] x_step, real[] y_step) {
+    int indx;
+    real value;
+    indx = max_indx(t, x_step);
+    if (indx == 0)
+      value = 0;
+    else
+      value = y_step[indx];
+    return value;
+  }
+
+  real[] two_comp (real t, real[] y, real[] theta, real[] x_r, int[] x_i) {
+    real dydt[2];
+    real c = force_constant(t, {0.0,3,5,7,9}, {0.0,3,0,3,0});
+    dydt[1] = -theta[1] * y[1] + c;
+    dydt[2] = theta[1] * y[1] - theta[2] * y[2];
+    return dydt;
+  }
+}
+data {
+  int<lower=1> T;
+  int<lower=1> E;
+  int<lower=1> S;
+  int<lower=1> K;
+  real y0[E,S];
+  real t0[E];
+  real ts[T];
+  int start[E];
+  int end[E];
+  real y[T,S];
+}
+transformed data {
+  real x[0];
+  int x_int[0];
+}
+parameters {
+  real theta[K];
+  real<lower=0> sigma;
+}
+model {
+  real y_hat[T,S];
+  // 1. stop the integrator at t
+  // 2. take the last value of y_hat (value at t-1) as the new y0 (initial state)
+  // 3. combine y0 with the event at t
+  // 4. restart the integrator
+  for (i in 1:E) {
+    real y_state[S];
+    if (i == 1)
+      y_state = y0[1,];
+    else
+      y_state = to_array_1d(to_vector(y0[i,]) + to_vector(y_hat[end[i-1],]));
+    y_hat[start[i]:end[i],] = integrate_ode_bdf(two_comp, y_state, t0[i], ts[start[i]:end[i]], theta, x, x_int);
+  }
+  for (t in 1:T)
+    target+= normal_lpdf(y[t] | y_hat[t], sigma);
+  target+= normal_lpdf(theta[1] | 5, 1);
+  target+= normal_lpdf(theta[2] | 0, 1);
+  target+= normal_lpdf(sigma | 0, 0.1);
+}