Vtyjkyuk

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')

First a reformulation of the question, then a pseudo-algorithm to solve it.

Reformulation: Let $ngeqslant1$ and assume that $X=(X_t)_1leqslant tleqslant n$ is a Markov chain of length $n$ on the state space $0,1$ with transition probability matrix $beginpmatrix1-a & a\ b& 1-bendpmatrix$ for some $a$ and $b$ in $(0,1)$.

Fix $kleqslant n$ and consider the set $Gamma_k^n$ of paths $gamma=(x_t)_1leqslant tleqslant 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 $tleqslant 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=fracab(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=gammamid X textvisits 1 textexactly k texttimes)=A_k^ncdot 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 $Tsubseteq1,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$.

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.