MATH 565 — Fall 2026

Markov Chain Monte Carlo
Owen, Chapters 11–12

Fred J. Hickernell

August 24, 2026

Course Map

Course decks

• Introduction
• Generating Samples
• Markov Chain Monte Carlo
• Improving Efficiency
• Selected Topics

In this deck

• From Direct Sampling to Markov Chains
• Metropolis–Hastings Sampling
• Kernel Discrepancy
• Discrepancy as Integration Error
• Applications of Markov Chains
• Big Ideas

From Direct Sampling to Markov Chains

(Quasi-)Random
Number Generation

• When direct sampling is difficult
• The Markov property
• Examples of Markov chains

When direct sampling is difficult

(Quasi-)Random
Number Generation

• The quantile transform works directly for distributions with tractable quantiles
• Affine transformations generate correlated Gaussian vectors
• Acceptance–rejection sampling needs a proposal that bounds the target and may reject many proposals
• Markov chain Monte Carlo (MCMC) moves through the sample space using only an unnormalized target density

MCMC replaces independent draws by a dependent sequence whose distribution approaches the target

The Markov property

Probability

Let \(\mathcal X\) be the state space for \(X_0,X_1,\ldots\)

The sequence is a Markov chain when, for every measurable \(A\subseteq\mathcal X\),

\[\begin{align*} &\Prob(X_{i+1}\in A\mid X_i=x_i,\ldots,X_0=x_0)\\ &\qquad=\Prob(X_{i+1}\in A\mid X_i=x_i) \end{align*}\]

Where the chain goes next depends only on where it is now, not on how it arrived

We use Markov chains to

• generate samples whose distribution converges to a target
• model systems such as queues
• study long-run or asymptotic behavior

Examples of Markov chains

Probability

Assume \(Z_0,Z_1,\ldots\ \IIDsim\ \Norm(0,1)\)

Markov

\[\begin{align*} X_0&=0,& X_{i+1}&=X_i+Z_i\\ X_0&=0,& X_{i+1}&=i+Z_i \end{align*}\]

The second chain is not time homogeneous: its transition law changes with \(i\)

Not Markov as written

\[ X_0=X_1=0, \qquad X_{i+1}=X_{i-1}+Z_i \]

The pair \(Y_i=(X_{i-1},X_i)\) is Markov because its state retains the needed history

Metropolis–Hastings Sampling

(Quasi-)Random
Number Generation

• The Metropolis–Hastings algorithm
• Why the acceptance ratio works
• The Metropolis algorithm
• Metropolis practice
• A random-walk Metropolis example
• Tuning and diagnostics

The Metropolis–Hastings algorithm

(Quasi-)Random
Number Generation

Suppose

• the target has known, possibly unnormalized density \(\varrho\)
• \(\varrho_{\brandnew\mid\mathrm{old}}(z\mid x)\) proposes a new state \(z\) from the current state \(x\)

Given \(X_0\), repeat for \(i=0,\ldots,n-1\):

1. Generate \(Z\sim\varrho_{\brandnew\mid\mathrm{old}}(\cdot\mid X_i)\) and \(U\sim\Unif[0,1]\)

2. Compute

   \[ \alpha(X_i,Z)=\min\!\left\{1, \frac{\varrho(Z)\varrho_{\brandnew\mid\mathrm{old}}(X_i\mid Z)} {\varrho(X_i)\varrho_{\brandnew\mid\mathrm{old}}(Z\mid X_i)}\right\} \]

3. Set \(X_{i+1}=Z\) if \(U\le\alpha(X_i,Z)\); otherwise set \(X_{i+1}=X_i\)

Every iteration accepts either the proposal or the current state

Why the acceptance ratio works

Probability

The ratio compares the two directions of a proposed transition:

\[ \frac{\varrho(Z)\varrho_{\brandnew\mid\mathrm{old}}(X_i\mid Z)} {\varrho(X_i)\varrho_{\brandnew\mid\mathrm{old}}(Z\mid X_i)} = \frac{\text{target weight at $Z$}\times\text{return proposal}} {\text{target weight at $X_i$}\times\text{forward proposal}} \]

• Moves toward higher target density are favored
• Asymmetric proposals are corrected by the proposal-density ratio
• The normalizing constant in \(\varrho\) cancels

The transition kernel satisfies detailed balance with the target distribution, making the target invariant under suitable conditions

The Metropolis algorithm

(Quasi-)Random
Number Generation

If the proposal is symmetric,

\[ \varrho_{\brandnew\mid\mathrm{old}}(z\mid x) =\varrho_{\brandnew\mid\mathrm{old}}(x\mid z), \]

then the acceptance probability simplifies to

\[ \alpha(x,z)=\min\!\left\{1,\frac{\varrho(z)}{\varrho(x)}\right\} \]

For a random-walk proposal \(Z=X_i+\epsilon_i\) with symmetric \(\epsilon_i\):

• accept every move to a state with at least as much target density
• accept some downhill moves, which keeps the chain from becoming trapped immediately
• repeat \(X_i\) when a proposal is rejected

Metropolis practice

Practice

Use the uniform numbers in order and choose a random order of the candidate \(x\) values

\(i\)	0	1	2	3	4	5	6	7
\(U_i\)	0.6	0.1	0.3	0.5	0.9	0.2	0.1	0.4

\(x\)	\(-1.2\)	\(-0.8\)	\(-0.5\)	\(0.6\)	\(1.1\)	\(1.4\)	\(2.0\)	\(2.2\)
\(\varrho(x)\)	0.1	0.8	2.0	1.6	0.5	2.6	2.4	1.0

\(\exstar\) Starting at \(x=-1.2\), follow several Metropolis steps and record each accepted or repeated state

A random-walk Metropolis example

(Quasi-)Random
Number Generation

Target a bimodal Gaussian mixture using the symmetric proposal

\[ Z\mid X_i=x\sim\Norm(x,1.4^2) \]

Tuning and diagnostics

Statistics

MCMC is attractive because it

• samples from complicated, unnormalized densities
• is often straightforward to implement

The proposal and starting point require care:

• proposals that are too small explore slowly and produce highly dependent samples
• proposals that are too large have low acceptance rates and many repeats
• multimodal targets may trap a single chain in one mode

Useful responses include pilot runs, multiple chains, and parallel tempering

Diagnose both exploration of the full state space and concentration in high-density regions

Kernel Discrepancy

Low
Discrepancy

Discrepancy
Measures

• How can we measure sample quality?
• What discrepancy measures
• Symmetric positive definite kernels
• Discrepancy between distributions
• Empirical distributions
• Unbiased estimates of squared discrepancy

How can we measure sample quality?

Low
Discrepancy

Discrepancy
Measures

For samples \(\{x_i\}_{i=0}^{m-1}\) and \(\{z_j\}_{j=0}^{n-1}\), a kernel \(K\) defines

\[\begin{multline*} D^2\!\left(\{x_i\}_{i=0}^{m-1},\{z_j\}_{j=0}^{n-1};K\right) =\frac{1}{m^2}\sum_{i,j=0}^{m-1}K(x_i,x_j)\\ -\frac{2}{mn}\sum_{i=0}^{m-1}\sum_{j=0}^{n-1}K(x_i,z_j) +\frac{1}{n^2}\sum_{i,j=0}^{n-1}K(z_i,z_j) \end{multline*}\]

• Also called maximum mean discrepancy (MMD) in machine learning
• Direct computation requires \(\Order(\max(m,n)^2)\) kernel evaluations
• The answer depends on the chosen kernel and its parameters

What discrepancy measures

Low
Discrepancy

Discrepancy
Measures

Kernel discrepancy measures

• the difference between two probability distributions
• the worst-case integration error over the unit ball of a reproducing kernel Hilbert space
• the root mean squared integration error when the integrand is modeled as a Gaussian process with covariance kernel \(K\)

The kernel determines which differences between distributions matter most

Symmetric positive definite kernels

Statistics

A kernel \(K:\mathcal X\times\mathcal X\to\reals\) is symmetric positive definite when every Gram matrix

\[ \mK=(K(x_i,x_j))_{i,j=0}^{n-1} \]

is symmetric positive definite for distinct \(x_0,\ldots,x_{n-1}\)

Examples include

\[\begin{align*} K(t,x)&=\exp\!\left(-\norm{t-x}^2/h^2\right), &&\mathcal X=\reals^d,\\ K(t,x)&=\prod_{\ell=1}^d\left[1+\frac{\gamma_\ell^2}{2} \left(\abs{t_\ell-\tfrac12}+\abs{x_\ell-\tfrac12}-\abs{t_\ell-x_\ell}\right)\right], &&\mathcal X=[0,1]^d \end{align*}\]

The first is a squared exponential kernel; the second yields centered discrepancy

Discrepancy between distributions

Low
Discrepancy

Discrepancy
Measures

For independent \(X,X'\sim F_X\) and \(Z,Z'\sim F_Z\),

\[\begin{align*} D^2(F_X,F_Z;K) &=\Ex[K(X,X')]-2\Ex[K(X,Z)]+\Ex[K(Z,Z')]\\ &=\int_{\mathcal X^2}K(x,x')\,\dif(F_X-F_Z)(x)\,\dif(F_X-F_Z)(x')\\ &\ge0 \end{align*}\]

If both distributions have densities,

\[ D^2(F_X,F_Z;K) =\int_{\mathcal X^2}K(x,x') \{\varrho_X(x)-\varrho_Z(x)\} \{\varrho_X(x')-\varrho_Z(x')\}\,\dif x\,\dif x' \]

Empirical distributions

Low
Discrepancy

Discrepancy
Measures

For the empirical distribution of \(z_0,\ldots,z_{n-1}\),

\[\begin{align*} D^2(F_X,\{z_i\}_{i=0}^{n-1};K) &=\int_{\mathcal X^2}K(x,x')\,\dif F_X(x)\,\dif F_X(x')\\ &\quad-\frac{2}{n}\sum_{i=0}^{n-1}\int_{\mathcal X}K(x,z_i)\,\dif F_X(x)\\ &\quad+\frac{1}{n^2}\sum_{i,j=0}^{n-1}K(z_i,z_j) \end{align*}\]

The same formula covers one target distribution and one sample or two samples

Unbiased estimates of squared discrepancy

Statistics

The squared discrepancy is nonnegative, but an unbiased estimator may be negative

Replace within-sample sums by off-diagonal averages:

\[\begin{multline*} D_{\mathrm{unb}}^2\!\left(\{x_i\}_{i=0}^{m-1},\{z_j\}_{j=0}^{n-1};K\right) =\frac{1}{m(m-1)}\sum_{\substack{i,j=0\\i\ne j}}^{m-1}K(x_i,x_j)\\ -\frac{2}{mn}\sum_{i=0}^{m-1}\sum_{j=0}^{n-1}K(x_i,z_j) +\frac{1}{n(n-1)}\sum_{\substack{i,j=0\\i\ne j}}^{n-1}K(z_i,z_j) \end{multline*}\]

Unbiased describes the estimator’s expectation, not the sign of each realization

Discrepancy as Integration Error

Estimation

• A reproducing kernel Hilbert space
• Worst-case integration error
• The discrepancy identity
• Average-case integration error
• What the kernel controls
• Centered discrepancy and smoothness

A reproducing kernel Hilbert space

Statistics

Every symmetric positive definite kernel \(K\) defines a Hilbert space \(\mathcal H_K\) of functions on \(\mathcal X\) with the reproducing property

\[ f(x)=\langle K(\cdot,x),f\rangle_{\mathcal H_K} \]

Assume \(\int_{\mathcal X^2}K(x,x')\,\dif F_X(x)\,\dif F_X(x')\) is finite

Then both

\[\begin{align*} f&\longmapsto\int_{\mathcal X}f(x)\,\dif F_X(x),\\ f&\longmapsto\int_{\mathcal X}f(x)\,\dif F_X(x)-\frac1n\sum_{i=0}^{n-1}f(z_i) \end{align*}\]

are bounded linear functionals on \(\mathcal H_K\)

Worst-case integration error

Estimation

By the Riesz representation theorem, some \(\zeta\in\mathcal H_K\) satisfies

\[ \int_{\mathcal X}f(x)\,\dif F_X(x)-\frac1n\sum_{i=0}^{n-1}f(z_i) =\langle\zeta,f\rangle_{\mathcal H_K} \]

Therefore, by Cauchy–Schwarz,

\[ \left\lvert\int_{\mathcal X}f(x)\,\dif F_X(x)-\frac1n\sum_{i=0}^{n-1}f(z_i)\right\rvert \le\norm{\zeta}_{\mathcal H_K}\norm{f}_{\mathcal H_K} \]

This separates sample quality, \(\norm{\zeta}_{\mathcal H_K}\), from integrand difficulty, \(\norm{f}_{\mathcal H_K}\)

The discrepancy identity

Low
Discrepancy

Discrepancy
Measures

The representer is

\[ \zeta(x')=\int_{\mathcal X}K(x,x')\,\dif F_X(x) -\frac1n\sum_{i=0}^{n-1}K(z_i,x') \]

Its squared norm is exactly the squared discrepancy:

\[\begin{align*} \norm{\zeta}_{\mathcal H_K}^2 &=\int_{\mathcal X^2}K(x,x')\,\dif F_X(x)\,\dif F_X(x')\\ &\quad-\frac2n\sum_{i=0}^{n-1}\int_{\mathcal X}K(x,z_i)\,\dif F_X(x) +\frac1{n^2}\sum_{i,j=0}^{n-1}K(z_i,z_j)\\ &=D^2(F_X,\{z_i\}_{i=0}^{n-1};K) \end{align*}\]

Thus discrepancy is the worst-case error over \(\norm{f}_{\mathcal H_K}\le1\)

Average-case integration error

Statistics

Suppose \(f\) is modeled as a zero-mean Gaussian process with covariance kernel \(K\)

For deterministic \(z_0,\ldots,z_{n-1}\),

\[\begin{align*} &\Ex_f\!\left[ \left\{\int_{\mathcal X}f(x)\,\dif F_X(x)-\frac1n\sum_{i=0}^{n-1}f(z_i)\right\}^2 \right]\\ &\quad=\int_{\mathcal X^2}K(x,x')\,\dif F_X(x)\,\dif F_X(x') -\frac2n\sum_{i=0}^{n-1}\int_{\mathcal X}K(x,z_i)\,\dif F_X(x) +\frac1{n^2}\sum_{i,j=0}^{n-1}K(z_i,z_j)\\ &\quad=D^2(F_X,\{z_i\}_{i=0}^{n-1};K) \end{align*}\]

Discrepancy is also the root mean squared error under the Gaussian-process model

What the kernel controls

Statistics

• The numerical value of discrepancy depends on \(K\) and its hyperparameters
• Multiplying \(K\) by \(c^2\) multiplies discrepancy by \(\abs{c}\)
• Kernel structure encodes assumptions about
  • domain and smoothness
  • periodicity
  • coordinate importance
• In the average-case setting, hyperparameters may be estimated by empirical Bayes or maximum likelihood

For IID \(Z_i\sim F_X\),

\[ \Ex\!\left[D^2(F_X,\{Z_i\}_{i=0}^{n-1};K)\right] =\frac1n\left\{ \int K(x,x)\,\dif F_X(x)-\int K(x,x')\,\dif F_X(x)\,\dif F_X(x') \right\} \]

Centered discrepancy and smoothness

Low
Discrepancy

Discrepancy
Measures

For \(\mathcal X=[0,1]^d\), the centered-discrepancy kernel is

\[ K(t,x)=\prod_{\ell=1}^d\left[1+\frac{\gamma_\ell^2}{2} \left(\abs{t_\ell-\tfrac12}+\abs{x_\ell-\tfrac12}-\abs{t_\ell-x_\ell}\right)\right] \]

Its norm combines the function value at the center with weighted square-integrals of mixed derivatives:

\[\begin{multline*} \norm{f}_{\mathcal H_K}^2 =\abs{f(1/2,\ldots,1/2)}^2 +\norm{\frac{\partial f(x_1,1/2,\ldots)}{\gamma_1\partial x_1}}_2^2\\ +\norm{\frac{\partial^2 f(x_1,x_2,1/2,\ldots)} {\gamma_1\gamma_2\partial x_1\partial x_2}}_2^2 +\cdots +\norm{\frac{\partial^d f(x)} {\gamma_1\cdots\gamma_d\partial x_1\cdots\partial x_d}}_2^2 \end{multline*}\]

The weights \(\gamma_\ell\) express coordinate importance

Applications of Markov Chains

Probability

Estimation

• Maximum likelihood estimation
• Bayesian inference
• Queueing systems
• Event-driven queue simulation
• Tracking customer time

Maximum likelihood estimation

Estimation

For parameter \(\theta\) and observed data \(y\), the likelihood is

\[ L(\theta;y)=\varrho_{Y\mid\theta}(y\mid\theta) \]

The maximum likelihood estimator maximizes \(L(\theta;y)\)

For \(Y_0,\ldots,Y_{n-1}\ \IIDsim\ \Norm(\mu,\sigma^2)\),

\[\begin{align*} L(\mu,\sigma^2;y) &=(2\mpi\sigma^2)^{-n/2} \exp\!\left[-\frac1{2\sigma^2}\sum_{i=0}^{n-1}(y_i-\mu)^2\right],\\ \hmu_{\mathrm{MLE}}&=\frac1n\sum_{i=0}^{n-1}y_i,& \hsigma^2_{\mathrm{MLE}}&=\frac1n\sum_{i=0}^{n-1}(y_i-\hmu_{\mathrm{MLE}})^2 \end{align*}\]

Bayesian inference

Probability

Estimation

Combine a prior density with the likelihood:

\[\begin{align*} \varrho_{\mathrm{post}}(\theta\mid y) &=\frac{L(\theta;y)\varrho_{\mathrm{prior}}(\theta)}{\varrho_Y(y)}\\ &\propto L(\theta;y)\varrho_{\mathrm{prior}}(\theta) \end{align*}\]

• The evidence \(\varrho_Y(y)\) depends only on the observed data
• MCMC needs only the unnormalized posterior
• Posterior samples estimate means, quantiles, credible intervals, and other summaries
• Optimization of the unnormalized posterior yields a maximum a posteriori estimate

Bayesian inference is a primary reason to sample from complicated unnormalized densities

Queueing systems

Probability

Consider a single-server queue observed at event times

\[\begin{align*} A_i&=\text{time until the next arrival},& A_i&\sim F_A,\\ N_i&=\text{number of customers in the system},\\ S_i&=\text{remaining service time},& S_i&\sim F_S \end{align*}\]

The state \((A_i,N_i,S_i)\) is Markov when arrival and service times are generated independently from their specified laws

At each event, either

• the next customer arrives before service finishes, or
• the current service finishes before the next arrival

Event-driven queue simulation

Practice

Initialize A[0] ~ F_A, N[0] = 0, S[0] ~ F_S
while the stopping rule is not met:
    if N[i] == 0:
        advance to the next arrival
        N[i+1] = 1
        generate A[i+1] ~ F_A
    elif A[i] < S[i]:
        advance to the next arrival
        N[i+1] = N[i] + 1
        S[i+1] = S[i] - A[i]
        generate A[i+1] ~ F_A
    else:
        advance to the next service completion
        N[i+1] = N[i] - 1
        A[i+1] = A[i] - S[i]
        generate S[i+1] ~ F_S

Events replace fixed time steps, avoiding work when nothing changes

Tracking customer time

Practice

Add

\[\begin{align*} T_i&=\text{total elapsed time},\\ W_j&=\text{waiting plus service time for customer $j$} \end{align*}\]

At an arrival after elapsed time \(\Delta t\):

• add \(\Delta t\) to every customer currently in the system
• initialize the arriving customer’s clock
• update the queue length and residual event times

At a service completion:

• add the completed service interval to all waiting customers
• record the departing customer’s \(W_j\)
• begin service for the next customer

\(\exstar\) Which state variables are needed if service times are not memoryless?

Big Ideas

Probability

Statistics

Estimation

• A Markov chain’s next-state law depends on its current state
• Metropolis–Hastings corrects a proposal mechanism so that the target distribution is invariant
• MCMC trades independent samples for access to complicated unnormalized densities
• Kernel discrepancy compares distributions and measures both worst-case and average-case integration error
• The kernel specifies which distributional differences and integrand features matter
• Markov chains model computational algorithms as well as dynamic systems such as queues