Non linear differential equation in stan

Modeling
1

Hello everyone,

I am trying to use Stan to model the behaviour of the physics pendulum and modelling the different source of friction. The ODE I am trying to use is:

equation

I am able to successfully simulate and recover the parameters when \gamma = \beta = 0. However, the solver is having a lot of trouble when including the constant friction terme \gamma. The main issue is that the sign is linked to the second derivative of \theta so it’s always opposing the motion.

So I guess the solver has a lot of trouble because of that discontinuity.

Here is the R code I used :


sho_fun_body <- "
  vector[2] dy_dt;
  dy_dt[1] = y[2];
  
  real  s  = - w * sin(y[1]);
        s += - alpha * y[2];
       
  if (y[2] > 0 )
        s += -1* gamma;
  else if  (y[2] < 0 )
        s += + gamma;
  else
        s += 0;
        

  dy_dt[2] = s ;

  return(dy_dt);
"
sho_loglik_body <- "
  real loglik = 0.0;
  for(n in 1:N) {
    loglik += normal_lpdf(y1_obs[n] | y_sol[n][1], sigma);
  }
  return(loglik);
"

N <- stan_dim("N", lower=1)
D <- stan_dim("D")
y0 <- stan_vector("y0", length = D)
y1_obs <- stan_vector("y1_obs", length=N)

w <- stan_param(stan_var("w", lower=0), "inv_gamma(5, 1)")
alpha <- stan_param(stan_var("alpha", lower=0), "inv_gamma(5, 1)")
gamma <- stan_param(stan_var("gamma", lower=0), "inv_gamma(5, 1)")

sigma <- stan_param(stan_var("sigma", lower=0), prior="normal(0, 2)")

sho_post <- ode_model(N,
                 odefun_vars = list(w,alpha,gamma), 
                 odefun_body = sho_fun_body,
                 odefun_init = y0,
                 loglik_vars = list(sigma, y1_obs),
                 loglik_body = sho_loglik_body)

data <- read.csv(file = 'simulation.csv') 
data<- data[2:225,]
sho_fit_post <- sho_post$sample( 
  t0 = 0.0, 
  t = data$T,
  data = list(y0 = c(1, 0), D = 2, y1_obs = data$z1),
  refresh = 100,
  solver = midpoint(4),
  max_treedepth = 20,
  parallel_chains = getOption("mc.cores", 4)
)

I tried to use different solver but it’s doesn’t seems to be working. By that, I mean, that even after 30 minutes, none of the chains has done more than 10 iterations. (When considering only the alpha coefficient, the 2000 iterations for 4 chains are done in 10 minutes.)

I tried replacing the absolute value in the ODE with the inverse logit function :


sho_fun_body <- "
  vector[2] dy_dt;
  dy_dt[1] = y[2];
  
  real  s  = - w * sin(y[1]);
        s += - alpha * y[2];
        s += gamma* (-1 + 2*inv_logit(y[2]));
       
      
  dy_dt[2] = s ;

  return(dy_dt);
"

But I am wondering if that’s really the best way to go around that issue. Inv logit is a good approximation of the sign function for large acceleration, but it’s quite bad for acceleration between -5 and 5 which would be the typical range of value taken in our experiment. what do you think ?

untitled
untitled2880×1446 38.4 KB

Or do you have any suggestion ?

edit: Using the inverse logit seems to be allowing Stan to sample data. When using the model on simulated data with \alpha > 0, \gamma = 0,stan is able to recover nicely the parameters. However, when gamma > 0 in the simulated data, it tend to largely underestimate gamma and overestimate alpha. I guess because of the use of the logit function instead of the sign function.

Thanks a lot,
Edouard