Problem in GQ block: binomial_logit_log: Probability parameter is nan, but must be finite!

[Colin Olito]

Hello all,

Apologies if this is just a stupid error on my part, but I was wondering if I could get some more eyes to help me figure out why my generated quantities block is breaking the sampler.  I am trying to fit a random intercept/slope logistic regression with covariance matrix for grouped bernoulli trial data. I *think* the model itself is working ok, but the minute I added a generated quantities block to calculate the log_lik for subsequent use with 'loo', the sampler stopped working (the exact same model file runs fine and returns reasonable looking results if I comment out the GQ block). With the GQ block included, the model still compiles fine, but immediately throws an error and returns no samples. At the same time, nearly the same GQ block (just without the [grp[]] indexing) works fine for a reduced model with no covariance matrix. Can anybody help me identify what I'm doing here that would give the sampler such a hard time and how to fix? I've copied the model and error message below. 

Thanks in advance for any help or advice.
- C

THE MODEL:

data {
   int<lower=1> N;               // Number of observations
   int<lower=0> P;               // Number of predictors
   int<lower=0> J;               // Number of groups (Runs)
   int<lower=0> K;               // Number of group-level predictors
   int<lower=1,upper=J> grp[N];  // Group index for each observation
   int<lower=0> nT[N];           // Number of trials
   int<lower=0> nS[N];           // Number of successes
   matrix [N,P] X;               // Predictor matrix
   row_vector[K] Z[N];           // Group predictor matrix (row_vector for efficiency)
}

parameters {
   cholesky_factor_corr[K] L_run;   // Cholesky factor of Run r.e. correlations
   vector<lower=0>[K] tau_run;      // standard deviations of unconditional Run r.e. dist
   vector[K] u[J];                  // spherical Run random effects
   vector [P] beta;                 // Vector of fixed effect estimates
   real<lower=0> sigma_y;           // error scale
}

transformed parameters {
   matrix[K,K] corrs;
   corrs <- tcrossprod(L_run);      // for monitoring Run correlations
}

model {

   // Local variables
   vector[N] mu;
   matrix[K,K] Lambda_run;
   vector[K] gamma[J];

   // Priors for covariance matrix
   tau_run  ~  cauchy(0,5);
   L_run    ~  lkj_corr_cholesky(2);
   Lambda_run  <-  diag_pre_multiply(tau_run,L_run);
   for (i in 1:J) {
      u[i] ~ normal(0,1);
      gamma[i] <- Lambda_run * u[i];
   }

   // Priors
   beta     ~  normal(0,3);
   sigma_y  ~  cauchy(0,5);

   // Likelihood
  for (i in 1:N) {
      mu[i]  <-  X[i] * beta + Z[i] * gamma[grp[i]];
  }
  nS  ~  binomial_logit(nT, mu);
}

generated quantities {
    vector[N] log_lik;
    vector[K] gamma[J];
 
    for (i in 1:N) {
      log_lik[i]  <-  binomial_logit_log(nS[i], nT[i], X[i]*beta + Z[i]*gamma[grp[i]]);
    }
}

THE ERROR MESSAGE:

> ## Call to STAN
> mat4 <- stan(data = data.list,
                          file = './Stan/mat-logistic-1Z-cov.stan',
                           chains = 1,
                           iter = numSavedSteps,
                           thin = thinSteps,
                           save_dso = TRUE
                           )
COMPILING THE C++ CODE FOR MODEL 'mat-logistic-1Z-cov' NOW.
starting worker pid=28067 on localhost:11870 at 15:11:51.074
starting worker pid=28076 on localhost:11870 at 15:11:51.213
starting worker pid=28085 on localhost:11870 at 15:11:51.354
starting worker pid=28094 on localhost:11870 at 15:11:51.493

SAMPLING FOR MODEL 'mat-logistic-1Z-cov' NOW (CHAIN 1).
Error : Exception thrown at line 64:
stan::math::binomial_logit_log: Probability parameter is nan, but must be finite!
error occurred during calling the sampler; sampling not done
here are whatever error messages were returned
[[1]]
Stan model 'mat-logistic-1Z-cov' does not contain samples.

Warning message:In .local(object, ...) : some chains had errors; consider specifying chains = 1 to debug
>





I also tried providing intial values for the sampler. Again, this worked fine with the GQ commented out, but gave the same error if I included the GQ.

inits  <-  list(
  list(L_run    =  matrix(runif(16^2,-0.5,0.5), nrow=16,ncol=16)),
  list(tau_run  =  runif(16,0.1,2)),
  list(u        =  matrix(runif((8*16),-1,1),nrow=8,ncol=16)),
  list(beta     =  runif(2,-1,1)),
  list(sigma_y  =  runif(1,0.1,2))
  )

[Jonah Gabry]

Looks like you're declaring 

vector[K] gamma[J]; 

in generated quantities and then using gamma without defining it. The gamma from the model block isn't accessible in generated quantities so you'll either need to define gamma again in generated quantities or move it from the model block to transformed parameters, in which case it will be accessible in generated quantities without having to redefine it. 

Hope that helps,

Jonah

[Colin Olito]

Hi Jonah,

Ah!  Ok.  Re-declaring gamma and Lambda_run in the generated quantities block like this seems to have worked:

generated quantities {

  // Local variables
  vector[N] log_lik;

  matrix[K,K] Lambda_run;
  vector[K] gamma[J];

 
  Lambda_run  <-  diag_pre_multiply(tau_run,L_run);

  for (i in 1:J) {

    gamma[i] <- Lambda_run * u[i];
  }

  for (i in 1:N) {
    log_lik[i]  <-  binomial_logit_log(nS[i], nT[i], X[i]*beta + Z[i]*gamma[grp[i]]);
   }
}

However, I kept running into declaration problems when I tried to move gamma to transformed parameters. Gave up since the above generated quantities block seemed to be working.

So, just to make sure I'm clear on how the pieces fit together:
1) Sampling statements can only occur in the model block, but are accessible in the generated quantities.
2) Local variables in model block are not accessible in generated quantities.
3) I can re-define variables like Lambda_run and gamma in the generated quantities block if needs be (although this results in redundant code).
4) Or I can declare variables I'd like to have in both model and generated quantities blocks as transformed parameters... but only if it doesn't require sampling statements or  result in declaration conflicts.

am I on the right track? Thanks for the help, I really appreciate it.

[Jonah Sol Gabry]

Hi Colin, some comments below:

On Monday, March 28, 2016, Colin Olito <colin...@gmail.com> wrote:

So, just to make sure I'm clear on how the pieces fit together:
1) Sampling statements can only occur in the model block, but are accessible in the generated quantities. 

If by sampling statements you mean statements involving "~", then those go only in the model block. All they do is increment the total log probability accumulator.

 

2) Local variables in model block are not accessible in generated quantities. 

Right. Also, everything in the model block is done once per leapfrog step and everything in GQ is once per iteration (much more efficient).

 

3) I can re-define variables like Lambda_run and gamma in the generated quantities block if needs be (although this results in redundant code). 

To reduce cutting and pasting I recommend creating some functions in the functions block and using the functions where needed throughout the rest of the Stan program.

4) Or I can declare variables I'd like to have in both model and generated quantities blocks as transformed parameters... but only if it doesn't require sampling statements or  result in declaration conflicts.

Right, you can't use "~" in transformed parameters. The transformed parameters block is for defining (deterministic) functions of the unknown parameters. The primitive parameters then get priors in the model block.

And regarding declaration conflict, if you use transformed parameters then you shouldn't need to declare anything with the same name(s) anywhere else so you should be ok. 

am I on the right track? Thanks for the help, I really appreciate it. 

No problem, happy to help. Can you share the code for the model that wasn't working when using transformed parameters? It will be easier to debug that way.

Best,

Jonah

[Marco Inacio]

On 28/03/2016 18:48, Colin Olito wrote:

Hi Jonah,

Ah!  Ok.  Re-declaring gamma and Lambda_run in the generated quantities block like this seems to have worked:

generated quantities {

  // Local variables
  vector[N] log_lik;
  matrix[K,K] Lambda_run;
  vector[K] gamma[J];
 
  Lambda_run  <-  diag_pre_multiply(tau_run,L_run);
  for (i in 1:J) {
    gamma[i] <- Lambda_run * u[i];
  }
  for (i in 1:N) {
    log_lik[i]  <-  binomial_logit_log(nS[i], nT[i], X[i]*beta + Z[i]*gamma[grp[i]]);
   }
}

However, I kept running into declaration problems when I tried to move gamma to transformed parameters. Gave up since the above generated quantities block seemed to be working.

You are probably trying to assign values for gamma inside model block while it's declared as transformed parameter.

Try this:

transformed parameters {
   matrix[K,K] corrs;

   vector[N] mu;


   matrix[K,K] Lambda_run;
   vector[K] gamma[J];



   corrs <- tcrossprod(L_run);      // for monitoring Run correlations



   Lambda_run  <-  diag_pre_multiply(tau_run,L_run);
   for (i in 1:J)
     gamma[i] <- Lambda_run * u[i];
  

   for (i in 1:N)


      mu[i]  <-  X[i] * beta + Z[i] * gamma[grp[i]];

}
model {


   // Priors for covariance matrix
   tau_run  ~  cauchy(0,5);
   L_run    ~  lkj_corr_cholesky(2);


   for (i in 1:J)


      u[i] ~ normal(0,1);




   // Priors
   beta     ~  normal(0,3);
   sigma_y  ~  cauchy(0,5);

   // Likelihood


   nS  ~  binomial_logit(nT, mu);
}
generated quantities {
    vector[N] log_lik;


    for (i in 1:N)


      log_lik[i]  <-  binomial_logit_log(nS[i], nT[i], mu[i]);
}

[Colin Olito]

Hello Jonah, Marco,

Thanks for the advice Jonah - think I'm starting to get the hang of it. And yes, Marco, that's exactly what I was doing (trying to assign values to gamma in the model block instead of the transformed parameters block). I realized it as soon as I looked at your code. I've now incorporated the changes you suggested, and everything runs beautifully.  Thanks again.  I've posted the final model code, in case this is helpful for anybody else, or in case anybody sees any other problems.

/*
 * Basic, flexible logistic ME Regression model
 * using design matrices to define model structure.
 * Covariance matrices estimated.
 */

data {
   int<lower=1> N;               // Number of observations
   int<lower=0> P;               // Number of predictors
   int<lower=0> J;               // Number of groups (Runs)
   int<lower=0> K;               // Number of group-level predictors
   int<lower=1,upper=J> grp[N];  // Group index for each observation
   int<lower=0> nT[N];           // Number of trials
   int<lower=0> nS[N];           // Number of successes
   matrix [N,P] X;               // Predictor matrix
   row_vector[K] Z[N];           // Group predictor matrix (row_vector for efficiency)
}

parameters {
   cholesky_factor_corr[K] L_run;   // Cholesky factor of Run r.e. correlations
   vector<lower=0>[K] tau_run;      // standard deviations of unconditional Run r.e. dist
   vector[K] u[J];                  // spherical Run random effects
   vector [P] beta;                 // Vector of fixed effect estimates
   real<lower=0> sigma_y;           // error scale
}

transformed parameters {

  matrix[K,K] corrs;              
  vector[N] mu;                  
  matrix[K,K] Lambda_run;
  vector[K] gamma[J];

  corrs <- tcrossprod(L_run);       // For monitoring correlations

  Lambda_run  <-  diag_pre_multiply(tau_run,L_run);
  for (i in 1:J)

    gamma[i] <- Lambda_run * u[i];  // Random-effect coefficients
 
  for (i in 1:N)                    // Loop to vectorize likelihood in model block

    mu[i]  <-  X[i] * beta + Z[i] * gamma[grp[i]];
}

model {

  tau_run  ~  cauchy(0,5);          // Priors for covariance matrix

  L_run    ~  lkj_corr_cholesky(2);
  for (i in 1:J)
    u[i] ~ normal(0,1);

  beta     ~  normal(0,3);          // Priors on fixed effects
  sigma_y  ~  cauchy(0,5);          // Prior on error scale

  nS  ~  binomial_logit(nT, mu);    // Likelihood
}

generated quantities {
  vector[N] log_lik;  // Log-likelihood