where In order to compute the entropy rate of HMPs based on (141), the key is to estimate the probability distribution of L$_{n}$ , P$_{L}$$_{n}$ . In [22] for binary Markov input and the binary symmetric channel, L$_{n}$ is considered as a 1-dim M -state Markov chain by quantizing the dynamic system expressed in (134). Then the distribution of L$_{∞}$ is the stationary distribution of the quantized Markov process and can be computed easily through eigenvector solving method. In this paper because the number of states of the Markov input, m is larger than 2 and the HMPs is observed through the AWGN channel, directly quantizing the dynamic system (134) will generate a quantized Markov chain with M m − 1 states, which is very difficult to deal with when large M is selected for good estimation precision. According to (134), since L$_{n}$$_{+1}$ is only dependent on L$_{n}$ and Y$_{n}$$_{+1}$ , { L$_{n}$ } can be considered as a Markov process. In order to compute the stationary probability distribution P$_{L}$$_{∞}$ , we can evolve the distribution of L$_{n}$ based on (134) from any initial distribution P$_{L}$$_{0}$ . When n is large enough, the distribution P$_{L}$$_{n}$ converges to P$_{L}$$_{∞}$ . A Monte Carlo algorithm for approximating h ( Y ) is introduced as follows: When M is very large, histogram method can be used to describe { L$_{N,k}$ } and reduce the computational load. 1) Idealized duty cycle constraint ( q, 0) : One implication of Theorem 1 is that directly computing the capacity-achieving input distribution requires solving an optimization problem with infinite variables 28 B. Numerical Results 1) Initialize M particles { L$_{0}$$_{,}$$_{1}$ , · · · , L$_{0}$$_{,M}$ } , L$_{0}$$_{,k}$ can be simply sampled from the ( m − 1) -dim Uniform distribution with each dimension on [ − max( X ( i $^{)}$) , max( X ( i $^{)}$)] . 2) for n = 0 , 1 , 2 , · · · , N , iteratively evolve the particles { L$_{0}$$_{,}$$_{1}$ , · · · , L$_{0}$$_{,M}$ } based on (134), where each y$_{n}$$_{+1}$$_{,k}$ is sampled according to r ( y, L$_{n,k}$ ) . 3) when N is large enough, { L$_{N,k}$ } can be used to estimate h ( Y ) as h ( Y ) ≈ − 1 M M ∑ k =1 ∫ r ( y, L$_{N,k}$ ) log r ( y, L$_{N,k}$ ) dy. (143) r ( y, L$_{n}$ ) = m − 1 ∑ i =0 φ ( y − x ( i $^{)}$)$^{m}$$^{−}$$^{1}$ ∑ k =0 e L ( i ) n ∑ m − 1 i =0 e L ( i ) n P$_{X}$$_{2}$$_{|}$$_{X}$$_{1}$ ( x ( i $^{)}$| x ( k $^{)}$) . (142)