# Estimating Parameters for Chaotic Systems (Lotka-Volterra Redux)

In response to my question Hierarchical Lotka Volterra on estimating parameters, it was pointed out that estimating parameters for chaotic systems was hard

I have put together a blog post showing how I think it’s possible to model the parameters in a chaotic system (Lorenz) as variables undergoing a stochastic process (Brownian motion as it happens) and get apparently good estimates for the “parameters”.

Here are 3 paths of the estimated stochastic process for 3 different but close initial values:

And taking the average of e.g. the last 50 observations for each path gives me a good estimate of the actual parameter value.

I used libbi to model the system. Is it possible to do something similar in Stan (in particular model the parameter as a stochastic process)?

1 Like

This is cool!

Yeah, you can totally estimate this sorta model in Stan I think.

So what you have is a noise model for your measurements:

``````y ~ normal(x, sigma);
``````

And some sort of random transition model for your states (pseudo-coding it):

``````x[1] ~ normal(0, sigma_x);
for (i in 2:N) {
target += normal_lpdf(x[i] | ode(x[i - 1], param), sigma_x);
}
``````

And you also have some sort of state space progression going on with your parameters:

``````param[1] ~ normal(0, sigma_param);
for (i in 2:N) {
target += normal_lpdf(x[i] | param[i - 1], param_x);
}
``````

So you’ll need to estimate x (your state), param (the state of your parameters), sigma, sigma_x, sigma_param (your noise terms). I’d be curious to see the Stan results for this model.

Out of curiosity, did using a single parameter value not give you reliable estimates?

I was trying a similar thing for a model misspecification situation and it produced tantalizing results for me, but it was kindof a toy problem situation.

Thanks very much for this - I think this could be a better way to estimate “parameters” for differential equations (especially chaotic systems) - I hope to try this today.

No because (and here I might be displaying my ignorance) a small change in the parameter results in a very large change in behaviour so my thinking based on the above comments was that any sampler would really struggle. I guess I should try it.

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")
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)

1000 transitions using 10 leapfrog steps per transition would take 116.52 seconds.

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)

1000 transitions using 10 leapfrog steps per transition would take 151.28 seconds.

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)

1000 transitions using 10 leapfrog steps per transition would take 113.12 seconds.

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)

1000 transitions using 10 leapfrog steps per transition would take 103.36 seconds.

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)
``````

I’d have to look at the vectorization of the lpdfs more closely, but I think what you have is correct.

If you’re gonna put a really loose tolerance on the ODE solve like that, then it might be better to just do a fixed timestep forward Euler. Loose tolerance on the solvers can lead to problems with HMC cause I think the differentiability becomes a bit shaky (it’s variable timestep). The fixed timestep should stay differentiable, and if you don’t need much accuracy it can probly be okay since you’re integrating over short times? Maybe?

Hmm, so the reason I let my parameter vary is because I felt like at certain times the model I was using would be correct and at certain times it would be wrong. It was a linearized pendulum (x’’ = -a x) fit to data generated from a non-linear pendulum (x’’ = -a sin(x)). So there I let the parameter a vary cause at some times I’d expect the model to be good and sometimes bad.

But here you’re fitting the correct model the whole way through. I think fitting the ODE in pieces the way you’re doing is letting you get around the long time issues with this model. I think you’d expect the same parameter to hold for all your data though. Mind giving that a try? I’m curious what results you’d get.

Sadly I get errors which stop the model running - I’ll post these under a separate topic.

I’m not sure what you mean by loose tolerance here. Do you mean

I don’t think that’s very loose. I can do an explicit Euler in Stan but my real problems are stiff and I don’t really want to implement something like BDF in Stan. But perhaps you mean I am setting the time step for RK45 to a large value? I could use BDF then it will calculate its own (variable) time steps.

I will certainly give it a try but first I’d like to get the current model (parameters as stochastic variables) working and then we can compare the two for chaotic and non-chaotic systems (it would be an interesting blog post).

Whoops! I take it back. I thought the 1e-3 was the tolerance for the ODE solver. My mistake. Why not estimate the measurement error?

I also think the 1 parameter would be much easier to fit :P.

Eh, well this isn’t your real problem so there’s freedom to play :P?

Just trying to make things as simple as possible. Maybe it doesn’t help?

Agreed but it’s not chaotic. Ok so we can start with non-chaotic system. But then we shouldn’t need to model the parameter as a stochastic process; we can actually model it as a parameter. Maybe I can do a compare and contrast (parameter vs stochastic process) for the pendulum to start with. Can you post your code and data? And then I can have a go.

Indeed - I choose this because it’s well known for being chaotic but if you have another chaotic system you’d rather use, I am always happy to learn 😀

Eh, if it’s a known 1e-3, just leave it, it’s probly fine. Otherwise I’d estimate it.

I think it’s the fitting things in pieces with latent states that handles the chaotic-ness here. By resetting the ODE, then none of the integrations are long time, so there’s not the long time sensitivity issues that come with chaotic systems.

You asked in the lotka volterra thread about ctsem and I forgot to reply, this made for a nice test of some of the bits that have gone in lately. Rough around the edges but maybe helpful / interesting! I tried out the simple ODE form, the time varying form you propose (though with the dynamics estimated instead of just random walk), and an SDE formulation. With 100 data points I couldn’t get much success with any of them, though I didn’t a) let differential evolution run for very long (via steptol and itermax args), or b) try HMC. With 200 points, the two stochastic forms seem to work ‘fairly’ reliably (failure is very obvious!). The ODE form was hard, I think it worked once, but required measurement error fixed and it didn’t work when I just tried it :)

Oh – and the general structure here is that the latent states are integrated over using an unscented Kalman filter, so different to the implementation you have. Most situations where I’ve tried sampling states in Stan I’ve had a lot more trouble than with integration, but would be interested to get this form working also…

If you try it I’m curious how HMC goes with the various forms!

ldat.csv (21.4 KB)

Thanks for this - I am tied up with other things at the moment but I hope to look at this again next week in Helsinki (perhaps you will be there?).

Dunno if that “you” is directed at me or @Charles_Driver or both of us, but I’ll be around in Helsinki to yammer about this if you want. See you there!

I’m looking forward to Helsinki yes, see you both there :)

1 Like

Great - I am really looking forward to meeting you both :-)

1 Like

You need to declare `x` with a lower bound of zero if you’re going to give it a lognormal prior (for multiplicative noise scale). Why is that x lognormal and the other `x` noise scalees are additive (have just a normal)?

I couldn’t follow the use of `y` being the first column of `x[, 1]`, but presumably that’s right.

Much of the dev team will be there. It’s going to be a big crowd!