I have been working on a model with two success probabilities and a change point. The model recovers well when I marginalize out the change point as a latent discrete parameter. However, in the final model I would like it if I don’t have to assume one population change point, but rather change points that can vary between subjects. I have read some examples of random change point models, but the one implementation I found in Stan does not model the change point as a discrete parameter. This brings me to my question: is it possible to have a hierarchical discrete latent parameter or can it only be hierarchical if the change point is modeled as a continuous parameter?

```
// data:
data{
int<lower=1> L; // number of participants
int<lower=1> N; // number of datapoints
int<lower=0,upper=1> Y[N]; // outcome (cancelled or not)
int<lower=1,upper=L> ll[N]; // cases going from int 1 to L
vector[N] X; // X-coordinates
}
// transformed data: (necessary for latent discrete parameter)
transformed data{
real boundary[7];
// possible boundary values are -0.75, -0.50, -0.25, 0, 0.25, 0.50, 0.75
for(i in 1:7){
boundary[i] = -0.75 + (0.25 * (i-1));
}
}
// parameters:
parameters{
// population-level theta's in log odds:
real theta_left_logit;
real theta_right_logit;
// participant-level theta's in log odds:
vector[L] theta_left_l;
vector[L] theta_right_l;
// population sd of theta by participant:
real<lower=0> sigma_l;
real<lower=0> sigma_r;
}
// transformed parameter:
// RECORD Log probability FOR EACH BOUNDARY LEVEL:
transformed parameters{
vector[7] lp; // lp = log probability
lp = rep_vector(-log(7), 7); // uniform prior on all boundary values
// for each possible boundary:
for(b in 1:7){
for(n in 1:N){
lp[b] = lp[b] + bernoulli_logit_lpmf(Y[n]|X[n]<boundary[b] ? theta_left_l[ll[n]] : theta_right_l[ll[n]]);
}
}
}
// model:
model{
// priors for population level: (hyperpriors)
theta_left_logit ~ normal(0, 2);
theta_right_logit ~ normal(0, 2);
// priors for population-level sigma
// (by lb they exclude values lower than zero):
sigma_l ~ normal(0.05, 0.025);
sigma_r ~ normal(0.05, 0.025);
// priors for participant level parameters:
// should be drawn with population parameter as the mean:
// and population sigma as sd:
theta_left_l ~ normal(theta_left_logit, sigma_l);
theta_right_l ~ normal(theta_right_logit, sigma_r);
target += log_sum_exp(lp);
}
// extra code that allows to infer the probability of the discrete parameter:
generated quantities{
vector[7] pState;
pState = exp(lp - log_sum_exp(lp));
}
```