Dirichlet-Multinomial for Transition Matrix Estimation

Dear Stan users,
I’m attempting to model my 3-by-3 transition matrix (discrete-time Markov chain) using a Dirichlet-Multinomial conjugate pair system. My data is a “raw” transition matrix: that is, instead of a given entry in the matrix representing the probability of a transition, it represents the number of counts of a transition. As such, the likelihood is inherently trinomial (not a categorical, because order doesn’t matter thanks to the Markov property). However, the likelihood calculation and generation of a posterior Dirichlet (a K=3 simplex) needs to occur three times per iteration — once for each row of my matrix.
So, I need a way to generate a Dirichlet prior (that doesn’t depend on my raw transition matrix), while still having my multinomial likelihood concentrate the information presented in the data into the single “theta” parameter that the multinomial sampling statement affords (14.1 in manual). I’m not sure which block to declare my Dirichlet prior (generated parameters or model?), and which block to define my single Multinomial parameter. I’m also very new to matrix representation in Stan, so I apologize for the presence of syntactical errors. Any tips on this would be greatly appreciated by this incipient Stan user!
Cheers


data {
	// Prior Parameters
   int<lower = 1> K; 		// COLS. for Dirichlet & txn mat
   int<lower = 1> N; 		// ROWS. for txn mat
   real<lower = 0> alpha; 	// for Dirichlet

        // Experimental Parameters
   matrix[N,K] T; 		// input RAW txn matrix, the data
 }

generated quantities { 
   vector[K] theta = dirichlet_rng(rep_vector(alpha, K));  // prior
 }

parameters {
  matrix[N,K] P;         // parameter of interest
}

model {
     theta ~ dirichlet(alpha);	// prior?
    for(i in 1:N) {
    P[i,] ~ multinomial T[i,];        // likelihood
    }
 }

Sorry, your question fell through a bit. I don’t immediately see what is the problem with your model, but if I understand it correctly, you want P to be a transition matrix. In that case, I think you should define it as an array of simplices (as that enforces that transition probabilities sum to 1). Additionally, T is counts, so probably should be int. Here is a rough sketch (not tested) how I would code the model:

data {
	// Prior Parameters
   int<lower = 1> K; 		// COLS. for Dirichlet & txn mat
   int<lower = 1> N; 		// ROWS. for txn mat
   real<lower = 0> alpha; 	// for Dirichlet

        // Experimental Parameters
   int T[N,K]; 		// input RAW txn matrix, the data
 }

parameters {
   simplex[K] P[N]; 
}

model {
  for(i in 1:N) {
    P[i] ~ dirichlet(rep_vector(alpha, K)); // Prior, Dirichlet with all arguments the same
    T[i,] ~ multinomial(P[i,]);        // likelihood
  }
}

Hope that clarifies some of the stuff and best of luck with your model!