How to set initial value to make probability parameter be in the interval [0, 1] in a stan model

Modeling
,
1

Dear stan community. I would appreciate any help to solve the problem that is causing the error below detailed in my stan model:

#error

SAMPLING FOR MODEL ‘m’ NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: binomial_lpmf: Probability parameter is 1.59857, but must be in the interval [0, 1] (in ‘model41bc3e0d5412_m’ at line 33)

#line 33

ncases ~ binomial(nn, pp_hat);

#model

data{
int N;  
int ncases[N];
int A[N]; 
int B[N];  
int nn[N]; 
}

parameters {
real alpha_0; 
real alpha_1;
real alpha_2; 
real alpha_3; 
}

transformed parameters {
vector[N] pp_hat;   
for (i in 1:N) {    
pp_hat[i] = exp(alpha_0) * (1 + alpha_1*A[i] + alpha_2*B[i] + alpha_3*A[i]*B[i]);
}
}

model {
alpha_0 ~ normal(0, 5);
alpha_1 ~ normal(0, 5);
alpha_2 ~ normal(0, 5);
alpha_3 ~ normal(0, 5);
ncases ~ binomial(nn, pp_hat);
}

I am not certain how to apply the ideas discussed here in my case in particular. Thank you in advance for any help.

2 
data {
  int N;
  /*  int ncases[N]; is not used */
  vector[N] A;
  vector[N] B;
  int nn[N];
}
transformed data {
  vector[N] AB = A .* B;
}
parameters {
  vector[3] alpha;
  real<lower = 0, upper = inv(max(1 + alpha[1] * A + alpha[2] * B + alpha[3] * AB))> 
    exp_alpha_0;
}
model {
  vector[N] right = 1 + alpha[1] * A + alpha[2] * B + alpha[3] * AB;
  real bound = inv(max(right));
  vector[N] pp_hat = exp_alpha_0 * right;
  ncases ~ binomial(nn, pp_hat);
  alpha ~ normal(0, 5);
  exp_alpha_0 ~ lognormal(0, 5) T[0, bound]; 
  // this prior is consistent with what you wrote but dubious given the constraints
}
3

Great. Thank you very much.

If you don’t mind, what prior would be more consistent given the constraints?

#inequality constraints for model 1

(1) alpha_0 ≤0;
(2) exp(- alpha_0) ≥ 1 + alpha_1 ≥0;
(3) exp(- alpha_0) ≥ 1 + alpha_2 ≥ 0; and
(4) exp(- alpha_0) ≥ 1 + alpha_1 + alpha_2 + alpha_3 ≥0

…

parameters {
  real<upper = 0> alpha_0;
  real<lower = -1, upper = exp(-alpha_0)> alpha_1;
  real<lower = -1, upper = exp(-alpha_0)> alpha_2;
  real<lower = -1 - alpha_1 - alpha_2, upper = exp(-alpha_0)> alpha_3;
}
4

I would start by just letting it be uniform on that bounded interval. Then maybe beta over that interval.

5

This is very helpful. Thanks for that @bgoodri. I will use uniform/beta on that interval.

Actually, I am getting the same error for the model 2. I thought it could be straightforward for me to translate the solution from model 1 to model 2, but it seems that I might still need your help to solve the error for model 2 too given the constraints and other considerations:

#inequality constraints for model 2

(1) 1 + beta_1 > 0;
(2) 1 + beta_2 > 0; and
(3) 1 + beta_1 + beta_2 + beta_3 > 0

parameters {
  real<lower = -1> beta_1;
  real<lower = -1> beta_2;
  real<lower = -1 - beta_1 - beta_2> beta_3;
}
6

I would reparameterize beta_0 into exp_beta_0 and do the same thing as before.

7

Thanks, @bgoodri . Let me see and will let you know my progress on model 2.

8

You just have to constrain things tightly enough so that your probabilities are between 0 and 1.

1 Like
9

Thanks, @bgoodri. I truly appreciate your help. Great!

11

What would be the change if the original pp_hat[i] = exp(alpha_0) * (1 + alpha_1A[i] + alpha_2B[i] + alpha_3*A[i]*B[i])

was actually

pp_hat[i] = exp(alpha_0 + alpha_1A[i] + alpha_2B[i] + alpha_3*A[i]*B[i])

thanks!