Past exam of the MSc in advanced computer science course of the University of Oxford / 2021 / Bayesian statistical probabilistic programming / 1

Assume $0<Z={\sum }_{x}r(x)<\infty$, and work on the support $S=\{x:r(x)>0\}$. In rejection sampling, choose a normalized proposal probability mass function $g$ and a known finite constant $M$ satisfying $r(x)\le Mg(x)$ everywhere. Draw independent random variables $Y\sim g$ and $U\sim \mathrm{Unif}(0,1)$, accepting $Y$ when $U\le r(Y)/(Mg(Y))$. The joint probability of proposing and accepting $x$ is $r(x)/M$, so the overall acceptance probability is $Z/M$ and the accepted value has probability mass function $r(x)/Z$. Fresh trials therefore produce independent and identically distributed random variables with the desired probability distribution. For a fixed positive number $N$ of accepted values, the sample average $N^{-1}{\sum }_{i=1}^{N}f(Y_i)$ is an unbiased estimator of the expectation when that expectation is finite. Its variance is ${\mathrm{Var}}_p(f)/N$ if the second moment is finite. Finding a valid envelope can be difficult, and a small acceptance probability makes the method expensive.

For the Metropolis–Hastings algorithm, choose a proposal transition matrix $q$ and, at the current state $x$, propose $y\sim q(x,\cdot )$. Accept with the Metropolis–Hastings acceptance probabilityA proposal outside $S$ is rejected; an impossible proposal need not be assigned a ratio. Otherwise remain at $x$. This constructs a Markov chain satisfying detailed balance with $p$, without computing $Z$ or a global envelope. Estimate the expectation by averaging $f(X_i)$ over iterations, including repeated states after rejections. Keeping only accepted jumps generally changes the stationary distribution of the recorded process.

The Metropolis–Hastings algorithm produces dependent samples, so the variance of its average includes autocovariances; starting outside stationarity can also cause finite-run bias. With Markov chain irreducibility, suitable recurrence and an integrable $f$, the countable-state Markov chain ergodic theorem gives consistency. Aperiodicity is additionally needed for the usual convergence of one-time probability distributions. In contrast, accepted rejection sampling values need no initial transient. Neither method is universally faster: envelope rejection competes with the proposal quality and mixing time of the Markov chain. Nonnegativity of $f$ alone does not guarantee a finite expectation or variance.

Here $q(x,y)=p(y)$, with fresh proposal randomness independent of the past. For $x,y\in S$, the Metropolis–Hastings acceptance probability isThus $P(x,y)=p(y)$ is independent of $x$, and $X_1,\dots ,X_N$ are independent and identically distributed random variables with probability distribution $p$, even if $X_0$ was chosen arbitrarily in $S$. The method becomes ordinary Monte Carlo method sampling from $p$. The average of $f(X_1),\dots ,f(X_N)$ is an unbiased estimator, with the variance stated in part (a). Including an arbitrary initial state would spoil this finite-sample property. This is the independent stationary proposal lemma.

Neither Markov chain irreducibility nor fast convergence follows from detailed balance. On two states, the identity transition matrix satisfies detailed balance for the uniform probability distribution, but neither state can reach the other.

Even after imposing Markov chain irreducibility and aperiodicity, considerIt satisfies detailed balance for $\pi =(1/2,1/2)$. Its nontrivial eigenvalue is $1-2\epsilon$, and, from state zero, its total variation distance from $\pi$ after $n$ steps is $\frac{1}{2}(1-2\epsilon {)}^n$. This approaches zero arbitrarily slowly as $\epsilon$ decreases. Detailed balance without rapid convergence separates correctness of a stationary distribution from reachability and speed.

Put $a=p(x)q(x,y)$ and $b=p(y)q(y,x)$. For $a+b>0$, the off-diagonal probability flux of the proposed rule isIf $a=b=0$, define both acceptance probabilities to be zero; both fluxes are zero. The diagonal detailed balance equations hold automatically, including the probability of rejected proposals. Summing the detailed balance equations over $x$ gives stationarity:The alternative rule therefore preserves $p$. The unknown normalizing constant cancels after replacing $p$ by $r$. This Barker acceptance probability is no larger than the Metropolis–Hastings acceptance probability for the same proposed move: if $t=b/a>0$, then $t/(1+t)\le \min (1,t)$. Preservation of $p$ still requires separate Markov chain irreducibility and convergence checks before it guarantees useful sampling.

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 isAll 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 giveswhere the symbol $u_n$ denotes the upward transition probability. Thus the complete transition matrix iswith no downward entry at zero and all other entries zero. In matrix form,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.

For each adjacent pair, the detailed balance equations follows directly:For nonadjacent distinct states both fluxes vanish, and diagonal fluxes agree trivially. Since ${\sum }_n{\pi }_n=e^{-\lambda }{\sum }_n{\lambda }^n/n!=1$, this is a genuine stationary distribution, rather than merely an invariant measure of infinite mass.

Every upward and downward transition probability between neighbouring states is positive for $\lambda >0$. Any two states can therefore be joined by a finite path of positive probability: the Markov chain is irreducible. Moreover,State zero has period one; Markov chain irreducibility makes every state aperiodic. Finally, the standard countable-state recurrence theorem states that an irreducible Markov chain admitting a stationary distribution of total mass one is positive recurrent, with that stationary distribution unique. The countable-state Markov chain convergence theorem then states that Markov chain irreducibility, positive recurrence and aperiodicity imply $P^k(i,j)\to {\pi }_j$ for all states $i,j$. The hypotheses just verified hold on the infinite state space ${\mathbb{N}}_0$, not merely on a numerical truncation.^[ref]

Consequently, the unique stationary distribution and the limiting distribution are both $\mathrm{Poisson}(\lambda )$, from every initial state. The countable-state Markov chain ergodic theorem also gives $N^{-1}{\sum }_{k=1}^{N}f(X_k)\to {\sum }_n{\pi }_{n}f(n)$ almost surely whenever ${\sum }_n{\pi }_n|f(n)|<\infty$.^[ref]