Looking for help of state space model and question about ode_rk45

Modeling
3

Hi Caesoma,

Thank you so much for your reply. They are all super helpful. I kept reading some related posts and papers these day. I understood more about the application and mechanism of the ode functions. Just would like to clarify my understanding. Please correct me if I am wrong.

The ode function take the targeted value of a time point t(i) and calculate the difference to the next one t(i+1). And the ode can estimate the parameter supplied in the function to correctly predict/estimate the targeted value at t(n) based on t(0), steps (n), and parameter. The parameters or interventions should also be a continuous process.

Blockquote
The way this is structured you seem to be simply using a parameter mu and fitting its values to the data, multiplying it by the growth/death rates (times c ) is not generating the effect you are probably expecting, it’s just scaling mu , but their values are probably meaningless.

Now I understand what you meant about the scaling mu with c. Thank you very much for pointed it out.

I also found your post about a SSM with a discrete output. I will try to understand that one as well.

Here is the updated model that separated the continuous process (handle by ode) and discrete disturbance event. This is just another try to make it closer to the cause of population fluctuation.

functions {
  vector dpop_dt(real t,               // time
                vector y, // popoulation size used in diffrential equation
                int N,
                // parameters
                real cp    // change rate = grwoth - death
                ) {
    // Declare differential equation type
    vector[1] dp_dt;
    
    // differential equations
    dp_dt[1] = cp * y[1]; 

    // Function product
    return dp_dt ;
  }
}
data {
  int<lower=0> N;      // number of measurement times
  array[N] real<lower=0> p;  // measured populations
  vector[N] c;
  array[N] real time;  // time steps
}
transformed data{
  real t0 = 0;
  array[N-1] real t;
  t = time[2:N];
}

parameters {
  real<lower=0> cp;      //  cp
  real<lower=0> sigma;      // measurement errors
  vector[1] y0;             // Initial population size
  real<lower=0> effectc;    //effect of cyclone
}
transformed parameters {
  // declare the estimated population size: pop
  // use ode_rk45 to estimate population size based on differential equation
  array[N] vector[1] pop;
  pop[1, 1] = y0[1];
  pop[2:N]= ode_rk45(dpop_dt, y0, t0, t, N, cp);
}
model {
  // priors
  cp ~ normal( 1 , 0.5 );      // gp
  sigma ~ exponential( 1 );    // measurment error
  y0 ~ lognormal( log(10) , 1 );    // initial population size
  effectc ~ lognormal(log(1), 0.5);  // effect of cyclone

  // observation model
  // connect latent population state to observed population size
  p[1] ~ lognormal(log(pop[1, 1]), sigma);
  for(i in 2:N){
    p[i] ~ lognormal(log(pop[i, 1] - effectc * c[i - 1] * pop[i, 1]) , sigma); // At the end of each time point, the effect of disturbance is introduced and substract a certain proportion of the population size
  }
}
generated quantities {
  // Directly produce posterior prediction in stan
  array[N] real p_pred;
  p_pred[1] = lognormal_rng(log(pop[1, 1]), sigma);
  for(i in 2:N){
    p_pred[i] = lognormal_rng(log(pop[i, 1] - effectc * c[i] *pop[i, 1),sigma);
  }
}