Error evaluating the log probability at the initial value

Modeling
1

Dear all,

I am trying to fit to my data differential equations models that describe the inter-relationship among susceptible cells, infected but not productive cells, infected and virus productive cells and the number of free viruses with Stan to estimate my model parameters.
My model has 5 parameters.
This is working fine with 4 parameters. Stan produces posterior distributions for all 4 parameters and initial value estimates. But I am constantly failing with the fifth (ie. parameter 1, rate constant characterizing infection). I got this error message:

SAMPLING FOR MODEL ‘VIR5param_fit’ NOW (CHAIN 1).
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: Exception: []: accessing element out of range. index 5 out of range; expecting index to be between 1 and 4; index position = 1params (in ‘model3e34b756a99_VIR5param_fit’ at line 15)
(in ‘model3e34b756a99_VIR5param_fit’ at line 57)

[1] “Error in sampler$call_sampler(args_list[[i]]) : "
[2] " Exception: Exception: []: accessing element out of range. index 5 out of range; expecting index to be between 1 and 4; index position = 1params (in ‘model3e34b756a99_VIR5param_fit’ at line 15)”
[3] " (in ‘model3e34b756a99_VIR5param_fit’ at line 57)"
[1] “error occurred during calling the sampler; sampling not done”

Modifying the “max_num_steps” is not solving the issue. I might have a problem related to this topic: Error at initial value. My parameter value is very small eg: 1e-7. I tried to replace:

  dydt[1] = - params[1] * y[1] * y[4];
  dydt[2] = params[1] * y[1] * y[4] - params[2] * y[2];

By

  dydt[1] = - pow(10,params[1]) * y[1] * y[4];
  dydt[2] = pow(10,params[1]) * y[1] * y[4] - params[2] * y[2];

And sample in log- transformed values in a normal distribution, but it also failling with the code: “Chain 1: Unrecoverable error evaluating the log probability at the initial value.”

Here is my Stan code:

  functions {
  
  real[] VIR(real ts,
  real[] y,
  real[] params,
  real[] x_r,
  int[] x_i) {
  
  real dydt[4];
  
  dydt[1] = - params[1] * y[1] * y[4];
  dydt[2] = params[1] * y[1] * y[4] - params[2] * y[2];
  dydt[3] = params[2] * y[2] - params[3] * y[3];
  dydt[4] = params[4] * y[3] - params[5] * y[4]  ;
  
  return dydt;
  }

  }
  
  data {
  int<lower = 1> n_obs; // Number of days sampled
  int<lower = 1> n_params; // Number of model parameters
  int<lower = 1> n_difeq; // Number of differential equations in the system
  int<lower = 1> n_obs_Pred; // This is to generate "predicted" data
  
  real y[n_obs]; // The data
  real t0; // Initial time point (zero)
  real ts[n_obs]; // Time points that were sampled

  real Pred_ts[n_obs_Pred]; // Predicted data
  
  
  }
  
  transformed data {
  real x_r[0];
  int x_i[0];
  }
  
  parameters {
  real<lower = 0> params[n_params]; // Model parameters
  real<lower = 100, upper = 100000000> S0; // Initial fraction of susceptible cells constrained to be between 100 and 1e8
  real<lower = 10, upper = 1000> V0; // Initial fraction of infectious viruses constrained to be between 10 and 1e3
  }
  
  transformed parameters{
  real y_hat[n_obs, n_difeq]; // Output from the ODE solver
  real y0[n_difeq]; // Initial conditions
  
  y0[1] = S0;
  y0[2] = 0;
  y0[3] = 0;
  y0[4] = V0;

  y_hat = integrate_ode_rk45(VIR, y0, t0, ts, params, x_r, x_i);
  
  }

  model {
  params[1] ~ normal(0.00001, 0.01); //constrained to be positive
  params[2] ~ normal(2, 1); //constrained to be positive
  params[3] ~ normal(6, 3); //constrained to be positive
  params[4] ~ normal(25000, 1000); //constrained to be positive
  params[5] ~ normal(10, 3); //constrained to be positive
  
  S0 ~ normal(10000000, 3000000); //constrained to be positive
  V0 ~ normal(100, 50); //constrained to be positive
  
  for (t in 1:n_obs) {y[t] ~ normal(y_hat[t,4], y_hat[t,4]/10);} // y_hat[,4] is the viremia fom the ODE solver
  
  }
  
  generated quantities {
  // Generate predicted data over the whole time series:
  
  real Pred_I[n_obs_Pred, n_difeq];
  
  Pred_I = integrate_ode_rk45(VIR, y0, t0, Pred_ts, params, x_r, x_i);
  
  }

I run it with rstan as follow:

stan_d = list(n_obs = 11,
n_params = 5,
n_difeq = 4,
n_obs_Pred = 16,
y = c(2000000,6000,941,240,50,10,0,0,0,0,0),
t0 = 0,
ts = c(5,7,8,9,10,11,12,13,14,15,16),
Pred_ts = c(1:max(sample_time)))

Test the model:

test = stan(“VIR5param_fit.stan”,
data = stan_d,
pars = params_monitor,
chains = 1, iter = 10)

This is new for me and I feel that the solution is just around the corner, but I cannot go through this! Clear examples would be appreciated :)

Thank you for your help!
Albin

1 Like
2

Looking at the error message alone, it is saying you have an indexing error on line 57, related to the function you define. It’s looking for the fifth element of some vector which has only four elements.

As far as I can see you mix supplying the length of vectors as data with a fixed index (params[5] on one hand and declaring params[n_params] on the other). Could something related to that be tripping you up?

1 Like
3

You have indeed spotted an issue! Thank you!!

I fix it, but it is still failing :( with this error message:

SAMPLING FOR MODEL ‘VIR5param_fit’ NOW (CHAIN 1).
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: Max number of iterations exceeded (1000000). (in ‘model29684ec4ea_VIR5param_fit’ at line 57)

[1] “Error in sampler$call_sampler(args_list[[i]]) : "
[2] " Exception: Max number of iterations exceeded (1000000). (in ‘model29684ec4ea_VIR5param_fit’ at line 57)”
[1] “error occurred during calling the sampler; sampling not done”

I don’t think the model is ill-defined because it is running fine with this set of parameter values.

In R:

Load required packages

##############
library(tidyverse)
library(deSolve)
library(dplyr)
library(ggplot2)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
##############

#Data:
df <- data.frame(Days=c(5,7,8,9,10,11,12,13,14,15,16), V=c(2e+06,6e+03,941,240,50,10,0,0,0,0,0))

within.host.model <- function (t, y, params) {

Ta <- y[1] #Starting value Target cells per mL
I1 <- y[2] #Starting value of infected but not productive cells
I2 <- y[3] #Starting value productively infected cells
V <- y[4] #Starting value free virus

beta <- params[1]
k <- params[2]
kappa <- params[3]
p <- params[4]
c <- params[5]

differential equations

dTa <- - betaTaV #Susceptible cells
dI1 <- betaTaV - kI1 # infected but not productive cells
dI2 <- kI1 - kappaI2 #productively infected cells
dV <- pI2 - c*V #free virus

list(c(dTa,dI1,dI2, dV))
}

####Parameter values
params <- c(
beta <-1/1e7,
k <- 2.1,
kappa <- 6,
p <- 25e3,
c <- 10)

times <- seq(from=0,to=17,by=1)
xstart <- c(dTa =1e7,dI1=0,dI2=0, dV=1e2 )

ode(
func=within.host.model,
y=xstart,
times=times,
parms=params,
method=“ode45”
) %>%
as.data.frame() -> out

out2 <- out %>% gather(variable,value,-time)
out2$variable <- str_replace(out2$variable, “dTa”, “Susceptible cells” )
out2$variable <- str_replace(out2$variable, “dI1”, “Infected but not productive cells” )
out2$variable <- str_replace(out2$variable, “dI2”, “Productively infected cells” )
out2$variable <- str_replace(out2$variable, “dV”, “Free virus (titre copy/mL)” )

pl2 <- ggplot(subset(out2, variable==“Free virus (titre copy/mL)”), aes(x=time,y=log10(value+1)))
pl2 <- pl2 + geom_vline(aes(xintercept = 5.8), color=“orange”, linetype=“longdash”)
pl2 <- pl2 + geom_line(size=1)
pl2 <- pl2 + geom_point(data=df, aes(x=Days, y=log10(V+1)), size=3)
pl2 <- pl2 + theme_minimal(base_size =20)
pl2 <- pl2 + guides(colour=FALSE) + xlim(0,17) + ylim(0,10)
pl2

4

Hmm. There is still something funky going on at line 57, where one of the datasets makes trouble and the other doesn’t. If I was debugging this, that’s where I’d start to try to figure out what is going wrong - what’s different between these two datasets at that line and not before.

Apart from that, your model is way over my head and mathematical skills, so I can’t really be of any help, but there might be others!

5

I messed with the stan code for a while and got a version that returns Rhats of 1.00 and neff of at least several hundred with 4 chains and iter = 1000. I made changes along two lines. First, looking at your R code, I noticed that the equivalent of params[1] was set to 1e-7 but the Stan prior for that was defined with lower=0, centered at 1e-5 and with a sigma of 0.01, so it would be orders of magnitude larger than 1e-7 most of the time. I changed params[1] to be centered at 0 with a sigma of 1e-7.
Second, I changed the parameterization to put priors of normal(0,1) on parameters I called paramsRaw, Sraw and Vraw and then transformed those to the params, S0, V0 parameters by multiplying by the desired sigma and added the desired mean. Section 22.7 of the 2.21 users guide talks about this.
I do not claim that this meets your needs because I do not understand what you are doing. I am particularly queasy about the bounds on S0 and V0. But it does run.
Edit: I forgot to mention there were errors like this
“SAMPLING FOR MODEL ‘8cd9afdf3c9c290dc3a532bf2a603908’ NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: normal_lpdf: Scale parameter is -4.85105e-09, but must be > 0”

MODEL <- "
functions {
  
  real[] VIR(real ts,
             real[] y,
             real[] params,
             real[] x_r,
             int[] x_i) {
    
    real dydt[4];
    
    dydt[1] = - params[1] * y[1] * y[4];
    dydt[2] = params[1] * y[1] * y[4] - params[2] * y[2];
    dydt[3] = params[2] * y[2] - params[3] * y[3];
    dydt[4] = params[4] * y[3] - params[5] * y[4]  ;
    
    return dydt;
  }
  
}

data {
  int<lower = 1> n_obs; // Number of days sampled
  int<lower = 1> n_params; // Number of model parameters
  int<lower = 1> n_difeq; // Number of differential equations in the system
  int<lower = 1> n_obs_Pred; // This is to generate predicted data
  
  real y[n_obs]; // The data
  real t0; // Initial time point (zero)
  real ts[n_obs]; // Time points that were sampled
  
  real Pred_ts[n_obs_Pred]; // Predicted data
  
  
}

transformed data {
  real x_r[0];
  int x_i[0];
}

parameters {
  real<lower = 0> paramsRaw[n_params]; // Model parameters
  real<lower = 0> Sraw; // Initial fraction of susceptible cells constrained to be between 100 and 1e8
  real<lower = 0> Vraw; // Initial fraction of infectious viruses constrained to be between 10 and 1e3
}

transformed parameters{
  real y_hat[n_obs, n_difeq]; // Output from the ODE solver
  real y0[n_difeq]; // Initial conditions
  real params[n_params]; // Model parameters
  real<lower = 100, upper = 100000000> S0; // Initial fraction of susceptible cells constrained to be between 100 and 1e8
  real<lower = 10, upper = 1000> V0; // Initial fraction of infectious viruses constrained to be between 10 and 1e3
  
  S0 = Sraw * 3000000 + 10000000; //constrained to be positive
  V0 = Vraw * 50 + 100; //constrained to be positive
  y0[1] = S0;
  y0[2] = 0;
  y0[3] = 0;
  y0[4] = V0;
  params[1] = paramsRaw[1] * 0.0000001 ; //constrained to be positive * 0.01+ 0.00001
  params[2] = paramsRaw[2] + 2; //constrained to be positive
  params[3] = paramsRaw[3] * 3 + 6; //constrained to be positive
  params[4] = paramsRaw[4] * 1000 + 25000; //constrained to be positive
  params[5] = paramsRaw[5] * 3 + 10; //constrained to be positive
  y_hat = integrate_ode_rk45(VIR, y0, t0, ts, params, x_r, x_i);
}

model {
  paramsRaw ~ normal(0, 1);
  Sraw ~ normal(0, 1);
  Vraw ~ normal(0, 1);
  for (t in 1:n_obs) {y[t] ~ normal(y_hat[t,4], y_hat[t,4]/10);} // y_hat[,4] is the viremia fom the ODE solver
  
}

//generated quantities {
  // Generate predicted data over the whole time series:
    
  //  real Pred_I[n_obs_Pred, n_difeq];
  
  //Pred_I = integrate_ode_rk45(VIR, y0, t0, Pred_ts, params, x_r, x_i);
  
//}
"

stan_d = list(n_obs = 11,
              n_params = 5,
              n_difeq = 4,
              n_obs_Pred = 16,
              y = c(2000000,6000,941,240,50,10,0,0,0,0,0),
              t0 = 0,
              ts = c(5,7,8,9,10,11,12,13,14,15,16),
              Pred_ts = c(1:16)) #max(sample_time)

  
  test = stan(model_code = MODEL,
               data = stan_d,
               #pars = params_monitor,
               chains = 4, iter = 1000)

> summary(test)$summary[52:63, c("mean", "se_mean", "n_eff", "Rhat")]
                   mean      se_mean     n_eff     Rhat
y0[1]      1.265964e+07 4.944756e+04 1478.3639 0.999103
y0[2]      0.000000e+00          NaN       NaN      NaN
y0[3]      0.000000e+00          NaN       NaN      NaN
y0[4]      1.389633e+02 8.193793e-01 1287.1200 1.002927
params[1]  6.455663e-08 1.274318e-09  725.1741 1.000889
params[2]  2.047288e+00 5.081819e-04  939.8764 1.003017
params[3]  7.652864e+00 4.753549e-02  900.3554 1.002701
params[4]  2.583748e+04 1.453733e+01 1740.8339 1.000074
params[5]  1.197085e+01 4.760474e-02 1105.1130 1.001607
S0         1.265964e+07 4.944756e+04 1478.3639 0.999103
V0         1.389633e+02 8.193793e-01 1287.1200 1.002927
lp__      -3.366496e+02 1.044012e-01  472.4575 1.006635
2 Likes
6

You are my hero :)

This is a nice troubleshooting, I couldn’t have done it myself for sure.
Puting priors of normal(0,1) and transform them in a second step was smart.
You are right, these bounds on S0 and V0 are arbitrary, I am going to “polish” this code now that I have a runing version.

Thank you for your help and for this Stan lesson. It is so nice to have support like this!

2 Likes