computing transition probabilites given only aggregate counts

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.

• A request for clarification: (1) The example matrix from the Chicago paper appears to show that M_{1,2}=19. However, I understood that by construction M_{t,j}=0 for all j > t. (2) The reported OLS estimates include p_{2,2} = 0.823. But I had understood that p{2,2} = 1.00 by construction, since State 2 is an absorbing state. – Arthur Small Dec 5 '12 at 4:52
• re (1)---the data for the chicago paper doesn't have the upper-triangulare constraint. re (2)---that paper also uses 1-based indexing for the states, so that's an inconsistency in my presentation; will fix. – ben w Dec 5 '12 at 16:08

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$.