I wonder if using reject statement in my STAN code would bias the estimates. The problem is that I analyze an interval-based data composed of pairs. Each pair is characterized by two symptom onsets o_1, and o_2, such that o_1 < o_2. Because the pairs are not independent and the intervals in the data are overlapping (I mean, O_{1L} < o_1 < O_{1R} and O_{2L} < o_2 < O_{2R}, but O_{2L} can be earlier than O_{1R}). It is difficult to impose restrictions directly and I use the reject statement:
if (o_1 >= o_2)
reject("Onset of the infector should not be later than onset
for the infectee, violated for pair ", n, "; onset1 = ", o_1, ";
onset2 = ", o_2);
Because of that, I would like to know whether it is bad practice and should not be used (for example, the posteriors would be biased), or it is kind of acceptible given the circumstances. I know that I would sacrifice at least the efficiency.
There may be a better way to parameterize your constraints, but here’s one way:
Marginalize over 3 possibillities. For notational simplicity, I omit the simple upper and lower bounds constraints describe above for each parameter.
o_1 < O_{2L}. In this rectangular region, we can parameterize simply with upper and lower bound constraints on each element.
o_1 > O_{2L} \quad \textrm{and} \quad o_2 > O_{1R}. In this rectangular region, we can again parameterize simply with upper and lower bound constraints on each element.
o_1 > O_{2L} \quad \textrm{and} \quad o_2 < O_{1R} \quad \textrm{and} \quad o_1 < o_2. In this right-triangular region, we can parameterize by declaring an ordered vector, transforming both elements into the unit interval using the upper-and-lower-bounds transform, and then shifting and scaling the triangle into position.
Yes, it is bad practice and should not be used. The fundamental problem with the rejection approach is that it causes Hamiltonian Monte Carlo to devolve to rejection sampling, which can be arbitrarily inefficient depending on how close you are to the constraint surface and how likely random jumps are to take you into violation.
if o1 and o2 are parameters, then just do this
parameters {
real o1;
real<upper=o1> o2;
}
This does not force the posterior interval for o1 to be strictly lower than the posterior interval for o2. For example, consider this simple program.
parameters {
real o1;
real<upper=o1> o2;
}
model {
o1 ~ normal(0, 1);
o2 ~ normal(0, 1);
}
Here is the relevant part of the posterior summary:
The posterior intervals are (-0.7, 2.0) and (-2.0, 0.8). These are what are called “order statistics” for the normal. You get the same result from drawing two standard normals and sorting them.
“Stan” is not an acronym—we named the system after Stanislaw Ulam.
I think the question is how to do this while also including separate fixed upper and lower bound constraints on o1 and o2. I’m unsure whether it would break anything to write real<upper = min(o1, o2L)> o2. Maybe this works out of the box (is it ok that the bound isn’t continuously differentiable?).
In the event that this doesn’t work, it’s possible to chuck this up into very straightforward shapes (2 rectangles and a triangle) and marginalize over those.
To provide a fully worked example, I’d use something like this if you have constraints lb1 <= o1 <= ub1, lb2 <= o2 <= ub2 and o2 < o1, you can just write:
The discontinuity of the derivatives at the point ub2 == o1 shouldn’t be too bad. Usually things are fine if there isn’t a huge energy barrier that would prevent sampling across the boundary. We do this kind of thing all the time with Stan.
Thank you very much for your replies, Jacob and Bob! Indeed I had similar thoughts regarding the rejection statement. I formulated my example in a bit simpler way, although, the main difficulty was the dependence of the pairs. Each pair is composed of individual 1 and individual 2, but individual 1 can be sometimes as individual 2 in some other pairs.