Estimating Parameters for Chaotic Systems (Lotka-Volterra Redux)

[idontgetoutmuch]

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

@betanalpha

With more than two populations, especially for long integration times, the ODE integration is unlikely to properly calculate the chaotic trajectories and it will be hard to say much of anything that depends on the precise details of the solution at those later times. You can try increasing the absolute and relative precision of the integration, but you’re fighting exponential errors so it’s a losing game. Chaotic systems are hard

@Bob_Carpenter

It’s the ODE, not the MCMC that’s causing the problem here. So SMC isn’t going to help if you still have to solve a tricky ODE.

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

[bbbales2]

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.

[idontgetoutmuch]

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.

@bbbales2

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

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.

[idontgetoutmuch]

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)

[bbbales2]

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?

@idontgetoutmuch

No because

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.

[idontgetoutmuch]

@bbbales2

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

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

@bbbales2

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.

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

@idontgetoutmuch

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

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.

@bbbales2

Mind giving that a try? I’m curious what results you’d get.

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

[bbbales2]

@idontgetoutmuch

I don’t think that’s very loose

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?

@idontgetoutmuch

I will certainly give it a try but first I’d like to get the current model

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

@idontgetoutmuch

but my real problems are stiff

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

[idontgetoutmuch]

@bbbales2

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?

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

@bbbales2

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

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.

@bbbales2

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

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 😀

[bbbales2]

@idontgetoutmuch

Just trying to make things as simple as possible.

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

@idontgetoutmuch

Agreed but it’s not chaotic.

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.

[Charles_Driver]

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!

https://gist.github.com/cdriveraus/dfa3f25fae9e66817a27425f5950c3ad

ldat.csv (21.4 KB)

[idontgetoutmuch]

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?).

[bbbales2]

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!

[Charles_Driver]

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

[idontgetoutmuch]

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

[Bob_Carpenter]

@idontgetoutmuch

as I get errors (see below)

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.

@Charles_Driver

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

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