Create a Markov chain that has a stationary distribution equal to the Zipf distribution using MCMC-MH.

Question

Background

I'm trying to learn about MCMC-MH by implementing one of the examples in a textbook, specifically problem 12.1.3. A description of how to solve the problem is given, but not an implementation.

The problem

Let $M > 2$ be an integer. An random variable $X$ has the zipf distribution with parameter $a > 0$ if its PMF is :

$$ P(X=k) = \frac{\frac{1}{k^a}}{\sum_{j=1}^{M} \frac{1}{j^a}} $$

for $k = 1, 2, ...$ (and $0$ other wise). Create a Markov chain $X_0, X_1, ...$ whose stationary distribution is the Zipf distribution, and such that $|X_{n+1} - X_n| \leq 1$ for all $n$.

The solution

We can use the Metropolis-Hastings algorithm, after coming up with a proposal distribution. The choice of proposal distribution is a random walk on $1, 2, ..., M$. From state $i$ with $i\neq 1$ or $i\neq M$, we move to state $i+1$ or $i-1$ with probability $0.5$ each. If we are at state $1$, we stay there or move to state $2$ each with $0.5$ probability. If we are in state M, we stay there or move to state $M-1$ each with probability $0.5$.

Let $P$ be the transition matrix for this chain and $X_0$ be the starting state. We can generate a chain $X_0, X_1, ...$ as follows, given that we are in state $i$:

1. Generate a proposal state, $j$, according to the proposal chain P.

2. Accept the proposal with probability $min(\frac{i^a}{j^a}, 1)$. If the proposal is accepted, we go to state $j$ and stay at $i$ otherwise.

The implementation (Python)

import numpy as np
from scipy.stats import zipf
import matplotlib.pyplot as plt
from matplotlib.pyplot import stem
import seaborn
from collections import Counter


def transition_prob():
    """
    :return: (np.ndarry). For example, when M=10:

        [[0.5 0.5 0.  0.  0.  0.  0.  0.  0.  0. ]
         [0.5 0.  0.5 0.  0.  0.  0.  0.  0.  0. ]
         [0.  0.5 0.  0.5 0.  0.  0.  0.  0.  0. ]
         [0.  0.  0.5 0.  0.5 0.  0.  0.  0.  0. ]
         [0.  0.  0.  0.5 0.  0.5 0.  0.  0.  0. ]
         [0.  0.  0.  0.  0.5 0.  0.5 0.  0.  0. ]
         [0.  0.  0.  0.  0.  0.5 0.  0.5 0.  0. ]
         [0.  0.  0.  0.  0.  0.  0.5 0.  0.5 0. ]
         [0.  0.  0.  0.  0.  0.  0.  0.5 0.  0.5]
         [0.  0.  0.  0.  0.  0.  0.  0.  0.5 0.5]]
    """
    P = np.array([0] * M * M, dtype=np.float32).reshape(M, M)
    P[0, 0] = 0.5
    P[M - 1, M - 1] = 0.5
    for i in range(M):  # 0 indexed python
        for j in range(M):
            if j == i + 1 or j == i - 1:
                P[i, j] = 0.5

    return P

def analytical(M, a):
    pmf = []
    for k in range(1, M + 1):
        pmf.append(zipf.pmf(k, a))
    return pmf


def process_markov_data(data):
    # count frequencies
    count = Counter(s_rec)
    ## normalise to 1
    count = {k: v / sum(count.values()) for k, v in count.items()}
    return count

def plot(data):
    # obtain and plot analytical solution
    ana = analytical(M, a)
    fig = plt.figure()
    ax1 = plt.subplot(121)
    ax1.set_title('Analytical')
    ax1.stem(ana, linefmt='k-.', markerfmt='ko')
    ax1.set_xlim(-1, M)
    ax1.set_ylim(0, 1)

    ax2 = plt.subplot(122)
    ax2.stem(data.keys(), data.values(), linefmt='k-.', markerfmt='ko')
    ax2.set_title('Simulated')
    ax2.set_xlim(-1, M)
    ax2.set_ylim(0, 1)
    seaborn.despine(fig=fig, top=True, right=True)

    plt.show()



if __name__ == '__main__':
    n_iterations = 100000   # number of iterations
    burnin = 1000           # number of iterations discarded
    M = 10                  # range of M for k
    a = 2                   # Zipf parameter a

    P = transition_prob()   # get transition probability
    si = 1                  # initial state for markov chain

    s_rec = []              # initialize a vector for storing the results

    for iteration in range(n_iterations):
        # retrieve probabilities relevant for state i and propose state j
        pmf = P[si,]
        sj = np.random.choice(range(pmf.shape[0]), p=pmf)

        # compute the acceptance probability
        aij = min(1, sj * P[si, sj] / si * P[si, sj])

        # generare a random number r~U(0, 1)
        r = np.random.uniform()

        # if the random number r is smaller than acceptance probability
        # accept and go to state j
        if r < aij:
            si = sj

        # Only record transitions after burn-in phase
        if iteration > burnin:
            s_rec.append(si)

    # turn simulated data into probabilities
    s_rec = process_markov_data(s_rec)

    # plot simulated vs analytical
    plot(s_rec)

Results

Questions

Why does my simulated PMF not resemble the Analytical PMF?

Edit: Implementing the suggestions in the comments

Just to follow up from for anybody else reading this, the line:

    aij = min(1, sj * P[si, sj] / si * P[si, sj])

should be

    aij = min(1, (si+1)**a / (sj+1)**a)

Results with the change

Answer 1

Your acceptance probability in the code does not seem consistent with what you wrote above. Also keep in mind that since the indexing is zero-based on Python, if you want to do any calculations involving the actual value of the state $\texttt{si}$, you should use $\texttt{si + 1}$.