I am trying to run the Skew-QBS ACD model for right skewed duration data in stan. when running these code i am getting an error

Modeling
1 
qbsacd2="functions{

real eta(real Q, real alpha){

real eta;

eta=alpha*Q+((alpha*Q)^2+4)^(1/2);

return eta;

}

vector a1(vector y,real Q, real alpha, vector E){

vector[num_elements(y)]a1;

for(i in 1:num_elements(y)){

a1[i]=1/alpha*((y[i]*(eta(Q,alpha))^2/(4*E[i]))^(1/2)-(4*E[i]/(y[i]*(eta(Q,alpha))^2))^(1/2));

}

return a1;

}

vector log_a2(vector y,real Q, real alpha, vector E){

vector[num_elements(y)]log_a2;

for(i in 1:num_elements(y)){

log_a2[i]=log(1/(2*alpha*y[i])*((y[i]*(eta(Q,alpha))^2/(4*E[i]))^(1/2)-(4*E[i]/(y[i]*(eta(Q,alpha))^2))^(1/2)));

}

return log_a2;

}

real sqcd_lpdf(vector y,real Q,real lambda, real alpha, vector E){

vector[num_elements(y)]lliksqcd;

real prob;

lliksqcd=log(2)+std_normal_lpdf(a1(y,Q,alpha,E))+normal_lcdf(lambda*a1(y,Q,alpha,E)|0,1)+log_a2(y,Q,alpha,E);

prob=sum(lliksqcd);

return(prob);

}

}

data{

int N;

vector<lower=0>[N]y;

real Q;

}

parameters{

real<lower=0>alpha;

real<lower=0>lambda;

real<lower=0,upper=1>omega;

real<lower=0,upper=1>rho;

real<lower=0,upper=1-rho>sigma;

}

transformed parameters{

vector<lower=0>[N]E;

E[1]=1;

for(i in 2:N){

E[i]=exp(omega+rho*log(E[i-1])+(sigma*y[i-1])/E[i-1]);

}

}

model{

alpha~cauchy(0,25)T[0,];

lambda~inv_gamma(0.001,0.001);

omega~normal(0,5)T[0,1];

rho~normal(0,5)T[0,1];

sigma~normal(0,5)T[0,1-rho];

y~sqcd(Q,lambda,alpha,E);

}

"

require(rstan)

N=length(mydata)

N

dacd=list(y=mydata,N=N,Q=1.141535)

Qbsacd=stan(model_code=qbsacd2,data=dacd,iter=5000,chains=2)


When i run this code stan rejecting the initial values, and the error is given below

Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
[1] "error occurred during calling the sampler; sampling not done"

Does anyone know what is going wrong with my code, please help

[edit: formatted code]

2

I escaped the code in the original post, but am posting a reformatted version here.

functions {
  real eta(real Q,  real alpha)  {
    real eta;
    eta = alpha * Q + ((alpha * Q) ^2 + 4) ^(1 / 2);
    return eta;
  }
  vector a1(vector y, real Q,  real alpha,  vector E)  {
    vector[num_elements(y)] a1;
    for(i in 1:num_elements(y))  {
      a1[i] = 1 / alpha * ((y[i]  * (eta(Q, alpha)) ^2 / (4 * E[i])) ^(1 / 2)  - (4 * E[i]  / (y[i]  * (eta(Q, alpha)) ^2)) ^(1 / 2));
    }
    return a1;
  }
  vector log_a2(vector y, real Q,  real alpha,  vector E)  {
    vector[num_elements(y)] log_a2;
    for(i in 1:num_elements(y))  {
      log_a2[i] = log(1 / (2 * alpha * y[i])  * ((y[i]  * (eta(Q, alpha)) ^2 / (4 * E[i])) ^(1 / 2)  - (4 * E[i]  / (y[i]  * (eta(Q, alpha)) ^2)) ^(1 / 2)) );
    }
    return log_a2;
  }
  real sqcd_lpdf(vector y, real Q, real lambda,  real alpha,  vector E)  {
    vector[num_elements(y)] lliksqcd;
    real prob;
    lliksqcd = log(2) + std_normal_lpdf(a1(y, Q, alpha, E)) + normal_lcdf(lambda * a1(y, Q, alpha, E) | 0, 1) + log_a2(y, Q, alpha, E);
    prob = sum(lliksqcd);
    return(prob);
  }
}
data {
  int N;
  vector<lower = 0>[N] y;
  real Q;
}
parameters {
  real<lower = 0> alpha;
  real<lower = 0> lambda;
  real<lower = 0, upper = 1> omega;
  real<lower = 0, upper = 1> rho;
  real<lower = 0, upper = 1 - rho> sigma;
}
transformed parameters {
  vector<lower = 0> [N] E;
  E[1] = 1;
  for(i in 2:N)  {
    E[i] = exp(omega + rho * log(E[i - 1]) + (sigma * y[i - 1])  / E[i - 1]);
  }
}
model {
  alpha ~ cauchy(0, 25) T[0, ];
  lambda ~ inv_gamma(0.001, 0.001);
  omega ~ normal(0, 5) T[0, 1];
  rho ~ normal(0, 5) T[0, 1];
  sigma ~ normal(0, 5) T[0, 1 - rho];
  y ~ sqcd(Q, lambda, alpha,E);
}

The error message says what’s going wrong:

Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: 
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

To paraphrase, it’s trying to initialize randomly uniform(-2, 2) on the unconstrained scale and it’s running into values that are out of support. That is the log density evaluates to zero (or negative infinity on the log scale).

It might help to start with a simpler model and then add to it to build up to something this complex. You can throw prints in of target(), i.e., just print("checkpoint 3, target = ", target()) to print the log pdf at every stage of the calculation.

I wasn’t clear on why the density defined by sqcd have both a standard normal for a1 and a scaled standard normal for a1?

Some coding tips:

  • Use single spacing in a standalone file—then your error messages have meaningful line numbers and you can read more of the code at once
  • You don’t need truncations on densities with constant parameters—the truncation is constant, so not needed for MCMC.
  • The functions can be vectorized to one-liners.
  • Unless you need E in the output, define it as a local variable in the model block.
1 Like
3

Thank you, sir, sorry for the late replay, I used reparameterisation for this problem.

data{
int N;
vector<lower=0>[N]y;
real Q;
}
parameters{
real <lower=0>u;
real<lower=0>lambda;
real<lower=0>v;
real<lower=0>w;
real<lower=0> sigmav;
real<lower=0> sigmaw;
}
transformed parameters{
real alpha;
vector<lower=0>[N]E;
vector<lower=1>[N]theta;
theta[1]=1;
for(i in 2:N){
theta[i]=theta[i-1]+v;
}
for(i in 1:N){
E[i]=exp(theta[i]+w);
}
alpha=1/sqrt(u);
}
model{
v~inv_gamma(5,0.16);
w~inv_gamma(5,0.16);
sigmav~gamma(2,0.5);
sigmaw~gamma(2,0.5);
u~cauchy(0,25)T[0,];
lambda~inv_gamma(5, 0.16);
y~sqcd(Q,lambda,alpha,E);
}
But again, it shows the same error.
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can’t start sampling from this initial value.
Chain 1:
Chain 1: Initialization between (-2, 2) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”