 # 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!