Markov chain simulation I'm wondering whether there is an algorithm to simulate a discrete Markov chain with a specific number of occurrences of state knowing the transition matrix way.
For example, how to simulate in R a Markov chain of length $n$ with $p$ occurrences ($p < n$)
TransitionMatrix<- matrix(c(0.7, 0.3, 0.4, 0.6),byrow=TRUE, nrow=2)
colnames(TransitionMatrix) <- c('0','1') row.names(TransitionMatrix) <- c('0','1')

 A: First a reformulation of the question, then a pseudo-algorithm to solve it.
Reformulation: Let $n\geqslant1$ and assume that $X=(X_t)_{1\leqslant t\leqslant n}$ is a Markov chain of length $n$ on the state space $\{0,1\}$ with transition probability matrix $\begin{pmatrix}1-a & a\\ b& 1-b\end{pmatrix}$ for some $a$ and $b$ in $(0,1)$. 
Fix $k\leqslant n$ and consider the set $\Gamma_k^n$ of paths $\gamma=(x_t)_{1\leqslant t\leqslant n}$ of length $n$ visiting $k$ times the state $1$. For every $x$ and $y$ in $\{0,1\}$ and every path $\gamma$ in $\Gamma_k^n$, let $N_\gamma(xy)$ denote the number of times $t\leqslant n$ such that $(x_t,x_{t+1})=(x,y)$. Finally, let $M(\gamma)=N_\gamma(01)$. 
Then, forgetting a possible discrepancy of $\pm1$ here and there, one sees that $N_\gamma(10)=M(\gamma)$, $N_\gamma(00)=n-k-M(\gamma)$ and $N_\gamma(11)=k-M(\gamma)$, hence, for every $\gamma$ in $\Gamma_k^n$,
$$
P(X=\gamma)=
(1-a)^{N_\gamma(00)}a^{N_\gamma(01)}b^{N_\gamma(10)}(1-b)^{N_\gamma(11)},
$$
is also, since $n$ and $k$ are fixed,
$$
P(X=\gamma)= (1-a)^{n-k-M(\gamma)}a^{M(\gamma)}b^{M(\gamma)}(1-b)^{k-M(\gamma)}\propto c^{M(\gamma)},
$$
where
$$
c=\frac{ab}{(1-a)(1-b)}.
$$
Thus, there exists some $A_k^n$ independent of $\gamma$ such that, for every $\gamma$ in $\Gamma_k^n$,
$$
P(X=\gamma\mid X\ \text{visits}\ 1\ \text{exactly}\ k\ \text{times})=A_k^n\cdot c^{M(\gamma)},
$$
and the question is to simulate a random path in $\Gamma_k^n$ following this distribution. 
Pseudo-algorithm: This assumes one is able to generate uniformly subsets $T$ of $\{1,2,\ldots,n\}$ of size $k$.


*

*Choose a path $\gamma$ uniformly in $\Gamma_k^n$, that is, choose uniformly a subset $T\subseteq\{1,2,\ldots,n\}$ of size $k$ and define $\gamma=(x_t)$ by $x_t=1$ if $t$ is in $T$ and $x_t=0$ otherwise.

*Compute $M(\gamma)$.

*Accept $\gamma$ with probability proportional to $c^{M(\gamma)}$.

*Repeat until a path $\gamma$ is accepted. Return $\gamma$.


To generate $T$, the following procedure might prove useful:


*

*Let $T=\varnothing$.

*Choose $x$ uniformly in $\{1,2,\ldots,n\}$. If $x$ is not already in $T$, add $x$ to $T$.

*Repeat until $T$ has size $k$. Return $T$.

A: Suppose you have transition probabilities $p_1 \leq p_2 \dots \leq p_k$, $\sum p_k = 1$. Define the distribution function in pseudo code:
function cdf(state) { a = math.random() % randomly generated number in U(0,1)
if 0 <= a < p_1, return 1
if p_1 <= a < p_1 + p_2, return 2
...
if p_1 + ... + p_{k-1} <= a < 1, return k

In the above, I have assumed that the transition probabilities are all the same, i.e. independent of the state. Otherwise, the $p_i$ will have to depend on state. 
The above code is just saying that if I have numbers $p_1$, ..., $p_n$ such that $\sum p_i = 1$, then I can partition the unit interval $[0,1]$ with these numbers. The area of the rectangle [$p_1 + ... + p_{k-1}, p_1 + ... + p_k$] with unit height is $p_k$, and that's the probability we want for the $k^{th}$ transition. This is the standard way of generating probability distributions, since if $X$ has cdf $F$, then $F(X) \sim U(0,1)$.
Now pick an initial state (say $s = 1$) and a state you want to measure, say $r = 4$. 
count = 0
for i = 1 : n
   s = cdf(s) // transition with the appropriate probabilities
   if s == 4, count++
   end
end

I should mention this algorithm is the worst. 
