I want to fit a mixture model hierarchically but have hit a wall. Currently, I have the parameters of a Dirichlet distribution describing the probability for a model. The 4 clusters that are thought to compose my data are 3 von Mises distributions (one centered at mu, one centered at swapLoc, and one centered at swapCLoc), and 1 uniform guessing distribution. The problem I have is that I want to get both individual-level and group-level estimates for these mixing proportions. The below code is a working non-hierarchical mixture model (the transformed parameters section I took from here). Note that all the data is is circular (in radians), which is why I use von Mises distributions.
I assume to implement hierarchy, I need to have different p[c]'s for each individual and v needs to be v ~ beta(a,b) instead of beta(1,1), where a and b are drawn from some larger distribution. But turning p[c] into a matrix like p[i,c] has presented me with a flurry of errors, and I’m not sure how to specify this larger distribution for v ~ beta(a,b) (if that’s even what I need to do to get this working). I am at a loss for how to proceed. All help appreciated.
data {
int NumSubj; //44
int NumData; //9635
int C; //num of clusters
vector[NumData] ID; // subject ID (1-44)
int NumResp[NumSubj]; // subject's num of responses
int CumResp[NumSubj]; // cumulative num of responses
vector[NumData] errors;
vector[NumData] swapLoc;
vector[NumData] swapCLoc;
}
parameters {
real<lower=-pi()/6,upper=pi()/6> mu;
real<lower=0,upper=100> kappa;
real<lower=0,upper=100> kappa2;
real<lower=0,upper=1> v[C];
}
transformed parameters {
simplex[C] p; // mixing proportions for each of C clusters (sums to 1)
p[1] = v[1];
// stick-breaking Dirichlet Process based on The BUGS book Chapter 11 (p.294)
for (l in 2:(C-1)) {
p[l]= v[l]*(1-v[l-1])*p[l-1]/v[l-1];
}
p[C]=1-sum(p[1:(C-1)]); // to make a simplex
}
model {
real ps[C];
v ~ beta(1,1);
kappa ~ normal(10,100);
kappa2 ~ normal(10,100);
mu ~ normal(0,5);
for (i in 1:NumSubj) {
for (j in 1:NumResp[i]) {
for (c in 1:C) {
if (c == 1)
ps[c] = log(p[c])+von_mises_lpdf(errors[CumResp[i]-(j-1)] | swapLoc[CumResp[i]-(j-1)], kappa2); // Swap
else if (c == 2)
ps[c] = log(p[c])+von_mises_lpdf(errors[CumResp[i]-(j-1)] | swapCLoc[CumResp[i]-(j-1)], kappa2); // Swap-Comparison
else if (c == 3)
ps[c] = log(p[c])+uniform_lpdf(errors[CumResp[i]-(j-1)] | -pi(), pi()); // Guess
else
ps[c] = log(p[c])+von_mises_lpdf(errors[CumResp[i]-(j-1)] | mu, kappa); // Target
}
target += log_sum_exp(ps);
}
}
}