I’ve simulated some data from the following discrete dynamical system

A_{n+1} = (1-r_1)A_n + r_4C_n

B_{n+1} = (1-r_2)B_n + r_1A_n + r_3 C_n

C_{n+1} = (1-r_3-r_4)C_n +r_2 B_n

This can be more compactly written as

\mathbf{x}_{n+1} = T \mathbf{x}_n = T ^n\mathbf{x}_0

Here, T is a matrix containing the transition rates. The columns of T sum to 1 and since the next state only depends on the current state I *believe* this constitutes a Markov chain.

My question concerns an appropriate likelihood for this model. As of now, I use a multinomial likelihood, where the multinomial parameter, \theta, is obtained through the matrix multiplication above at each time step. Below is my current Stan model

```
data{
//Number of timesteps
int n;
// Number of states
int p;
//Matrix of state observations, one column per state
array[n, p] int X;
}
transformed data {
//Initial state observation
int X0[p] = X[1, :];
simplex[p] theta0 = 1.0 .* to_vector(X0) ./ sum(X0);
}
parameters {
real<lower=0, upper=1> r[4];
}
transformed parameters {
matrix[p, p] T = rep_matrix(0.0, p, p);
array[n] simplex[p] theta;
theta[1] = theta0;
T[1, 1] = 1.0 - r[1];
T[2, 1] = r[1];
T[2, 2] = 1-r[2];
T[3, 2] = r[2];
T[1, 3] = r[4];
T[2, 3] = r[3];
T[3, 3] = 1.0 - r[3] - r[4];
for(i in 2:n){
theta[i] = T*theta[i-1];
}
}
model {
r ~ beta(10, 90);
for(i in 2:n){
target += multinomial_lpmf(X[i] |theta[i]);
}
}
generated quantities {
array[n, p] int Xppc;
Xppc[1] = X0;
for(i in 2:n){
Xppc[i, 1:p] = multinomial_rng(theta[i], sum(X0));
}
}
```

While the multinomial does a good job of estimating the rates, it isn’t the right likelihood to use. The randomness is in *how many units in each state leave the state*, and this must be an integer since the dynamics are discrete. For example:

- Suppose there are 100 units in C
- Units leave C to go to A at a rate of r_4=0.01 and leave C to go to B at a rate of r_3=0.03.
- Units enter C from B at a rate of r_2=0.02.
- So the randomness is in determining how many units are going to leave the state, as opposed to randomness from measurement error.
- To determine how many leave, I draw from a multinomial to determine
*where units are headed*. - Here is a python script to simulate the process as I understand it working.

The result is that the multinomial likelihood is less autororrelated than the processes I simulate.Here is a simulation from the processes, and here is a draw from the posterior predictive.

How can I better specify my likelihood so that my posterior predictive is more autocorrelated? The answer might be “its hard” or “you can’t”, which is fine.