Estimating Parameters for Chaotic Systems (Lotka-Volterra Redux)

Modeling
4

Here’s my first attempt although I don’t think it’s correct as I get errors (see below) and some of the chains take a very long time (if fact I killed them off).

functions {
  real[] ds_dt(real t,       // time
               real[] s,     // system state
               real[] theta, // parameters
               real[] x_r,   // data
               int[] x_i)    // unused integer data
  {
    real x = s[1];
    real y = s[2];
    real z = s[3];

    real alpha = theta[1];

    real rho = x_r[1];
    real beta = x_r[2];

    real dx_dt = alpha * (y - x);
    real dy_dt = x * (rho - z) - y;
    real dz_dt = x * y - beta * z;

    return { dx_dt, dy_dt, dz_dt };
  }
}

data {
  int<lower = 1> N; // Number of measurements
  real ts[N];       // Times of the measurements
  real y[N];        // The measurements themselves
}

transformed data {
  real rho  = 45.92;
  real beta =  4.00;
  real rdummy[2];
  rdummy[1] = rho;
  rdummy[2] = beta;
}

parameters {
  real x[N,3];
  real ln_alpha[N];
  real mu_alpha;
  real<lower = 0> sigma_alpha;
}

model {
  real parms[1];
  real deltaT;
  x[1,] ~ lognormal(log(1.0), 0.2);
  ln_alpha[1] ~ normal(mu_alpha, sigma_alpha);
  for (i in 2:N) {
    parms[1] = exp(ln_alpha[i-1]);
    target += normal_lpdf(x[i,] | (integrate_ode_rk45(ds_dt, x[i-1,], ts[i-1], ts[i:i], parms, rdummy, rep_array(0, 0)))[1], 1.0E-3);
    deltaT = ts[i] - ts[i-1];
    target += normal_lpdf(ln_alpha[i] | ln_alpha[i - 1] - deltaT * sigma_alpha * sigma_alpha / 2, sqrt(deltaT) * sigma_alpha);
  }
  y ~ normal(x[,1], 0.2);
}

Here’s the data:

stan_synthetic.csv (9.5 KB)

I get this output from

model <- stan_model("LorenzStan.stan")
read_df = read.csv("stan_synthetic.csv")
stan_data <- list(N=length(read_df$X.time), ts=read_df$X.time, y=read_df$X_obs.value)
fit <- sampling(model, data = stan_data, iter = 400, cores = 4, chains = 4, seed = 123)

SAMPLING FOR MODEL 'LorenzStan' NOW (CHAIN 1).
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[3] is -0.150929, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.00966144, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.338497, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -0.249718, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[3] is -1.051, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.634166, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.156463, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -1.75997, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.229222, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -0.849689, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -1.28355, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -1.70489, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -1.7801, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -1.58441, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -0.000635361, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -1.47918, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[3] is -0.675783, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -1.70179, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.84971, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -0.453571, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)


Gradient evaluation took 0.011652 seconds
1000 transitions using 10 leapfrog steps per transition would take 116.52 seconds.
Adjust your expectations accordingly!


Iteration:   1 / 400 [  0%]  (Warmup)

SAMPLING FOR MODEL 'LorenzStan' NOW (CHAIN 2).
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.421135, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[3] is -0.936374, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -1.65254, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -1.49812, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -1.47934, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.592281, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -1.31991, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)


Gradient evaluation took 0.015128 seconds
1000 transitions using 10 leapfrog steps per transition would take 151.28 seconds.
Adjust your expectations accordingly!


Iteration:   1 / 400 [  0%]  (Warmup)

SAMPLING FOR MODEL 'LorenzStan' NOW (CHAIN 3).
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -0.809475, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)


Gradient evaluation took 0.011312 seconds
1000 transitions using 10 leapfrog steps per transition would take 113.12 seconds.
Adjust your expectations accordingly!


Iteration:   1 / 400 [  0%]  (Warmup)

SAMPLING FOR MODEL 'LorenzStan' NOW (CHAIN 4).
Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[3] is -1.03763, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.389829, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.601032, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[3] is -0.336315, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.122167, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -1.81669, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[3] is -1.48924, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.734302, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -1.8106, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -1.22877, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -1.20762, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -1.94996, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[3] is -1.85899, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -1.86504, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.609341, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -1.2936, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -1.39079, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.727495, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -0.626018, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -1.90859, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.651495, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[2] is -0.405731, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[3] is -0.732877, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[3] is -0.0372504, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: lognormal_lpdf: Random variable[1] is -0.0707845, but must be >= 0!  (in 'modela6b82e60a007_LorenzStan' at line 49)


Gradient evaluation took 0.010336 seconds
1000 transitions using 10 leapfrog steps per transition would take 103.36 seconds.
Adjust your expectations accordingly!


Iteration:  40 / 400 [ 10%]  (Warmup)
Iteration:   1 / 400 [  0%]  (Warmup)
Iteration:  40 / 400 [ 10%]  (Warmup)
Iteration:  40 / 400 [ 10%]  (Warmup)
Iteration:  40 / 400 [ 10%]  (Warmup)
Iteration:  80 / 400 [ 20%]  (Warmup)
Iteration: 120 / 400 [ 30%]  (Warmup)
Iteration:  80 / 400 [ 20%]  (Warmup)
Iteration: 160 / 400 [ 40%]  (Warmup)
Iteration:  80 / 400 [ 20%]  (Warmup)
Iteration: 200 / 400 [ 50%]  (Warmup)
Iteration: 201 / 400 [ 50%]  (Sampling)
Iteration:  80 / 400 [ 20%]  (Warmup)
Iteration: 120 / 400 [ 30%]  (Warmup)
Iteration: 240 / 400 [ 60%]  (Sampling)
Iteration: 280 / 400 [ 70%]  (Sampling)
Iteration: 160 / 400 [ 40%]  (Warmup)
Iteration: 200 / 400 [ 50%]  (Warmup)
Iteration: 201 / 400 [ 50%]  (Sampling)
Iteration: 320 / 400 [ 80%]  (Sampling)
Iteration: 240 / 400 [ 60%]  (Sampling)
Iteration: 360 / 400 [ 90%]  (Sampling)
Iteration: 400 / 400 [100%]  (Sampling)

 Elapsed Time: 9.46237 seconds (Warm-up)
               22.6419 seconds (Sampling)
               32.1042 seconds (Total)

Iteration: 280 / 400 [ 70%]  (Sampling)
Iteration: 320 / 400 [ 80%]  (Sampling)
Iteration: 360 / 400 [ 90%]  (Sampling)
Iteration: 400 / 400 [100%]  (Sampling)

 Elapsed Time: 17.4205 seconds (Warm-up)
               71.1707 seconds (Sampling)
               88.5913 seconds (Total)

Iteration: 120 / 400 [ 30%]  (Warmup)
Iteration: 120 / 400 [ 30%]  (Warmup)
Iteration: 160 / 400 [ 40%]  (Warmup)
Iteration: 160 / 400 [ 40%]  (Warmup)
Iteration: 200 / 400 [ 50%]  (Warmup)
Iteration: 200 / 400 [ 50%]  (Warmup)
Iteration: 201 / 400 [ 50%]  (Sampling)
Iteration: 201 / 400 [ 50%]  (Sampling)