Choosing the right likelihood for a discrete dynamical system

Modeling
1

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.

2

If I understand correctly, the determinism implied here is not desired, and what you want instead is

a_to_b[i] ~ Binomial(A[i], r1)
b_to_c[i] ~ Binomial(B[i], r2)
leave_c[i] ~ Binomial(C[i], r3 + r4)
c_to_a[i] ~ Binomial(leave_c[i], r4/(r3 + r4))

subject to the constraint that r3 + r4 <= 1. At least I’m pretty sure this is what you want.

The obvious challenge here is that none of the left-hand-sides is ever observed, so the question is whether there is an efficient way to marginalize over all the possibilities. Maybe @Bob_Carpenter can spot a useful marginalization trick here.

1 Like