= Solution
Use the <Poisson distribution> with <probability mass function> $\pi_n=e^{-\lambda}\lambda^n/n!$ on $\mathbb N_0$. A symmetric nearest-neighbour proposal with a reflecting boundary is
$$
q(n,n+1)=\tfrac12\quad(n\geq0),\qquad q(n,n-1)=\tfrac12\quad(n\geq1),\qquad q(0,0)=\tfrac12.
$$
All other entries vanish. In particular, the boundary proposal does not force a move to one: retaining the <probability> $q(0,0)=1/2$ makes $q(0,1)=q(1,0)$.
The adjacent ratios of the <Poisson distribution> are $\pi_{n+1}/\pi_n=\lambda/(n+1)$ and $\pi_{n-1}/\pi_n=n/\lambda$. The <Metropolis–Hastings acceptance probability> consequently gives
$$
u_n=\tfrac12\min\{1,\lambda/(n+1)\},\qquad d_n=\tfrac12\min\{1,n/\lambda\}\ (n\geq1),\qquad d_0=0,
$$
where the symbol $u_n$ denotes the upward <transition probability>. Thus the complete <transition matrix> is
$$
\boxed{P(n,n+1)=u_n,\quad P(n,n-1)=d_n,\quad P(n,n)=1-u_n-d_n,}
$$
with no downward entry at zero and all other entries zero. In matrix form,
$$
P=\begin{pmatrix}
1-u_0&u_0&0&0&\cdots\\
d_1&1-u_1-d_1&u_1&0&\cdots\\
0&d_2&1-u_2-d_2&u_2&\cdots\\
\vdots&\vdots&\ddots&\ddots&\ddots
\end{pmatrix}.
$$
Every row is nonnegative and sums to one. An executable <Python> step of this <Poisson Metropolis chain> is:
``
def poisson_step(n, lam, rng):
if lam <= 0 or n < 0 or int(n) != n:
raise ValueError("require lam > 0 and an integer n >= 0")
y = n + 1 if rng.random() < 0.5 else max(0, n - 1)
ratio = lam / (n + 1) if y > n else n / lam if y < n else 1.0
return y if rng.random() < min(1.0, ratio) else n
``
Supply a <pseudorandom number generator> such as `random.Random(2021)` and retain each returned state, including repeats. No <normalizing constant> or evaluation of a <factorial> is needed by the algorithm.
Back to article page
