Constrain a probability parameter in the binomial likelihood

Modeling
1

Hi! I ran a bayesian hierarchical model to my dataset and manually edit log likelihood function for the probability parameter. I ran into error because some of the probability (p[I]) is 1, thus, became 0 after log.
Is there a way to constrain the p[I] to never exceed 0.99 in the model block? or shall I transformed this p[I] in the transformed parameters block first?


data{
    int<lower=0> M; // data length
    int<lower=0> N; 
    int ca[M]; 
    int isol[M]; 
    int s0[M]; 
    int i0[M];  
    real dt[M];
    matrix[M,N] timeck;
    int pens;
    real logdelta;
}
//
parameters{
    vector[pens] z; //intercept of cluster
    real loga_bar; // log transmission rate
}
model{
    vector[M] p;  // probability variable in binomial likelihood
    vector[M] et;  
    loga_bar ~ normal(-10,10);  // Intercept of entire population
    z ~ normal( 0 ,1); //Intercept of each cluster
    et=et_f(timeck,logdelta,M,N);  //
    for ( i in 1:M ) {
      //first calculate the log-rate of being infected
      p[i]  = loga_bar+z[isol[i]] - logdelta + log((1-exp(-exp(logdelta)*dt[i]))*et[i] + i0[i]*(((exp(-exp(logdelta)*dt[i])-1)/exp(logdelta))+dt[i]));
      p[i] = 1-exp(-exp(p[i]));            
    }
   ca ~ binomial( s0 , p );
}
2

Theoretically p<1 as long as dt>0 but when the calculation has multiple nested log() and exp() intermediate results may overflow. It’s better to use the composed functions expm1() and log1m_exp() whenever possible. Likewise, logit-parameterized binomial is more stable than the usual parameterization.

model {
    vector[M] p;
    vector[M] et;
    loga_bar ~ normal(-10,10);  // Intercept of entire population
    z ~ normal( 0 ,1); //Intercept of each cluster
    et=et_f(timeck,logdelta,M,N);  //
    for ( i in 1:M ) {
      //first calculate the log-rate of being infected
      p[i]  = loga_bar+z[isol[i]] - logdelta + log(-exp1m(-exp(logdelta)*dt[i])*et[i] + i0[i]*((exp1m(-exp(logdelta)*dt[i])/exp(logdelta))+dt[i]));
      p[i] = exp(p[i]) + log1m_exp(-exp(p[i])); // logit of the probability variable
    }
    ca ~ binomial_logit( s0 , p );
}

If that still doesn’t work, or you want upper limit 0.99 specifically, there’s not much other option than

      p[i] = 0.99*p[i];

after all other transforms.

3

Dear @nhuurre Thank you very much!. I tried both methods and it work fine.