# Optimal decision process to estimate Markov chain limiting distribution

Suppose there is a irreducible, reversible Markov chain with known states $1,\ldots,N$ and unknown transition matrix $T_{ij}$ and unknown limiting distribution $\pi_i$. I am able to repeatedly initialize the system in an arbitrary state $s$ then observe a random transition $s\rightarrow t$ where $t$ is a sample from $P(S_{k+1} | S_k=s)$ in the Markov chain. I would like an algorithm to choose the initial states $s_k$ dependent on the previous results $(s_1 \rightarrow t_1), \ldots, (s_{k-1} \rightarrow t_{k-1})$, to optimally estimate $\log \pi_2 - \log \pi_1$.

A naive method would be to choose $s_k = t_{k-1}$ to obtain a sample of length $M$ from the Markov chain. Then, the estimator $\log \frac{1+\sum_k (s_k = 2)}{1+M} - \log \frac{1+\sum_k (s_k = 1)}{1+M}$ would converge to $\log \pi_2 - \log \pi_1$ as $M \rightarrow \infty$. This method is inefficient when either $\pi_1$ or $\pi_2$ is very small, which is the case I am interested in.

Is there a better, or even optimal, scheme to choose $s_k$ dependent on $(s_1 \rightarrow t_1), \ldots, (s_{k-1} \rightarrow t_{k-1})$ to estimate $\log \pi_2 - \log \pi_1$ for a fixed number of observed transitions $M$?

## Notation / preliminaries

I use hats ($\hat \theta$) for estimates, except when I say otherwise or forget. I try to use OP's notation: $N$ states, tpm $T$, and stationary dist $\pi$.

The maximum likelihood estimate of $T_{i\rightarrow j}$ is $m_{i\rightarrow j} / m_i$, the number of $i$-to-$j$ transitions divided by the number of transitions from state $i$.

## A more flexible estimation tactic

Instead of letting the chain run forever, another way to estimate the stationary distribution is to estimate the transition matrix and use its leading eigenvector (scaled to sum to 1). Indeed, one way to calculate the eigenvector is via the power method: repeat $v^T\leftarrow v^T \hat T$. Almost any initial guess will work, but if you use $v_0 = [1, 0, 0, 0,...]$, then the your strategy and mine start to look very similar. You run the chain forever, starting from state 1. I calculate the successive marginal distributions sampled by your strategy -- but with errors, because my estimate $\hat T$ won't be perfect.

Estimating the TPM and using the top eigenvector gives you the freedom to sample from the states where the transition probabilities are most uncertain (rare states) or most useful (states that feed towards your target states). You can also initialize flexibly: for example, you can start 100 times from every state, rather than letting the chain drag you around to all the stupid/popular/touristy states. You can use an advantageous estimation method if you think the TPM has useful structure (toeplitz, banded, low rank, nonzero, smooth...).

The question then becomes:

• What's the expected information gain from a single sample starting in state $i$? Or at least, what's the argmax?
• How quickly can I compute this argmax relative to taking a single sample? If it's slow, then I'd compute (say) 10000 draws starting from state argmax(infogain) before recomputing argmax(infogain).

## Defining information gain

Suppose for the moment that $T$ is known. I will define information gain in terms of $T$, then figure out what to plug in later.

Your estimand can be written as $\theta(T)$. I'll define the information as

$$I_\theta (m_1, m_2, ...) = Var [\theta( \hat T_{m_1, m_2, ...})] ^{-1}$$, where $\hat T_{m_1, m_2, ...}$ is a random matrix obtained by sampling $m_i$ times from state $i$. Then I'll define the information gain with batch size $B$ as $$IG_\theta(i, B) = I_\theta (m_1,..., m_i + B, ...) - I_\theta (m_1, m_2, ...)$$.

This suggests a plug-in estimation strategy for $IG(i, B)$. Suppose you've sampled $m_1, m_2, ...$ times. Here's R code to estimate the information gain. This is slow; I recommend you try to re-use samples and vectorize everything.

info = function( to_sample = 1, batch_size = 10000, boot_niter = 1000){
for( n_bootstrap in 1:boot_niter ){
for( from_i in 1:N ){
for( to_j in 1:N ){
size = m[i]
if( from_i == to_sample_next ){ size = size + batch_size }
T_boot[ from_i, to_j ] = rbinom(size = size,
prob = T_current[from_i, to_j]) / size
}
}
limit_dist = T_boot[1, ]
# Get eigenvector
for( power_method in 1:10){
limit_dist = limit_dist %*% T_boot
}
your_parameter[n_boostrap] = log( limit_dist[1] ) - log( limit_dist[2] )
}
}
return( 1 / var( your_parameter ) )
}
info_gain = function( to_sample_next ){
info( to_sample_next = to_sample_next, batch_size = 1000) -
info( to_sample_next = 1,              batch_size = 0)
}


If you choose to use this, it's crucial to avoid getting stuck in traps. For example, due to sampling variability, you could end up with a TPM that has multiple disconnected regions. If one of these doesn't contain your target states, then sampling more from it will never decrease the bootstrap variance of your target parameter, and the info gain will never tell you to look there. I suggest you add pseudo-counts to the MLE or start with a non-fancy strategy until your estimated TPM has enough nonzeroes for any state to reach any other.

Also, I don't actually know how to bootstrap properly, and I think my variance estimate might be biased low or otherwise shitty. Maybe you should use the parameter estimate from T_current instead of the bootstrap mean (which I used implicitly in my call to var).

## Without resampling

Another option is to use the delta method:

$$sd(\theta) \approx \sum_{i,j,k}\frac{d\theta}{d\pi_k}\frac{d\pi_k}{dT_{ij}}sd(T_{ij})$$ (I've omitted the hats, but everything here is an estimate, not the ground truth.)

The standard deviation of the maximum likelihood estimator given above is roughly $\sqrt{ m_i^{-1}\hat T_{ij} (1-\hat T_{ij}) }$ (plugging in the estimate $\hat T_{ij}$ as you would for a simple binomial distribution). With $B$ additional samples, one would expect roughly $\sqrt{ (B+m_i)^{-1}\hat T_{ij} (1-\hat T_{ij}) }$.

The obstacle is $\frac{d\pi_k}{dT_{ij}}$, which is the derivative of an eigenvector with respect to the matrix it came from. I don't know how to compute that, but it's definitely possible.

• Thanks for a great answer, especially to such an old question. Your answer reminds me of Bayesian hyperparameter optimization. Since my space in a comment is limited, I added another answer where I tried to work out a Bayesian version. I think it makes some of the technical problems go away. – John Jumper Mar 27 '17 at 4:45
• My pleasure! There definitely is a possible Bayesian approach lurking here. If you have a (conjugate) Dirichlet prior on the columns of the TPM, the posterior is also Dirichlet, and the parameters are a simple function of the prior params plus counts. Then you can somehow propagate uncertainty from the posterior over the TPM to the posterior over the limiting distribution. – eric_kernfeld Mar 27 '17 at 18:19
• Yes, although I am not terribly worried about working with conjugate distributions. I have had a lot of success before using Stan to sample complex posteriors, and I may want the freedom on the prior. My application is in protein molecular simulation, where sampling a single transition can take minutes to hours, so time spent picking the starting state for the simulation is not important. – John Jumper Mar 27 '17 at 19:39
• Yikes! Wow! Is there any sense in modeling the simulation times as a function of the starting state and oversampling the states where the computation is minutes instead of hours? You could measure information gain per hour instead of information gain per sample. – eric_kernfeld Mar 27 '17 at 20:38
• That's certainly an interesting thought. In some sense, the expected transition matrix over the prior already incorporates that information. If self-transitions are very high probability, they would likely have low expected information gain. Ultimately, parallelism helps as you can sample many transitions in parallel on cluster (see folding.stanford.edu/home/faq/faq-simulation for an example of Markov modeling in a simulation context). – John Jumper Mar 27 '17 at 20:44

Suppose we have a prior $p(T)$ on transition matrices and a sample $\{T^{(k)}\}_{k=1,\ldots,K}$ from that prior. We can estimate the prior-averaged transition matrix as $\bar{T} = \frac{1}{K} \sum_k T^{(k)}$. The entry $\bar{T}_{s\rightarrow t}$ is the marginal probability of observing state $t$ for a single transition starting from state $s$ with a transition matrix drawn from the prior distribution $p$.
Let $f$ be some scalar function of the transition matrix, where $f(T) = \log \pi_2(T) - \log \pi_1(T)$ in the original question. The information gain of observing a transition $s\rightarrow t$ can be defined as $\operatorname{KL}(q_{s\rightarrow t}(f),p(f))$ where $q_{s \rightarrow t}(T)$ is defined as the posterior probability of T after observing the transition $s\rightarrow t$.
We can estimate the prior distribution of $f$ using a kernel density estimate on $\{f(T^{(k)})\}_{k=1,\ldots,K}$. To approximate the posterior distribution on $T$, we can use a reweighting estimate. Define the reweighted probability $q^{(k)}_{s\rightarrow t} = T^{(k)}_{s\rightarrow t}/K \bar{T}_{s\rightarrow t}$. The information gain for observing $s \rightarrow t$ is just $$I_{s\rightarrow t} = \operatorname{KL}(q_{s\rightarrow t}(f),p(f)) \approx \sum_k q^{(k)}_{s\rightarrow t} \log \left(\frac{\sum_{k'} q^{(k')}_{s\rightarrow t} \operatorname{kernel}(f(T^{(k)}),f(T^{(k')}))}{ \sum_{k'} \frac{1}{K}\operatorname{kernel}(f(T^{(k)}),f(T^{(k')}))} \right)$$
The expected information gain from observing a single transition from state $s$ is $I_s = \sum_t \bar{T}_{s\rightarrow t} I_{s\rightarrow t}$. We can choose the initial state $s$ to maximize the information gain. The process can be iterated by replacing the prior distribution with the new posterior after sampling $s\rightarrow t$.