1
$\begingroup$

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

enter image description here

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

enter image description here

$\endgroup$
  • 2
    $\begingroup$ Your acceptance probability does not seem consistent with what you wrote above. $\endgroup$ – kccu May 3 at 13:32
  • 1
    $\begingroup$ Also keep in mind that since the indexing is zero-based in Python, if you want to do any calculations involving the actual value of the the state $\texttt{si}$, you should use $\texttt{si} + 1$. $\endgroup$ – kccu May 3 at 13:43
  • $\begingroup$ Perfect, thanks for pointing that out. Feel free to post an answer and I'll accept. $\endgroup$ – CiaranWelsh May 3 at 14:45
1
$\begingroup$

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}$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.