Introduction

This short method note presents a systematic framework for constructing Markov chains with transition matrices (or kernels) that achieve specific stationary distributions, including the Gaussian, heavy-tailed Cauchy, extreme-value Gumbel, and Bose–Einstein (Planck’s law) distributions.

Using the Metropolis–Hastings algorithm, we derive detailed balance conditions for each distribution, ensuring convergence to the target distributions. For the Gaussian distribution, a symmetric proposal (e.g., Gaussian) is used, while the heavy-tailed Cauchy distribution employs a symmetric proposal to capture its undefined-mean behavior. The Gumbel distribution utilizes an asymmetric proposal tailored to its tail behavior, and Planck’s law is derived via a birth-death process mimicking photon occupation numbers. The derivations are validated through rigorous checks of detailed balance, ergodicity (irreducibility and aperiodicity), and continuous limits (e.g., Fokker–Planck equations for Gaussian).

Applications span statistics (robust modeling with Cauchy), risk assessment (Gumbel for extreme events), and physics (Planck’s law for blackbody radiation).

Keywords: Markov chain Monte Carlo (MCMC), Metropolis–Hastings algorithm, detailed balance, Gaussian distribution, Cauchy distribution, Gumbel distribution, Planck’s law, extreme value theory, robust statistic

a. Deriving the Gaussian Distribution

Step 1: Define the Discrete State Space

Define the state space as discrete positions:

$x_i = i \cdot \Delta \quad \text{for } i \in \mathbb{Z},$

where $\Delta$ is a small step size. In the limit $\Delta \to 0$, the chain becomes “continuous.”

Step 2: Specify the Target Stationary Distribution

We want the chain’s equilibrium probability (or invariant measure) on state $x_i$ to be:

$\pi(x_i) \propto \exp\left(-\frac{2}{\sigma^2} x_i^2\right).$

Step 3: Impose Detailed Balance

The Markov chain will have transition probabilities $P(x_i \to x_j)$ given by the elements of the matrix $B$. Detailed balance requires that for any pair of neighboring states:

$\pi(x_i) P(x_i \to x_j) = \pi(x_j) P(x_j \to x_i)$

For nearest-neighbor transitions, from a state $x_i$, one may only jump to $x_{i+1}$ or $x_{i-1}$.

Step 4: Choose a Symmetric Proposal and Determine the Acceptance Factor

Start with a symmetric proposal probability $Q(x_i, x_j)$ (for example, $Q(x_i, x_j) = \frac{1}{2}$ if $x_j$ is a neighbor of $x_i$). Apply an acceptance rule $A(x_i, x_j)$ so that:

$P(x_i \to x_j) = Q(x_i, x_j) A(x_i, x_j).$

To satisfy detailed balance:

$A(x_j \to x_i) A(x_i \to x_j) = \frac{\pi(x_i)}{\pi(x_j)}.$

A standard choice (as in the Metropolis algorithm) is:

$A(x_i \to x_j) = \min\left{1, \frac{\pi(x_i)}{\pi(x_j)}\right}.$

Step 5: Recover the Continuous Limit

In the limit where $\Delta \to 0$, the resulting Markov chain dynamics can be described by a Fokker–Planck equation:

$\frac{\partial}{\partial t} p(x, t) = -\frac{\partial}{\partial x} [D_1(x) p(x, t)] + \frac{\partial^2}{\partial x^2} [D_2(x) p(x, t)].$

For a process obeying:

$p_{\text{stat}}(x) \propto \exp\left[\int \frac{D_2}{D_1}(x) \, dx\right].$

If one chooses $D_1 \propto -x$ and $D_2$ constant, then:

$p_{\text{stat}}(x) \propto \exp\left(-\frac{2}{\sigma^2} x^2\right),$

which is exactly a Gaussian distribution.

Step 6: Summarize the Role of Matrix B

The matrix $B$ representing the transition probabilities is designed such that:

1. It allows moves only between “neighboring” states.
2. The acceptance probabilities satisfy detailed balance with respect to

$\pi(x_i) \propto \exp\left(-\frac{2}{\sigma^2} x_i^2\right)$

1. The chain is ergodic (irreducible and aperiodic), ensuring convergence to the stationary Gaussian.

---

b. Deriving Planck’s Law via a Markov Chain Model

Step 1: Identify the Quantized Energy States

For a particular photon mode of frequency $\nu$, the energy is quantized. Denote the state by an occupation number $n = 0, 1, 2, \ldots$, with energy:

$E_n = n h \nu.$

Step 2: Specify the Target Stationary Distribution

We wish the stationary probability for state $n$ to follow the Boltzmann factor:

$\pi(n) \propto \exp(-\beta n h \nu), \quad \beta = \frac{1}{kT}.$

This can be normalized to yield a geometric distribution:

$\pi(n) = (1 - q) q^n, \quad q = \exp(-\beta h \nu).$

Step 3: Construct the Transition Matrix B

Consider transitions between states $n$ and $n+1$:

• Allow “abs
• Allow “absorption” (n → n + 1) and “emission” (n + 1 → n) processes.
• Assume symmetric trial steps with probabilities governed by an underlying proposal ( Q(n, n+1) = Q(n+1, n) ).

To enforce detailed balance:

$\pi(n) P(n \to n+1) = \pi(n+1) P(n+1 \to n).$

If we set:

$P(n \to n+1) = Q(n, n+1) A(n, n+1), \quad P(n+1 \to n) = Q(n+1, n) A(n+1, n),$

and taking $Q$ symmetric, the acceptance probabilities must satisfy:

$A(n+1, n) A(n, n+1) = \exp(-\beta h \nu).$

A common Metropolis choice is:

$A(n \to n+1) = \min\left{1, \exp(-\beta h \nu)\right}, \quad A(n+1 \to n) = \min\left{1, \exp(\beta h \nu)\right}.$

Given that $\exp(\beta h \nu) > 1$, one would typically have:

$A(n+1 \to n) = 1 \quad \text{and} \quad A(n \to n+1) = \exp(-\beta h \nu).$

Step 4: Obtain the Stationary Distribution

Under the imposed detailed balance condition, the stationary distribution must be:

$\pi(n) \propto \exp(-\beta n h \nu).$

Normalization leads to:

$\pi(n) = (1 - \exp(-\beta h \nu)) \exp(-\beta h \nu \cdot n).$

Step 5: Connecting to Planck’s Law

The average energy of a single mode is:

$\langle E \rangle = \sum_{n=0}^{\infty} (n h \nu) \pi(n).$

Let’s compute this sum:

$\langle E \rangle = h \nu \cdot \sum_{n=0}^{\infty} n \cdot (1 - q) q^n,$

where $q = \exp(-\beta h \nu)$

We know that:

$\sum_{n=0}^{\infty} n q^n = \frac{q}{(1 - q)^2}.$

So:

$\langle E \rangle = h \nu \cdot (1 - q) \cdot \frac{q}{(1 - q)^2} = h \nu \cdot \frac{q}{1 - q}.$

But $q = \exp(-k T h \nu)$ , so:

$\frac{1 - q}{q} = \frac{1 - \exp(-k T h \nu)}{\exp(-k T h \nu)} = \frac{\exp(k T h \nu) - 1}{1}.$

Therefore:

$\langle E \rangle = \frac{h \nu}{\exp(k T h \nu) - 1}.$

Now, incorporating the density of states factor (number of modes per unit frequency interval). For a cavity radiator, the density of states $g(\nu)$ in three dimensions is:

$g(\nu) = \frac{8 \pi \nu^2}{c^3}.$

Multiplying this by the average energy per mode gives us Planck’s spectral energy density:

$u(\nu) = g(\nu) \cdot \langle E \rangle = \frac{8 \pi \nu^2}{c^3} \cdot \frac{h \nu}{\exp(\frac{h \nu}{k T}) - 1}.$

Which simplifies to:

$u(\nu) = \frac{8 \pi h \nu^3}{c^3} \cdot \frac{1}{\exp(\frac{h \nu}{k T}) - 1}.$

Thus, through a Markov chain that satisfies detailed balance and produces the correct geometric (or Bose–Einstein) stationary distribution for the occupation numbers, the equilibrium energy distribution for photons—Planck’s law—can be recovered.

c. Deriving The Heavy-tailed Cauchy Distribution

Step 1: Define the State Space and Target Distribution

• State Space: Continuous state space ( \mathbb{R} ).
• Target Distribution: [ \pi(x) \propto \frac{1}{1 + x^2}. ]

Step 2: Construct the Transition Kernel

• Symmetric Proposal: Use a symmetric proposal distribution ( Q(x, y) ), such as a Gaussian: [ Q(x, y) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left( -\frac{(y - x)^2}{2\sigma^2} \right). ]
• Acceptance Probability: [ A(x, y) = \min\left{ 1, \frac{\pi(y)}{\pi(x)} \right} = \min\left{ 1, \frac{1 + x^2}{1 + y^2} \right}. ]
• Transition Kernel: [ P(x, dy) = Q(x, y) A(x, y) dy. ]

Step 3: Verify Detailed Balance

• Detailed Balance Condition: [ \pi(x) P(x, dy) = \pi(y) P(y, dx). ]
• Substitution: [ \pi(x) Q(x, y) A(x, y) = \pi(y) Q(y, x) A(y, x). ]
• Symmetry of ( Q ): ( Q(x, y) = Q(y, x) ).
• Acceptance Ratio: [ A(x, y) = \min\left{ 1, \frac{\pi(y)}{\pi(x)} \right} = \min\left{ 1, \frac{1 + x^2}{1 + y^2} \right}. ] This ensures detailed balance.

Step 4: Ensure Ergodicity

• Irreducibility: Gaussian proposal allows moves across ( \mathbb{R} ), ensuring irreducibility.
• Aperiodicity: Gaussian proposal does not create cycles, ensuring aperiodicity.

Conclusion for Cauchy:

The transition kernel ( P(x, dy) ) constructed using the Metropolis-Hastings approach with a symmetric Gaussian proposal satisfies detailed balance with the Cauchy distribution as its stationary distribution. The chain is ergodic, ensuring convergence to the Cauchy distribution, modeling systems with large fluctuations and undefined means. It is useful in robust statistics, physics, and finance.

d. Deriving Extreme-Value Gumbel Distribution

Step 1. State Space and Target Distribution

State Space: Continuous ( \mathbb{R} ).

Target Distribution:
[ \pi(x) \propto \exp\left( -x - e^{-x} \right). ]

Step 2. Transition Kernel Construction

Asymmetric Proposal:
Define ( Q(x, y) ) with distinct normalization constants for upward and downward moves: [ Q(x, y) = \begin{cases} \frac{1}{\sigma_1} \exp\left( -\frac{y - x}{\sigma_1} \right), & \text{if } y > x, \ \frac{1}{\sigma_2} \exp\left( \frac{y - x}{\sigma_2} \right), & \text{if } y \leq x, \end{cases} ] where ( \sigma_1 \neq \sigma_2 ) ensures asymmetry.

Step 3. Acceptance Probability

Detailed Balance Condition:
Compute the proposal ratio ( \frac{Q(y, x)}{Q(x, y)} ): - For ( y > x ):
[ Q(y, x) = \frac{1}{\sigma_2} \exp\left( -\frac{x - y}{\sigma_2} \right) = \frac{1}{\sigma_2} \exp\left( \frac{y - x}{\sigma_2} \right), ] [ \frac{Q(y, x)}{Q(x, y)} = \frac{\frac{1}{\sigma_2} \exp\left( \frac{y - x}{\sigma_2} \right)}{\frac{1}{\sigma_1} \exp\left( -\frac{y - x}{\sigma_1} \right)} = \frac{\sigma_1}{\sigma_2} \exp\left( (y - x)\left( \frac{1}{\sigma_2} + \frac{1}{\sigma_1} \right) \right). ]

For ( y \leq x ):
[ Q(y, x) = \frac{1}{\sigma_1} \exp\left( -\frac{x - y}{\sigma_1} \right) = \frac{1}{\sigma_1} \exp\left( \frac{y - x}{\sigma_1} \right), ] [ \frac{Q(y, x)}{Q(x, y)} = \frac{\frac{1}{\sigma_1} \exp\left( \frac{y - x}{\sigma_1} \right)}{\frac{1}{\sigma_2} \exp\left( \frac{y - x}{\sigma_2} \right)} = \frac{\sigma_2}{\sigma_1} \exp\left( (y - x)\left( \frac{1}{\sigma_1} - \frac{1}{\sigma_2} \right) \right). ]

Acceptance Probability:
[ A(x, y) = \min\left{ 1, \frac{\pi(y) Q(y, x)}{\pi(x) Q(x, y)} \right} = \min\left{ 1, \frac{\pi(y)}{\pi(x)} \cdot \frac{Q(y, x)}{Q(x, y)} \right}. ] Substituting the proposal ratios:

For ( y > x ):
[ A(x, y) = \min\left{ 1, \frac{\pi(y)}{\pi(x)} \cdot \frac{\sigma_1}{\sigma_2} \exp\left( (y - x)\left( \frac{1}{\sigma_2} + \frac{1}{\sigma_1} \right) \right) \right}. ]

For ( y \leq x ):
[ A(x, y) = \min\left{ 1, \frac{\pi(y)}{\pi(x)} \cdot \frac{\sigma_2}{\sigma_1} \exp\left( (y - x)\left( \frac{1}{\sigma_1} - \frac{1}{\sigma_2} \right) \right) \right}. ]

Step 4. Balance Verification

The acceptance probability ensures detailed balance:
[ \pi(x) Q(x, y) A(x, y) = \pi(y) Q(y, x) A(y, x). ] This holds because the proposal ratio ( \frac{Q(y, x)}{Q(x, y)} ) explicitly accounts for the asymmetry in ( Q(x, y) ).

Step 5. Ergodicity

Irreducibility: The proposal allows transitions to any ( y \in \mathbb{R} ) from any ( x \in \mathbb{R} ), ensuring the chain can reach all states. Aperiodicity: The continuous proposal density prevents periodic behavior (e.g., fixed cycles).

Conclusion for Gumbel:

The transition kernel ( P(x, dy) ) constructed using the Metropolis-Hastings approach with an asymmetric exponential proposal satisfies detailed balance with the Gumbel distribution as its stationary distribution. The chain is ergodic, ensuring convergence to the Gumbel distribution, modeling systems dominated by extreme values. It is widely used in risk assessment, climatology, and reliability engineering to analyze phenomena like extreme weather events, financial market crashes, and system failure thresholds.

Notes

• Cauchy Distribution:
• Transition Kernel: Metropolis-Hastings with symmetric Gaussian proposal.
• Verification: Detailed balance confirmed analytically.
• Ergodicity: Chain is irreducible and aperiodic.

• Conclusion: Methodology sound and rigorously derived.

• Gumbel Distribution:

• Transition Kernel: Metropolis-Hastings with asymmetric exponential proposal.
• Verification: Detailed balance maintained through careful design of acceptance.
• Ergodicity: Chain ensures irreducibility and aperiodicity.
• Conclusion: Construction is valid with thorough analytical backing.

Cauchy Distribution: A Metropolis-Hastings algorithm with a symmetric Cauchy proposal preserves the heavy-tailed target distribution. Detailed balance is ensured by Metropolis acceptance rules, and ergodicity is maintained by proposal symmetry.

Gumbel Distribution: An asymmetric proposal tailored to the distribution’s tail behavior, adjusted via Metropolis-Hastings acceptance ratios, captures the skewed equilibrium. Ergodicity is ensured by non-symmetric jumps and acceptance probability corrections.

Convergence to Target Distributions Both transition matrices are designed to satisfy detailed balance and ensure ergodicity, validating their correctness in generating samples from the Cauchy and Gumbel distributions, respectively.

Cauchy Distribution

Transition Matrix / Markov Chain:
The constructed Markov chain uses a Metropolis-Hastings algorithm with a symmetric proposal (e.g., Gaussian). The transition kernel ensures detailed balance with respect to the Cauchy distribution:
[ \pi(x) \propto \frac{1}{1 + x^2}. ]
The chain’s stationary distribution is the Cauchy distribution itself.

Interpretation:
- Heavy-Tailed Behavior: The chain models systems with heavy-tailed fluctuations, such as financial returns or physical systems with long-range interactions.
- Robustness to Outliers: The Cauchy distribution’s lack of a defined mean makes it useful for modeling phenomena where extreme values are common.
- Equilibrium: The chain converges to a state where transitions balance the heavy-tailed nature of the distribution, ensuring the system’s long-term behavior matches the Cauchy distribution.

Gumbel Distribution

Transition Matrix / Markov Chain:
The constructed Markov chain uses a Metropolis-Hastings algorithm with an asymmetric proposal tailored to the Gumbel distribution’s tail behavior. The proposal density is:
[ Q(x, y) = \begin{cases} \frac{1}{\sigma_{\text{up}}} \exp\left( -\frac{y - x}{\sigma_{\text{up}}} \right), & y > x, \ \frac{1}{\sigma_{\text{down}}} \exp\left( \frac{y - x}{\sigma_{\text{down}}} \right), & y < x. \end{cases} ]
The acceptance probability ensures detailed balance:
[ A(x, y) = \min\left{1, \frac{\pi(y) Q(y, x)}{\pi(x) Q(x, y)} \right}. ]

Interpretation:

Extreme Value Theory: The chain models systems where the stationary distribution represents maxima or minima of processes, such as annual flood levels, stock market crashes, or extreme weather events.

Asymmetric Dynamics: The proposal’s asymmetry captures the Gumbel distribution’s skew, favoring transitions toward higher values (for the standard Gumbel) or lower values (for the reflected Gumbel).

Equilibrium: The chain converges to a state where transitions balance the accumulation of extreme values, ensuring the system’s long-term behavior matches the Gumbel distribution.