(Logistic growth) ODE model in Stan getting stuck before sampling begins

[Ben Lambert]

Hello all,

I am trying to solve the following ODE in Stan (I know that analytic solutions exist for this ODE, but it's a stepping stone for more complex models I want to fit):

 real[] logistic_diff(real t,real[] y,real[] theta,real[] x_r,int[] x_i) {

    real dydt[1];

    

    dydt[1] <- theta[1] * y[1] * (1 - theta[2] * y[1]);

    return dydt;

    }

The above defines a logistic growth model. So it has a couple of parameters, and is only a system of a single variable. I assume a normal sampling model, with another parameter: sigma.

I can generate data for this model (attached) and have been trying to get Stan to fit the three parameters of the model (theta[1],theta[2],sigma), where the true values are (0.5,0.05,2). 

I am having an issue whereby Stan appears to get 'stuck' before it has begun to sample - it's really strange; it doesn't appear to be due to poor convergence, as it hasn't even outputted ``iteration 1/2000'' for example. It doesn't happen for every run, but only for about 1/4 where it seems to break my rstudio, and I am forced to kill the process manually.

I don't think this is due to slow convergence - I have 100 simulated datasets, and I know that in general this simple system is well-identified. 

Does anyone know what is going on? (The full model is below.)

Best,

Ben

functions {

    real[] logistic_diff(real t,real[] y,real[] theta,real[] x_r,int[] x_i) {

    real dydt[1];

    

    dydt[1] <- theta[1] * y[1] * (1 - theta[2] * y[1]);

    return dydt;

    }

}

data {

    int<lower=1> T;

    int<lower=0> N;

    real t0;

    real ts[T];

    matrix[T,N] y;

}

transformed data {

    real x_r[0];

    int x_i[0];

}

parameters {

      real<lower=0> theta[2];

      real<lower=0> sigma;

      real<lower=0> y0[1];

}

model {

  real y_hat[T,1];

  sigma ~ cauchy(0,1);

  theta ~ normal(0,2);

  y0 ~ normal(5,2);

  y_hat <- integrate_ode(logistic_diff, y0, t0, ts, theta, x_r, x_i);

  for (i in 1:N){

    for (t in 1:T) {

          y[t,i] ~ normal(y_hat[t,1],sigma);

      }

  }

}

 

[Daniel Lee]

Hi Ben,

Thanks for the example. I've been trying to hunt down problems in a more complex problem, so it's useful seeing a stripped down example. 

I'll try to find out what's going on this week. 

Daniel

  y_hat <- integrate_ode(logitic_diff, y0, t0, ts, theta, x_r, x_i);

  for (i in 1:N){

    for (t in 1:T) {

          y[t,i] ~ normal(y_hat[t,1],sigma);

      }

  }

}

 

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ben Lambert]

Hi Dan,

That's great - thanks! If you do manage to have a look before Tuesday; please let me know - I was wanting to include this example in a lecture I'm giving then on the Stan language!

If not, no worries, but would be good to find out what's causing the issue nonetheless.

Best,

Ben

[Daniel Lee]

P.S. There's no way that program compiles. It has a spelling error. 

--

[Ben Lambert]

Hi Dan,

Sorry - that was an error of when I copied it across (I changed the name of the function to be something more meaningful!) I've changed this in the above now, so should compile fine.

Best,

Ben

[Sebastian Weber]

Hi!

Stan did freeze the integrator quite often due to trying to solve for extreme parameter values. This is solved in the about to appear 2.10 version. So just wait a few more days and try the new version. Please let us know if this solves your issues.

Best,
Sebastian

[Ben Lambert]

Great - I have put upper limits on my parameters for now, and that seems to solve the problem. (I will remove these when the new release comes out.) Thanks very much! Best, Ben

[Daniel Lee]

Ben, I was able to recover the parameters fine with RStan v2.9.0. I'm guessing there's something going on with how you generated data. Can you post your script? I'll try to replicate and figure out what's going on.

Here's mine:

theta <- c(0.5, 0.05)
sigma <- 2
y0 <- 5

N <- 10
T <- 20

library(deSolve)
t0 <- 0.0
ts <- 1:T

logistic_diff <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dy <- theta_1 * y * (1 - theta_2 * y)
    list(c(dy))
  })
}
parameters <- c(theta_1 = theta[1], theta_2 = theta[2])
out <- ode(y = c(y = y0), times = c(t0, ts), func = logistic_diff, parms = parameters)

out <- out[-1, 2]

y <- matrix(out, T, N)
y <- y + matrix(rnorm(T * N, 0, sigma), T, N)

data <- list(T = T, N = N, t0 = t0, ts = ts, y = y)

library(rstan)
fit <- stan("ode_example.stan", data = data)

On Sun, Jun 19, 2016 at 9:46 AM, 'Ben Lambert' via Stan users mailing list <stan-...@googlegroups.com> wrote:

Great - I have put upper limits on my parameters for now, and that seems to solve the problem. (I will remove these when the new release comes out.) Thanks very much! Best, Ben

--

[Daniel Lee]

I take it back... I'm seeing the behavior on different random seeds. I'll keep digging.

[Daniel Lee]

As soon as I said that, I couldn't recreate the issue. I'm trying different seeds and not seeing it. Please let me know if you have seeds to reproduce it.

Daniel

[Ben Lambert]

Hi Dan,

Thanks for trying here - I actually used Stan to generate the data (I called this using N=100 iterations using 'rstan'). Let me know how it goes, or if you need clarification on anything.

Best,

Ben

Functions {

    real[] sho(real t,real[] y,real[] theta,real[] x_r,int[] x_i) {

    real dydt[1];

    

    dydt[1] <- theta[1] * y[1] * (1 - theta[2] * y[1]);

    return dydt;

    }

}

data {

    int<lower=1> T;

    real y0[1];

    real t0;

    real ts[T];

    real sigma;

    real a1[2];

    real a2[2];

}

transformed data {

    real x_r[0];

    int x_i[0];

}

model {

}

generated quantities {

  real theta[2];

  real y_hat[T,1];

  

  theta[1] <- gamma_rng(a1[1],a1[2]);

  theta[2] <- gamma_rng(a2[1],a2[2]);

    

    y_hat <- integrate_ode(sho, y0, t0, ts, theta, x_r, x_i);

// add measurement error

    for (t in 1:T) {

        y_hat[t,1] <- y_hat[t,1] + normal_rng(0,sigma);

    }

}

[Daniel Lee]

Mind attaching the data? it'll make life easier. 

--

[Ben Lambert]

No problem - find it attached. I just tried running the aforementioned ODE using 24 chains for 1 iteration only, and found that about 3 of them got caught, and stan essentially froze.

Best,

Ben

[Bram Thijssen]

Dear Ben and Daniel,

Did you manage to find a solution for this, other than putting limits on the parameters?

I'm running into the same problem when trying to infer the parameters of an ODE system with Stan, even with the latest version 2.11 - the sampling freezes before the first iteration. I've put uniform priors on the parameters to restrict their range to reasonable values, but that doesn't seem to help in my case, so I'm curious whether you have been able to solve this?

Kind regards,
Bram

[Bob Carpenter]

Wold you mind sharing the model and data? Ideally
the simplest model and data that have the problem.

The new ODE solvers shouldn't go into infinite loops,
so if they do, it's a bug. It could be the sampler
taking many leapfrog steps for a slow solver.

You can put refresh=1 on the output to track every
iteration and see if it's really iteration 1 where things
are getting stuck.

The other culprit can be very small step sizes.

- Bob

[Bram Thijssen]

Hi Bob,

Of course, attached are the Stan and R files. The data was generated from the ODE model, using CVODE as well, but not using Stan.
The output is as follows - I've let it run for quite a while, but nothing happens anymore.

Thanks!
Bram

TRANSLATING MODEL 'tyson6var' FROM Stan CODE TO C++ CODE NOW.
successful in parsing the Stan model 'tyson6var'.
COMPILING THE C++ CODE FOR MODEL 'tyson6var' NOW.
OS: x86_64, mingw32; rstan: 2.11.1; Rcpp: 0.12.5; inline: 0.3.14 
Compilation argument:
 D:/PROGRA~1/R-32~1.4RE/bin/x64/R CMD SHLIB file89c4505188d.cpp 2> file89c4505188d.cpp.err.txt 
g++ -m64 -I"D:/PROGRA~1/R-32~1.4RE/include" -DNDEBUG   -I"D:/Programming/R-3.2.4revised/library/Rcpp/include/"  -I"D:/Programming/R-3.2.4revised/library/RcppEigen/include/"  -I"D:/Programming/R-3.2.4revised/library/RcppEigen/include/unsupported"  -I"D:/Programming/R-3.2.4revised/library/rstan/include/boost_not_in_BH" -I"D:/Programming/R-3.2.4revised/library/BH/include" -I"D:/Programming/R-3.2.4revised/library/StanHeaders/include/src/"  -I"D:/Programming/R-3.2.4revised/library/StanHeaders/include/"  -I"D:/Programming/R-3.2.4revised/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS    -I"d:/RCompile/r-compiling/local/local323/include"     -O2 -Wall  -mtune=core2 -c file89c4505188d.cpp -o file89c4505188d.o
g++ -m64 -shared -s -static-libgcc -o file89c4505188d.dll tmp.def file89c4505188d.o -LD:/Programming/R-3.2.4revised/library/StanHeaders/libs/x64 -lStanHeaders -Ld:/RCompile/r-compiling/local/local323/lib/x64 -Ld:/RCompile/r-compiling/local/local323/lib -LD:/PROGRA~1/R-32~1.4RE/bin/x64 -lR
D:/Programming/R-3.2.4revised/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function]

CHECKING DATA AND PREPROCESSING FOR MODEL 'tyson6var' NOW.

COMPILING MODEL 'tyson6var' NOW.

STARTING SAMPLER FOR MODEL 'tyson6var' NOW.

SAMPLING FOR MODEL 'tyson6var' NOW (CHAIN 1).

[Bob Carpenter]

[changed subject]

Sebastian --- mind taking a look at this to see why it's hanging?

- Bob

> <tyson6var.stan><runstan_6param.r>

[Bram Thijssen]

Hey Bob & Sebastian,

Did you have a chance to look at this?

I've noticed now that if I put limits on the parameters with the pointy bracket notation (e.g. real<lower=0.0015, upper=0.15> theta1 in the parameters block), instead of with a prior (e.g. theta1 ~ uniform(0.0015, 0.15) in the model block), then it does work normally. The sampler still has a lot of difficulty converging, but at least it's starting the sampling then. Does this give any clue what might be going on?

Bram

[Bob Carpenter]

> On Aug 22, 2016, at 10:26 AM, Bram Thijssen <brams...@gmail.com> wrote:
>
> Hey Bob & Sebastian,
>
> Did you have a chance to look at this?
>
> I've noticed now that if I put limits on the parameters with the pointy bracket notation (e.g. real<lower=0.0015, upper=0.15> theta1 in the parameters block), instead of with a prior (e.g. theta1 ~ uniform(0.0015, 0.15) in the model block), then it does work normally.

If you have theta1 ~ uniform(0.0015, 0.15) in the model block
then you are required to have the lower- and upper-bound
constraints real<lower=0.0015, upper=0.15>. The log density should
remain finite for any value of the parameters meeting the constraint.

So what's probably happening is that you're managing to initialize,
but then can't move because of rejections from parameter going
outside of (0.0015, 0.15).

> The sampler still has a lot of difficulty converging, but at least it's starting the sampling then. Does this give any clue what might be going on?

Yes :-)

Is there a natural constraint on theta1, such as theta1 in (0, 1)
or theta1 in (0, infinity)? If so, you may want to try a stronger
unconstrained prior, like normal(0.075, 0.07).

Another thing that will really help is if you can keep the parameters
on the unit scale---then Stan doesn't have to find the right scale
during warmup. So that'd be something like

theta_raw ~ normal(0, 1);

theta = 0.07 * theta_raw + 0.075;

It's an ugly hack, but it can be hard for Stan to adapt for complex
models during warmup.

- Bob