I’m trying to fit a mixture of transition matrices but I’m a bit stuck on how I can “order” the k matrices. I’ve basically just adapted the normal mixture from the documentation here but the ordered part is tripping me up. I’m sure this not the optimal model specification, but I could converge on a k=2 model with fake data but getting issues beyond that. For example, the data below will put all the data in one cluster.

```
// Pathways model -- Finite mixture of transition matrices
data {
int<lower=1> K; // the number of clusters
int<lower=1> N; // the number of datapoints (traversals)
int<lower=1> S; // the number of observable states
int<lower=0> X[N, S, S]; // the 3d-array of transitions (members * to-state * from-state)
}
parameters {
simplex[K] theta; // mixture distribution
simplex[S] transitions[K, S]; //transition matrix is array of simplices
}
model {
vector[K] log_theta = log(theta); //cache log calculation
for (n in 1:N) {
vector[K] lps = log_theta;
for(k in 1:K){
for (s in 1:S){
lps[k] += multinomial_lpmf(X[n,:,s]| transitions[k, s]);
}
}
// lps[k] += normal_lpdf(y[n]| mu[k], sigma[k]);
target += log_sum_exp(lps) ;
}
}
```

Here’s the fake data (R):

```
##Input data
n_clusters <- 3
n_states <- 6
mixture_dist <- c(.6, .20, .20)
N <- 10000
#cluster 1 -- middle (same vector for all states for simplicity)
p1 <- rep(1/n_states, n_states)
#Cluster 2 - increasing
inc <- 1:n_states/n_states
p2 <- inc/sum(inc)
#Cluster 3 - decreasing
p3 <- rev(p2)
X <- array(, dim=c(N, n_states, n_states))
for(i in 1:6000){
X[i, ,] <- rmultinom(6, 100, prob= p1)
}
for(i in 6001:8000){
X[i, ,] <- rmultinom(6, 100, prob= p2)
}
for(i in 8001:10000){
X[i, ,] <- rmultinom(6, 100, prob= p3)
}
##############
input_data <- list(
K=n_clusters,
N=N,
S=n_states,
X=X
)
stanmod1 <-
stan_model('pathways.stan',
model_name = 'sample'
)
fit1 <-
optimizing(
stanmod1, data=input_data
)
```