computing transition probabilites given only aggregate counts

Question

I asked something like this on math.stackexchange but it got no answers so I'm hoping for more enlightenment here (the other question is this one).

So: We have some number of people; each person i is in a state $X_i(t)$ at time $t$ and can move at time $t+1$ into any state $0…X_i(t)+1$ inclusive (unless they are in the final, accepting, state). At time $t=0$ everyone is in state $0$.

The individual state transitions are unobserved; all we observe are the counts of people in a given state $j$ at time $t$: $\mathbf{M}_{t,j}=\#\{i:X_i(t)=j\}$. Because everyone starts in state 0 at time 0, and because you can only move to at most state $k+1$ if you were formerly in state $k$, this generates an upper-triangular matrix of counts. The last state is accepting.

We want to know the probabilities of the possible state transitions $p_{k,j}=\Pr(X_i(t+1)=j∣X_i(t)=k)$ given the counts in $\mathbf{M}$.

I found a paper (here) about basically exactly this; this other report has some more info.

That first paper contains an example with this data:

$\begin{bmatrix}200 & 123 & 73 & 42 & \\ 0 & 58 & 84 & 90 & \\ 0 & 19 & 43 &68\end{bmatrix}$

And claims that (using OLS) $p_{00} = 0.606$, $p_{01} = 0.291$, $p_{11} = 0.823$. I managed to hack something in octave (and now have access to matlab) that gives this:

function ans = transitions(data)
     data = data / data(1,1)
     xs = data(:,1:columns(data)-1)'
     p = []
     for i = 1 : rows(data) -1
         p2 = lsqnonneg(xs, data(i, 2:columns(data))')
         p = [p, p2]
     endfor
     ans = p
endfunction

BUT function gives totally nuts results for the actual data (in particular the diagonals are screwy and the last state comes out to be not even close to accepting). I attempted to sanity-check it using two test matrices that I generated:

$\mathbf{T}_1 = \begin{bmatrix}1000 & 900 & 815 & 742.75 \\ 0 & 100 & 180 & 244 \\ 0 & 0 & 5 & 13 \\ 0 & 0 & 0 & .25\end{bmatrix}$

$\mathbf{T}_2 = \begin{bmatrix} 1000 & 903 & 818 & 745 \\ 0 & 97 & 177 & 241 \\ 0 & 0 & 5 & 13 \\ 0 & 0 & 0 & 1\end{bmatrix}$

The first one was just multiplied straight through by the pre-decided transition probability matrix, the second one deviates from that. The function gives back the correct transition probabilities for the first test matrix but, again, very odd (and very different) results for the second.

I don't have anything close to the knowledge necessary to directly implement what's in either of those papers, so any help on how to do this would be greatly appreciated.

Answer 1

This is not a complete answer, but I'll offer a couple suggestions.

First, you might try breaking your problem down into pieces, and solving iteratively, perhaps using a Bayesian updating procedure.

For example, in your test data set $T_1$, you have $M_{0,0} = 1000$, $M_{1,0} = 900$, and $M_{1,1} = 100$. Focusing just on this subset of the data, a MLE procedure will give you estimates of $p_{0,0} \approx 0.9$ and $p_{0,1} = 1-p_{0,0} \approx 0.1$.

By construction you already know that $p_{0,2} = p_{2,0} = p_{2,1} = 0$ and $p_{2,2}=1$. So you've got a lot of structure that should help reduce the dimensionality of the problem.

Expanding on this theme, you could estimate these transition probabilities using a Bayesian procedure, using a Beta distribution for the prior on $p_{0,0}$. Using only the data from the transition from $t=0$ to $t=1$, your posterior distribution over $p_{0,0}$ would be a Beta distribution with mean $0.9$. This posterior distribution could then be used as the prior for the next round, when you need to analyze data reflecting the transition from $t=1$ to $t=2$.

A final thought: you might look at maximum entropy estimation procedures. They tend to work very nicely when you need to estimate probabilities for under-determined problems.