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

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.

``````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
}
``````

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;
}
``````

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

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;
}
``````

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

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

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

1 Like

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

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!