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

## 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

• Your acceptance probability does not seem consistent with what you wrote above. – kccu May 3 '19 at 13:32
• 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$. – kccu May 3 '19 at 13:43
• Perfect, thanks for pointing that out. Feel free to post an answer and I'll accept. – CiaranWelsh May 3 '19 at 14:45

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