error occurred during calling the sampler; sampling not done

[Tommaso Guerrini]

Hi, I just started using STAN and unfortunately keep getting an error on initial values: 

"Rejecting initial value:"                                                                                                                    
[600] "  Error evaluating the log probability at the initial value."                                                                                
[601] "Initialization between (-2, 2) failed after 100 attempts. "

The model I'm trying to implement is in the image. 

data{
int<lower=1> N;  //Number of obs
int<lower=1> K;  //Number of fixed effects
int<lower=1> R;  //number of mixed effects
int<lower=1> R2; //number of single product mixed effects
int<lower=1> I; // number of products
int<lower=1> B; // number of brands
int<lower=1> S; // number of sugars/flavours
matrix[N,R2] Zeta; //matrix for single product mixed effects
matrix[N,K] X; //fixed effects matrix
matrix[N,R] Z; //mixed effects matrix (the same for brand and sugar)
int Y[N,I];  //Observations
vector[K] zerobeta; //prior mean for fixed effects
vector[3] zeroeta;  //prior mean for single product mix effects
vector[R] zeromix;  //prior mean for brand/flavour mix effects
matrix[K,K] Sigma_beta; //now I assume these are known
matrix[3,3] Sigma_eta; //(covariance matrix for eta)
matrix[R,R] Sigma_brand; //(" " " b1)
matrix[R,R] Sigma_sugar; //(" " " b2)
int<lower=1,upper=B> brand[I]; //brand for each prod
int<lower=1,upper=S> sugar[I]; //flavour/sugar for each prod


}


parameters{  
vector[K] beta; //fized effects
vector[3] eta[I]; // random, product specific effects
//in the first version I put R2 here, which is 141, so it was a mistake
vector[R] b1[B]; // brand raneffs
vector[R] b2[S]; // sugar raneffs


}


model {
 


 beta ~ multi_normal(zerobeta , Sigma_beta);


 for(i in 1:I){
 eta[i] ~ multi_normal(zeroeta , Sigma_eta);
 }
 for(j in 1:B){
 b1[j] ~ multi_normal(zeromix , Sigma_brand);
 }
 for(k in 1:S){
 b2[k] ~ multi_normal(zeromix , Sigma_sugar);
 }


 for(i in 1:N){
  for(l in 1:I){
 real lambda;
 ip = c(l,l+47,l+94);
 lambda = X[i,]*beta  + Z[i,]*b1[brand[l]] + Z[i,]*b2[sugar[l]] + Zeta[i,ip]*eta[l];


 Y[i,l] ~ neg_binomial_2(exp(lambda),exp(lambda/2));
}
}


}

Thank you for your help!

[Daniel Lee]

Hi Tommaso,

It'll help if you provide the data in this case. It's hard to debug blind. 

If you're running in RStan, can you try with one chain? I think the parallel mechanism hides some of the output.

And if you're able, try printing some of the parameters while you're running. Hopefully you can see which parts are not behaving well. If I were to debug this, I'd start by commenting out most of the model block (just leaving in the first multi_normal), verifying that works, then add the next bit of the model, until I find where it doesn't run.

Daniel

vector[R2] eta[I]; // random, product specific effects

vector[R] b1[B]; // brand raneffs
vector[R] b2[S]; // sugar raneffs


}


model {
 


 beta ~ multi_normal(zerobeta , Sigma_beta);


 for(i in 1:I){
 eta[i] ~ multi_normal(zeroeta , Sigma_eta);
 }
 for(j in 1:B){
 b1[j] ~ multi_normal(zeromix , Sigma_brand);
 }
 for(k in 1:S){
 b2[k] ~ multi_normal(zeromix , Sigma_sugar);
 }


 for(i in 1:N){
  for(l in 1:I){
 real lambda;
 ip = c(l,l+47,l+94);
 lambda = X[i,]*beta  + Z[i,]*b1[brand[l]] + Z[i,]*b2[sugar[l]] + Zeta[i,ip]*eta[l];


 Y[i,l] ~ neg_binomial_2(exp(lambda),exp(lambda/2));
}
}


}

Thank you for your help!

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Tommaso Guerrini]

Hi Daniel, thank you for your help. First of all, in the code above I sample eta[i] in the parameters vector is a vector[3] not a vector[R2] parameter.

Here the data: 

Y <- allts[,bibex[bgcp]]    //allts is a xts object, bibex[bgcp] are the 47 products I am considering
N <- 729                    //number of observations (2 years)  
I <- 47                     //number of products

X <- cbind(wdtim,os,indfest,indfestlag) //weekday, open/closed store, holiday or not, holiday lagged   (all 0,1)

Z <- cbind(allnets[,bibex[bgcp]],apply(stockouts[,bgcp],2,na.zero)) //matrix for mixed effects: price, availability or not of products (0/1)
R=94 //number of mixed effects
K=9  //number of fixed effects
R2 = 141 //size of Zeta (columns)
Zeta <- cbind(Z,alltslag[,bibex[bgcp]]) //matrix I need for mixed single product effects

brand <- as.numeric(brand)
sugar <- as.numeric(sugar)
Sigma_beta <- matrix(0.01,9,9)  //variance covariance matrix for each effect
Sigma_eta <- matrix(0.01,3,3)
Sigma_brand <- matrix(0.01,94,94)
Sigma_sugar <- matrix(0.01,94,94)

Here the call: 

fit1 <- stan(
  file = "C:/Users/t.guerrini/Desktop/Bennet/estrazione_dati_bennet/Stan code/multigroup.stan",  # Stan program
  data = list(N,K,R,I,B=3,S=3,X,Z,Zeta,R2,Y,
              brand,sugar,zerobeta=rep(0,K),
              zeroeta=rep(0,3),zeromix=rep(0,R),
              Sigma_beta,Sigma_eta,Sigma_brand,Sigma_sugar),    # named list of data
  chains = 1,             # number of Markov chains
  warmup = 1000,          # number of warmup iterations per chain
  iter = 2000,            # total number of iterations per chain
  cores = 2,              # number of cores (using 2 just for the vignette)
  refresh = 1000          # show progress every 'refresh' iterations
)

I'll try with just fixed effects and then increase complexity as you suggest in the meanwhile. Thank you again!

 

[Bob Carpenter]

That error comes when there are values of parameters that satisfy
the declared constraints but where the model block throws an
exception (due to illegal inputs) or returns negative infinity.

For efficiency, you need to:

* vectorize the multi-normals

For soundness, you should declare the data variables as
cov_matrix---that'll test positive definiteness when the
data's read in.

- Bob

> vector[R2] eta[I]; // random, product specific effects

> vector[R] b1[B]; // brand raneffs
> vector[R] b2[S]; // sugar raneffs
>
>
> }
>
>
> model {
>
>
>
> beta ~ multi_normal(zerobeta , Sigma_beta);
>
>
> for(i in 1:I){
> eta[i] ~ multi_normal(zeroeta , Sigma_eta);
> }
> for(j in 1:B){
> b1[j] ~ multi_normal(zeromix , Sigma_brand);
> }
> for(k in 1:S){
> b2[k] ~ multi_normal(zeromix , Sigma_sugar);
> }
>
>
> for(i in 1:N){
> for(l in 1:I){
> real lambda;
> ip = c(l,l+47,l+94);
> lambda = X[i,]*beta + Z[i,]*b1[brand[l]] + Z[i,]*b2[sugar[l]] + Zeta[i,ip]*eta[l];
>
>
> Y[i,l] ~ neg_binomial_2(exp(lambda),exp(lambda/2));
> }
> }
>
>
> }
>
>
>
> Thank you for your help!
>
>
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Tommaso Guerrini]

Thank you Bob, 

by specifying a positive definite matrix I don't get the error I was getting before and sampling starts. 

With vectorizing I guess you mean: 

for(i in 1:I){
 eta[i] ~ multi_normal(zeroeta , Sigma_eta);
 }

becoming

eta ~ multi_normal(zeroeta , Sigma_eta);

(also for the other multi-normals). 

Well, you solved my problem and I'll mark your answer as complete :)

Thank you very much!