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$:
Generate a proposal state, $j$, according to the proposal chain P.
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)