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$?
 A: 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.
A: This is an elaboration of Eric Kernfeld's answer.  Please upvote his instead of this one.
The information gain optimization of the problem is particularly natural in a Bayesian formulation.  The Bayesian formulation also avoids some of the problems with either invariant subspaces or bootstrapping the variance in the plugin estimator.
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$.
