Added Stan model templates.

simo_IDCL_sindy_stan_code_12.stan

+functions {
+real[] dy_dt(real t, real[] y, real[] xi, real[] params, int[] ncupoly) {
+  vector[ncupoly[1]] poly;
+  row_vector[ncupoly[1]] T;
+  real ut;
+  real dydt;
+
+  poly = to_vector(params);
+
+  /* Excitation polynomial: excite system in-loop */
+  T[1] = 1;
+  T[2] = t;
+  T[3] = pow(t, 2);
+  T[4] = pow(t, 3);
+
+  ut = dot_product(poly, T);
+
+  /* Notice that `w=xi[1]` */
+  dydt = xi[2] * yt + xi[3] * ut + xi[1];
+  return { dydt };
+}
+}
+data {
+  int<lower=0> N;               // simulation horizon
+  real<lower=0> ts[N-1];        // Sampling instances
+  int<lower=0> nu;              // number of system inputs
+  int<lower=0> ny;              // number of system outputs
+  int<lower=0> nXi;             // number of SINDy model coefficients
+  real z_init[1, ny];           // initial model state = `y_init`
+  real y[N, ny];                // state observations
+  real u[N-1, nu];              // control inputs
+  real<lower=0> sigma_x_init;   // initial state associated uncertainty
+  real<lower=0> sigma_y_init;   // initial measurement associated uncertainty
+  real Xistar[nXi];             // SINDyc-optimal `Xi` (absolute values)
+  real l_Xibounds[nXi];         // Guess on the lower bounds of `Xi`
+  real u_Xibounds[nXi];         // Guess on the upper bounds of `Xi`
+  int<lower=0> ocupolys;        // Excitation polynomial degree
+  real params[nu, ocupolys];    // General parameter vector
+}
+parameters {
+  real<lower=l_Xibounds[1],upper=u_Xibounds[1]> Xi_1; // The Stan-optimal `Xi_1`
+  real<lower=l_Xibounds[2],upper=u_Xibounds[2]> Xi_2; // The Stan-optimal `Xi_2`
+  real w[1];                                          // State dynamics error
+  real<lower=0> sigma_Xi[nXi];
+  real<lower=0> sigma_w[1];     // Dynamics error variance
+  real<lower=0> sigma_y[N];     // Measurement error variance
+}
+transformed parameters {
+  real Xi[nXi];
+  real z[N, ny];
+  real Xi_[nXi+1];      // Augmented model (incl. dynamics error)
+  real results[N-1, 1];
+
+  Xi[1] = Xi_1;
+  Xi[2] = Xi_2;
+  Xi_ = append_array(w, Xi);
+
+  for (n in 1:ny){              /* Loop over outputs */
+    z[1, n] = z_init[1, n];
+    results = integrate_ode_rk45(dy_dt, z_init[, n], 0, ts, Xi_, params[n, ],
+                                 rep_array(ocupolys, 1));
+    z[2:, n] = results[, 1];
+  }
+}
+model {
+  /* Use for evaluation something along:
+     plt.hist(np.random.lognormal(mean=-2, sigma=1, size=100));plt.show() */
+  sigma_Xi[1] ~ lognormal(-1, .1);
+  sigma_Xi[2] ~ lognormal(-1, .1);
+  sigma_w ~ lognormal(log(sigma_x_init), .1); // Dynamics error variance prior
+  sigma_y ~ lognormal(log(sigma_y_init), .1); // Measurement error variance prior
+
+  Xi_1 ~ normal(Xistar[1], sigma_Xi[1]);
+  Xi_2 ~ normal(Xistar[2], sigma_Xi[2]);
+  w ~ normal(0., sigma_w);       // Dynamics error
+  for (n in 1:ny) {              // loop over outputs
+    y[, n] ~ normal(z[, n], sigma_y);
+  }
+}
+generated quantities {
+  real log_lik[N, ny];           // Check Log-Likelihood
+  real y_hat[N, ny];             // Model-based output predictions
+  for (n in 1:ny){               // Loop over trajectories
+    for (m in 1:N){              // Loop over timestamps
+      log_lik[m, n] = normal_lpdf(y[m, n] | z[m, n], sigma_y[m]);
+      y_hat[m, n] = normal_rng(z[m, n], sigma_y[m]);
+    }
+  }
+}

simo_IDCL_sindy_stan_code_1234.stan

+functions {
+real[] dy_dt(real t, real[] y, real[] xi, real[] params, int[] ncupoly) {
+  vector[ncupoly[1]] poly;
+  row_vector[ncupoly[1]] T;
+  real ut;
+  real dydt;
+
+  poly = to_vector(params);
+
+  /* Excitation polynomial: excite system in-loop */
+  T[1] = 1;
+  T[2] = t;
+  T[3] = pow(t, 2);
+  T[4] = pow(t, 3);
+
+  ut = dot_product(poly, T);
+
+  /* Notice that `w=xi[1]` */
+  dydt = xi[2] * y[1] + xi[3] * ut + xi[4] * pow(y[1], 2) + xi[5] * y[1] * ut + xi[1];
+  return { dydt };
+}
+}
+data {
+  int<lower=0> N;               // simulation horizon
+  real<lower=0> ts[N-1];        // Sampling instances
+  int<lower=0> nu;              // number of system inputs
+  int<lower=0> ny;              // number of system outputs
+  int<lower=0> nXi;             // number of SINDy model coefficients
+  real z_init[1, ny];           // initial model state = `y_init`
+  real y[N, ny];                // state observations
+  real u[N-1, nu];              // control inputs
+  real<lower=0> sigma_x_init;   // initial state associated uncertainty
+  real<lower=0> sigma_y_init;   // initial measurement associated uncertainty
+  real Xistar[nXi];             // SINDyc-optimal `Xi` (absolute values)
+  real l_Xibounds[nXi];         // Guess on the lower bounds of `Xi`
+  real u_Xibounds[nXi];         // Guess on the upper bounds of `Xi`
+  int<lower=0> ocupolys;        // Excitation polynomial degree
+  real params[nu, ocupolys];    // General parameter vector
+}
+parameters {
+  real<lower=l_Xibounds[1],upper=u_Xibounds[1]> Xi_1; // The Stan-optimal `Xi_1`
+  real<lower=l_Xibounds[2],upper=u_Xibounds[2]> Xi_2; // The Stan-optimal `Xi_2`
+  real<lower=l_Xibounds[3],upper=u_Xibounds[3]> Xi_3; // The Stan-optimal `Xi_3`
+  real<lower=l_Xibounds[4],upper=u_Xibounds[4]> Xi_4; // The Stan-optimal `Xi_4`
+  real w[1];                                          // State dynamics error
+  real<lower=0> sigma_Xi[nXi];
+  real<lower=0> sigma_w[1];     // Dynamics error variance
+  real<lower=0> sigma_y[N];     // Measurement error variance
+}
+transformed parameters {
+  real Xi[nXi];
+  real z[N, ny];
+  real Xi_[nXi+1];      // Augmented model (incl. dynamics error)
+  real results[N-1, 1];
+
+  Xi[1] = Xi_1;
+  Xi[2] = Xi_2;
+  Xi[3] = Xi_3;
+  Xi[4] = Xi_4;
+  Xi_ = append_array(w, Xi);
+
+  for (n in 1:ny){              /* Loop over outputs */
+    z[1, n] = z_init[1, n];
+    results = integrate_ode_rk45(dy_dt, z_init[, n], 0, ts, Xi_, params[n, ],
+                                 rep_array(ocupolys, 1));
+    z[2:, n] = results[, 1];
+  }
+}
+model {
+  /* Use for evaluation something along:
+     plt.hist(np.random.lognormal(mean=-2, sigma=1, size=100));plt.show() */
+  sigma_Xi[1] ~ lognormal(-1, .1);
+  sigma_Xi[2] ~ lognormal(-1, .1);
+  sigma_Xi[3] ~ lognormal(-1, .1);
+  sigma_Xi[4] ~ lognormal(-1, .1);
+  sigma_w ~ lognormal(log(sigma_x_init), .1); // Dynamics error variance prior
+  sigma_y ~ lognormal(log(sigma_y_init), .1);   // Inform error scale
+
+  Xi_1 ~ normal(Xistar[1], sigma_Xi[1]);
+  Xi_2 ~ normal(Xistar[2], sigma_Xi[2]);
+  Xi_3 ~ normal(Xistar[3], sigma_Xi[3]);
+  Xi_4 ~ normal(Xistar[4], sigma_Xi[4]);
+
+  w ~ normal(0., sigma_w);       // Dynamics error
+  for (n in 1:ny) {              // loop over outputs
+    y[, n] ~ normal(z[, n], sigma_y);
+  }
+}
+generated quantities {
+  real log_lik[N, ny];           // Check Log-Likelihood
+  real y_hat[N, ny];             // Model-based output predictions
+  for (n in 1:ny){               // Loop over trajectories
+    for (m in 1:N){              // Loop over timestamps
+      log_lik[m, n] = normal_lpdf(y[m, n] | z[m, n], sigma_y[m]);
+      y_hat[m, n] = normal_rng(z[m, n], sigma_y[m]);
+    }
+  }
+}

simo_IDCL_sindy_stan_code_124.stan

+functions {
+real[] dy_dt(real t, real[] y, real[] xi, real[] params, int[] ncupoly) {
+  vector[ncupoly[1]] poly;
+  row_vector[ncupoly[1]] T;
+  real ut;
+  real dydt;
+
+  poly = to_vector(params);
+
+  /* Excitation polynomial: excite system in-loop */
+  T[1] = 1;
+  T[2] = t;
+  T[3] = pow(t, 2);
+  T[4] = pow(t, 3);
+
+  ut = dot_product(poly, T);
+
+  /* Notice that `w=xi[1]` */
+  dydt = xi[2] * y[1] + xi[3] * ut + xi[4] * y[1] * ut + xi[1];
+  return { dydt };
+}
+}
+data {
+  int<lower=0> N;               // simulation horizon
+  real<lower=0> ts[N-1];        // Sampling instances
+  int<lower=0> nu;              // number of system inputs
+  int<lower=0> ny;              // number of system outputs
+  int<lower=0> nXi;             // number of SINDy model coefficients
+  real z_init[1, ny];           // initial model state = `y_init`
+  real y[N, ny];                // state observations
+  real u[N-1, nu];              // control inputs
+  real<lower=0> sigma_x_init;   // initial state associated uncertainty
+  real<lower=0> sigma_y_init;   // initial measurement associated uncertainty
+  real Xistar[nXi];             // SINDyc-optimal `Xi` (absolute values)
+  real l_Xibounds[nXi];         // Guess on the lower bounds of `Xi`
+  real u_Xibounds[nXi];         // Guess on the upper bounds of `Xi`
+  int<lower=0> ocupolys;        // Excitation polynomial degree
+  real params[nu, ocupolys];    // General parameter vector
+}
+parameters {
+  real<lower=l_Xibounds[1],upper=u_Xibounds[1]> Xi_1; // The Stan-optimal `Xi_1`
+  real<lower=l_Xibounds[2],upper=u_Xibounds[2]> Xi_2; // The Stan-optimal `Xi_2`
+  real<lower=l_Xibounds[3],upper=u_Xibounds[3]> Xi_3; // The Stan-optimal `Xi_3`
+  real<lower=0> sigma_Xi[nXi];
+  real<lower=0> sigma_w[1];     // Dynamics error variance
+  real<lower=0> sigma_y[N];        // TODO error scale
+}
+transformed parameters {
+  real Xi[nXi];
+  real z[N, ny];
+  real Xi_[nXi+1];      // Augmented model (incl. dynamics error)
+  real results[N-1, 1];
+
+  Xi[1] = Xi_1;
+  Xi[2] = Xi_2;
+  Xi[3] = Xi_3;
+  Xi_ = append_array(w, Xi);
+
+  for (n in 1:ny){              /* Loop over outputs */
+    z[1, n] = z_init[1, n];
+    results = integrate_ode_rk45(dy_dt, z_init[, n], 0, ts, Xi_, params[n, ],
+                                 rep_array(ocupolys, 1));
+    z[2:, n] = results[, 1];
+  }
+}
+model {
+  /* Use for evaluation something along:
+     plt.hist(np.random.lognormal(mean=-2, sigma=1, size=100));plt.show() */
+  sigma_Xi[1] ~ lognormal(-1, .1);
+  sigma_Xi[2] ~ lognormal(-1, .1);
+  sigma_Xi[3] ~ lognormal(-1, .1);
+  sigma_w ~ lognormal(log(sigma_x_init), .1); // Dynamics error variance prior
+  sigma_y ~ lognormal(log(sigma_y_init), .1);   // Inform error scale
+
+  Xi_1 ~ normal(Xistar[1], sigma_Xi[1]);
+  Xi_2 ~ normal(Xistar[2], sigma_Xi[2]);
+  Xi_3 ~ normal(Xistar[3], sigma_Xi[3]);
+
+  w ~ normal(0., sigma_w);       // Dynamics error
+  for (n in 1:ny) {             // Loop over outputs
+    y[, n] ~ normal(z[, n], sigma_y);
+  }
+}
+generated quantities {
+  real log_lik[N, ny];           // Check Log-Likelihood
+  real y_hat[N, ny];             // Model-based output predictions
+  for (n in 1:ny){               // Loop over trajectories
+    for (m in 1:N){              // Loop over timestamps
+      log_lik[m, n] = normal_lpdf(y[m, n] | z[m, n], sigma_y[m]);
+      y_hat[m, n] = normal_rng(z[m, n], sigma_y[m]);
+    }
+  }
+}

simo_IDCL_sindy_stan_code_23.stan

+functions {
+real[] dy_dt(real t, real[] y, real[] xi, real[] params, int[] ncupoly) {
+  vector[ncupoly[1]] poly;
+  row_vector[ncupoly[1]] T;
+  real ut;
+  real dydt;
+
+  poly = to_vector(params);
+
+  /* Excitation polynomial: excite system in-loop */
+  T[1] = 1;
+  T[2] = t;
+  T[3] = pow(t, 2);
+  T[4] = pow(t, 3);
+
+  ut = dot_product(poly, T);
+
+  /* Notice that `w=xi[1]` */
+  dydt = xi[2] * ut + xi[3] * pow(y[1], 2) + xi[1];
+  return { dydt };
+}
+}
+data {
+  int<lower=0> N;               // simulation horizon
+  real<lower=0> ts[N-1];        // Sampling instances
+  int<lower=0> nu;              // number of system inputs
+  int<lower=0> ny;              // number of system outputs
+  int<lower=0> nXi;             // number of SINDy model coefficients
+  real z_init[1, ny];           // initial model state = `y_init`
+  real y[N, ny];                // state observations
+  real u[N-1, nu];              // control inputs
+  real<lower=0> sigma_x_init;   // initial state associated uncertainty
+  real<lower=0> sigma_y_init;   // initial measurement associated uncertainty
+  real Xistar[nXi];             // SINDyc-optimal `Xi` (absolute values)
+  real l_Xibounds[nXi];         // Guess on the lower bounds of `Xi`
+  real u_Xibounds[nXi];         // Guess on the upper bounds of `Xi`
+  int<lower=0> ocupolys;        // Excitation polynomial degree
+  real params[nu, ocupolys];    // General parameter vector
+}
+parameters {
+  real<lower=l_Xibounds[1],upper=u_Xibounds[1]> Xi_1; // The Stan-optimal `Xi_1`
+  real<lower=l_Xibounds[2],upper=u_Xibounds[2]> Xi_2; // The Stan-optimal `Xi_2`
+  real w[1];                                          // State dynamics error
+  real<lower=0> sigma_Xi[nXi];
+  real<lower=0> sigma_w[1];     // Dynamics error variance
+  real<lower=0> sigma_y[N];     // Measurement error variance
+}
+transformed parameters {
+  real Xi[nXi];
+  real z[N, ny];
+  real Xi_[nXi+1];      // Augmented model (incl. dynamics error)
+  real results[N-1, 1];
+
+  Xi[1] = Xi_1;
+  Xi[2] = Xi_2;
+  Xi_ = append_array(w, Xi);
+
+  for (n in 1:ny){              /* Loop over outputs */
+    z[1, n] = z_init[1, n];
+    results = integrate_ode_rk45(dy_dt, z_init[, n], 0, ts, Xi_, params[n, ],
+                                 rep_array(ocupolys, 1));
+    z[2:, n] = results[, 1];
+  }
+}
+model {
+  /* Use for evaluation something along:
+     plt.hist(np.random.lognormal(mean=-2, sigma=1, size=100));plt.show() */
+  sigma_Xi[1] ~ lognormal(-1, .1);
+  sigma_Xi[2] ~ lognormal(-1, .1);
+  sigma_w ~ lognormal(log(sigma_x_init), .1); // Dynamics error variance prior
+  sigma_y ~ lognormal(log(sigma_y_init), .1); // Measurement error variance prior
+
+  Xi_1 ~ normal(Xistar[1], sigma_Xi[1]);
+  Xi_2 ~ normal(Xistar[2], sigma_Xi[2]);
+  w ~ normal(0., sigma_w);       // Dynamics error
+  for (n in 1:ny) {              // loop over outputs
+    y[, n] ~ normal(z[, n], sigma_y);
+  }
+}
+generated quantities {
+  real log_lik[N, ny];           // Check Log-Likelihood
+  real y_hat[N, ny];             // Model-based output predictions
+  for (n in 1:ny){               // Loop over trajectories
+    for (m in 1:N){              // Loop over timestamps
+      log_lik[m, n] = normal_lpdf(y[m, n] | z[m, n], sigma_y[m]);
+      y_hat[m, n] = normal_rng(z[m, n], sigma_y[m]);
+    }
+  }
+}

simo_IDCL_sindy_stan_code_34.stan

+functions {
+real[] dy_dt(real t, real[] y, real[] xi, real[] params, int[] ncupoly) {
+  vector[ncupoly[1]] poly;
+  row_vector[ncupoly[1]] T;
+  real ut;
+  real dydt;
+
+  poly = to_vector(params);
+
+  /* Excitation polynomial: excite system in-loop */
+  T[1] = 1;
+  T[2] = t;
+  T[3] = pow(t, 2);
+  T[4] = pow(t, 3);
+
+  ut = dot_product(poly, T);
+
+  /* Notice that `w=xi[1]` */
+  dydt = xi[2] * pow(y[1], 2) + xi[3] * y[1] * ut + xi[1];
+  return { dydt };
+}
+}
+data {
+  int<lower=0> N;               // simulation horizon
+  real<lower=0> ts[N-1];        // Sampling instances
+  int<lower=0> nu;              // number of system inputs
+  int<lower=0> ny;              // number of system outputs
+  int<lower=0> nXi;             // number of SINDy model coefficients
+  real z_init[1, ny];           // initial model state = `y_init`
+  real y[N, ny];                // state observations
+  real u[N-1, nu];              // control inputs
+  real<lower=0> sigma_x_init;   // initial state associated uncertainty
+  real<lower=0> sigma_y_init;   // initial measurement associated uncertainty
+  real Xistar[nXi];             // SINDyc-optimal `Xi` (absolute values)
+  real l_Xibounds[nXi];         // Guess on the lower bounds of `Xi`
+  real u_Xibounds[nXi];         // Guess on the upper bounds of `Xi`
+  int<lower=0> ocupolys;        // Excitation polynomial degree
+  real params[nu, ocupolys];    // General parameter vector
+}
+parameters {
+  real<lower=l_Xibounds[1],upper=u_Xibounds[1]> Xi_1; // The Stan-optimal `Xi_1`
+  real<lower=l_Xibounds[2],upper=u_Xibounds[2]> Xi_2; // The Stan-optimal `Xi_2`
+  real w[1];                                          // State dynamics error
+  real<lower=0> sigma_Xi[nXi];
+  real<lower=0> sigma_w[1];     // Dynamics error variance
+  real<lower=0> sigma_y[N];     // Measurement error variance
+}
+transformed parameters {
+  real Xi[nXi];
+  real z[N, ny];
+  real Xi_[nXi+1];      // Augmented model (incl. dynamics error)
+  real results[N-1, 1];
+
+  Xi[1] = Xi_1;
+  Xi[2] = Xi_2;
+  Xi_ = append_array(w, Xi);
+
+  for (n in 1:ny){              /* Loop over outputs */
+    z[1, n] = z_init[1, n];
+    results = integrate_ode_rk45(dy_dt, z_init[, n], 0, ts, Xi_, params[n, ],
+                                 rep_array(ocupolys, 1));
+    z[2:, n] = results[, 1];
+  }
+}
+model {
+  /* Use for evaluation something along:
+     plt.hist(np.random.lognormal(mean=-2, sigma=1, size=100));plt.show() */
+  sigma_Xi[1] ~ lognormal(-1, .1);
+  sigma_Xi[2] ~ lognormal(-1, .1);
+  sigma_w ~ lognormal(log(sigma_x_init), .1); // Dynamics error variance prior
+  sigma_y ~ lognormal(log(sigma_y_init), .1); // Measurement error variance prior
+
+  Xi_1 ~ normal(Xistar[1], sigma_Xi[1]);
+  Xi_2 ~ normal(Xistar[2], sigma_Xi[2]);
+  w ~ normal(0., sigma_w);       // Dynamics error
+  for (n in 1:ny) {              // loop over outputs
+    y[, n] ~ normal(z[, n], sigma_y);
+  }
+}
+generated quantities {
+  real log_lik[N, ny];           // Check Log-Likelihood
+  real y_hat[N, ny];             // Model-based output predictions
+  for (n in 1:ny){               // Loop over trajectories
+    for (m in 1:N){              // Loop over timestamps
+      log_lik[m, n] = normal_lpdf(y[m, n] | z[m, n], sigma_y[m]);
+      y_hat[m, n] = normal_rng(z[m, n], sigma_y[m]);
+    }
+  }
+}

simo_IDCL_sindy_stan_code_345.stan

+functions {
+real[] dy_dt(real t, real[] y, real[] xi, real[] params, int[] ncupoly) {
+  vector[ncupoly[1]] poly;
+  row_vector[ncupoly[1]] T;
+  real ut;
+  real dydt;
+
+  poly = to_vector(params);
+
+  /* Excitation polynomial: excite system in-loop */
+  T[1] = 1;
+  T[2] = t;
+  T[3] = pow(t, 2);
+  T[4] = pow(t, 3);
+
+  ut = dot_product(poly, T);
+
+  /* Notice that `w=xi[1]` */
+  dydt = xi[2] * pow(y[1], 2) + xi[3] * y[1] * ut + xi[4] * pow(ut, 2) + xi[1];
+  return { dydt };
+}
+}
+data {
+  int<lower=0> N;               // simulation horizon
+  real<lower=0> ts[N-1];        // Sampling instances
+  int<lower=0> nu;              // number of system inputs
+  int<lower=0> ny;              // number of system outputs
+  int<lower=0> nXi;             // number of SINDy model coefficients
+  real z_init[1, ny];           // initial model state = `y_init`
+  real y[N, ny];                // state observations
+  real u[N-1, nu];              // control inputs
+  real<lower=0> sigma_x_init;   // initial state associated uncertainty
+  real<lower=0> sigma_y_init;   // initial measurement associated uncertainty
+  real Xistar[nXi];             // SINDyc-optimal `Xi` (absolute values)
+  real l_Xibounds[nXi];         // Guess on the lower bounds of `Xi`
+  real u_Xibounds[nXi];         // Guess on the upper bounds of `Xi`
+  int<lower=0> ocupolys;        // Excitation polynomial degree
+  real params[nu, ocupolys];    // General parameter vector
+}
+parameters {
+  real<lower=l_Xibounds[1],upper=u_Xibounds[1]> Xi_1; // The Stan-optimal `Xi_1`
+  real<lower=l_Xibounds[2],upper=u_Xibounds[2]> Xi_2; // The Stan-optimal `Xi_2`
+  real<lower=l_Xibounds[3],upper=u_Xibounds[3]> Xi_3; // The Stan-optimal `Xi_3`
+  real w[1];                                          // State dynamics error
+  real<lower=0> sigma_Xi[nXi];
+  real<lower=0> sigma_w[1];     // Dynamics error variance
+  real<lower=0> sigma_y[N];     // Measurement error variance
+}
+transformed parameters {
+  real Xi[nXi];
+  real z[N, ny];
+  real Xi_[nXi+1];      // Augmented model (incl. dynamics error)
+  real results[N-1, 1];
+
+  Xi[1] = Xi_1;
+  Xi[2] = Xi_2;
+  Xi[3] = Xi_3;
+  Xi_ = append_array(w, Xi);
+
+  for (n in 1:ny){              /* Loop over outputs */
+    z[1, n] = z_init[1, n];
+    results = integrate_ode_rk45(dy_dt, z_init[, n], 0, ts, Xi_, params[n, ],
+                                 rep_array(ocupolys, 1));
+    z[2:, n] = results[, 1];
+  }
+}
+model {
+  /* Use for evaluation something along:
+     plt.hist(np.random.lognormal(mean=-2, sigma=1, size=100));plt.show() */
+  sigma_Xi[1] ~ lognormal(-1, .1);
+  sigma_Xi[2] ~ lognormal(-1, .1);
+  sigma_Xi[3] ~ lognormal(-1, .1);
+  sigma_w ~ lognormal(log(sigma_x_init), .1); // Dynamics error variance prior
+  sigma_y ~ lognormal(log(sigma_y_init), .1);   // Inform error scale
+
+  Xi_1 ~ normal(Xistar[1], sigma_Xi[1]);
+  Xi_2 ~ normal(Xistar[2], sigma_Xi[2]);
+  Xi_3 ~ normal(Xistar[3], sigma_Xi[3]);
+  w ~ normal(0., sigma_w);       // Dynamics error
+  for (n in 1:ny) {             // Loop over outputs
+    y[, n] ~ normal(z[, n], sigma_y);
+  }
+}
+generated quantities {
+  real log_lik[N, ny];           // Check Log-Likelihood
+  real y_hat[N, ny];             // Model-based output predictions
+  for (n in 1:ny){               // Loop over trajectories
+    for (m in 1:N){              // Loop over timestamps
+      log_lik[m, n] = normal_lpdf(y[m, n] | z[m, n], sigma_y[m]);
+      y_hat[m, n] = normal_rng(z[m, n], sigma_y[m]);
+    }
+  }
+}