Exponential decay running extremely slow

Modeling
1

Hi all,

I am dealing with the following Stan code below, which I call through the following code in R:
ebllist ← list(N = obs, nsubjects = participants, subj_ids = id_participants,
y = y, ts = ts, t0 = 0)

defining the initial values

rAB0 = rep(0, participants)
rmuab = rep(0, participants)
r = cbind(rAB0, rmuab)
str(r)
initebl = function(){
list(logmuab = log(0.05),
logAB0 = 7, logsigma = 0, logsigmaAB0 = -3, logsigma_muab = -2,
rho_rhobit = 0.5, r = r)
}

Expdecay ← stan_model(“C:/Users/IGarciaFogeda/Documents/WP4/DRC/code/Mechanistic models/exponentialdecay.stan”)
fit_expdecay ← sampling(Expdecay,
data = ebllist, init = initebl,
iter = 1000,
chains = 2, warmup = 500,
seed = 0, control = list (stepsize = 0.1))

When I run this I get the message that the Gradient evaluation took 0.032 seconds. However, it is stuck for a while to even reach 10% and it is stuck in the warmup. Where could it possibly be causing this?

functions {
  real[] edecay(real t, real[] y, real[] theta,  real[] x_r, int[] x_i) {
    real dydt[1];
    dydt[1] = theta[1]*y[1];
    
    return dydt;
  }
}
data {
  int<lower=1> N;
  int<lower=1> nsubjects;
  int<lower=1> subj_ids[N];
  real <lower=0> y[N];
  real ts[N];
  real t0;
}
transformed data {
    real x_r[0];
    int x_i[0]; 
}
parameters {
  real logmuab;
  
  real logAB0;
  real logsigma;
  real logsigmaAB0;
  real logsigma_muab;
  real rho_rhobit;
  
  vector[2] r[nsubjects];
}
model {
  real new_mu[N,1];
  real mu[1,1];
  real eval_time[1];
  real theta[1];
  vector[nsubjects] rAB0;
  vector[nsubjects] rmuab;
  
  real inits[1]; 
  
  real muab = exp(logmuab);
  
  real AB0 = exp(logAB0);

  real sigma = exp(logsigma);
  real sigmaAB0 = exp(logsigmaAB0);
  real sigma_muab = exp(logsigma_muab);
  
  vector[2] mean_vec;
  matrix[2,2] Sigma; 
  real rho = (exp(rho_rhobit) - 1)/(exp(rho_rhobit) + 1);
  
  //prior distributions
  logmuab ~ normal(0, 10);
  
  logAB0 ~ normal(7, 2);

  logsigmaAB0 ~ normal(0, 1);
  logsigma_muab ~ normal(0, 1);
  rho_rhobit ~ uniform(-10,10);
  
  mean_vec[1] = -pow(sigmaAB0, 2.0)/2;
  mean_vec[2] = -pow(sigma_muab,2.0)/2;
  
  Sigma[1,1] = pow(sigmaAB0, 2.0);
  Sigma[2,2] = pow(sigma_muab, 2.0);
  Sigma[1,2] = rho*sigmaAB0*sigma_muab;
  Sigma[2,1] = rho*sigmaAB0*sigma_muab;
  
  for (subj_id in 1:nsubjects) {
        r[subj_id] ~ multi_normal(mean_vec, Sigma);
        rAB0[subj_id] = r[subj_id,1];
        rmuab[subj_id] = r[subj_id,2];
  }
  
  for (obs_id in 1:N){
    inits[1] = AB0*exp(rAB0[subj_ids[obs_id]]);
    
    theta[1] = muab*exp(rmuab[subj_ids[obs_id]]); 
     
    eval_time[1] = ts[obs_id];
    if (eval_time[1] == 0){
      new_mu[obs_id,1] = log(inits[1]) - pow(sigma, 2.0)/2;
      }
    else{
      mu = integrate_ode_rk45(edecay, inits, t0, eval_time, theta, x_r, x_i);
      new_mu[obs_id,1] = log(mu[1,1]) - pow(sigma, 2.0)/2;
    }  
    //likelihood function for AB levels
    y[obs_id] ~ lognormal(new_mu[obs_id,1], sigma); 
  }
}
generated quantities {
  real z_pred[N];
  real sigma2 = exp(logsigma);
  for (t in 1:N){
   z_pred[t] = lognormal_rng(log(y[t] + 1e-5), sigma2); 
  }
}
2

I’m pretty sure theta is positive so this would be exponential growth, not decay.

Is this uniform prior meant to approximate a Beta(0,0) prior for rho? Could you try a uniform prior for rho instead?

Anyway, this ODE is so simple that it can be solved analytically. Try and replace that line with

mu[1,1] = inits[1]*exp(-theta[1] * (eval_time[1] - t0));

(where I’ve added the minus sign for decay)

1 Like
3

Hi,

Thanks! That’s right, theta should have a minus. However, I do not entirely know what you mean by the prior of rho? I am already using a uniform distribution for rho?

4

No, you have a uniform prior for rho_rhobit. The prior for rho is highly concentrated near the extremes.

rho_rhobit <- runif(1000, -10, 10)
rho <- (exp(rho_rhobit)-1) / (exp(rho_rhobit) + 1)
hist(rho_rhobit)
hist(rho)

rho_rhobit
rho_rhobit691×726 20.4 KB

rho
rho691×726 19.5 KB

3 Likes