I have a situation where data from the following process is observed:
For $i = 1, \dots, n$ let $(X_{i,1}, \dots, X_{i,m_i})$ be a sequence of $m_i$ random variables coming from a discrete-space Markov chain with $N \times N$ transition matrix $P$. For simplicity's sake I'll assume that each $X_{ij} \in \{1, \dots, N\}$. For $i \neq i'$ assume that the chains $(X_{i,1}, \dots, X_{i,m_i})$ and $(X_{i',1}, \dots, X_{i',m_{i'}})$ are independent of each other and both have the same transition matrix $P$. In general $m_i \neq m_{i'}$.
A big complication: $m_i$ is really small. The largest $m_i$ is 9, the smallest is 3, and about 95% of the $m_i = 5$. I have around 50 of these chains (i.e. $n = 50$).
For a single chain the MLE $\hat P_i$ of $P$ is easy to compute: it's just the proportions of times that each transition was observed.
My question: how to get a single estimate of $P$ from these $\hat P_i$? Each $\hat P_i$ estimates the same thing, and due to differences in $m_i$ some do better than others. But because each $m_i$ is small, none of them do very well. There are lots of 0's and 1's which I don't think reflect reality. What should I do here?