Island Hopping with Metropolis Algorithm

From John Kruschke's book, Chapter 7, Pg. 120 (summarised for succinctness):

A politician is constantly travelling from island to island on a chain of islands... His goal is to visit all the islands proportionally to their relative population, so that he spends proportionally more time on the more populated islands...He doesn't even know exactly how many islands there are!.. His advisers can ask the mayor of the island they are on how many people are on the island. And when the politician proposes to visit an adjacent island (to the east or west, by flipping a coin), they can ask the mayor of that adjacent island how many people there are on that island.

I wrote some Python code, where A is the island chain and run_Metropolis(A, N, start, burn_in) does the computations. The code is intentionally dumbed-down to make it more human-readable.

import random
def run_Metropolis(A, N, start, burn_in):

''' A: the distribution to be predicted (a list).
N: the number of trials / moves
start: starting index in A
burn_in: the number of  initial moves to be discarded.
'''

Posterior = [0] * len(A) # to be returned
cur_pos = start # starting index

for i in range(N):

# proposal distribution is (-1, +1) cordinates, takng into account the list boundaries
if cur_pos<len(A)-1 and cur_pos>0:
proposed_positions = [cur_pos-1, cur_pos+1]
elif cur_pos == len(A)-1:
#proposed_positions =[cur_pos-1] # if this is enforced, the posterior does not look like A
proposed_positions = [0,cur_pos-1]
elif cur_pos == 0 :
#proposed_positions = [cur_pos+1] # if this is enforced, the posterior does not look like A
proposed_positions = [len(A)-1, cur_pos+1]
proposed_pos = random.choice(proposed_positions)

# decide whether you will move
if A[proposed_pos] > A[cur_pos]: # definitely move
cur_pos = proposed_pos
else: # move with a prob proportional to how close A[cur_pos] is to A[proposed_pos]

u = random.random() # sample from uniform distribution from 0 to 1

if float(A[proposed_pos])/A[cur_pos] > u: # move to proposed
cur_pos = proposed_pos
else: # stay at cur_pos
pass

# if past brun_in point, increment the location by 1 in the posterior distribution.
if i > burn_in:
Posterior[cur_pos] +=  1

return Posterior

A = [20,30,10,50,80,20,90] # populations of the islands
posterior_A = run_Metropolis(A, 4000, 3, 1000)
# Validate the sampling works:
[n/float(sum(A)) for n in A]
> [0.0666, 0.1, 0.0333, 0.166, 0.266, 0.066, 0.3]
[n/float(sum(posterior_A)) for n in posterior_A]
> [0.0593, 0.0910, 0.0393, 0.177, 0.263, 0.0666, 0.302]


So, posterior_A does reflect the population distribution of A. What Kruschke does not make clear is that in order for this to work, the island chain needs to be cyclic, i.e. you need to be able to travel from the 7th island to the 1st one and vice versa. Otherwise, posterior_A does not look like A. You can test this by changnig the boundary proposal distributions with the commented out ones, which don't assume cyclicity.

Questions:

1.Is this because the ergodicity assumption is being violated?

2.Can we modify the acceptance criteria (move decision) to account for the lack of symmetry for the easter-most and western-most islands? (I guess this is Metropolis-Hastings)

3.If so, can you please modify the code to illustrate?

-

1. The problem is not about ergodicity

No, this is not related to ergodicity. In the chain without cycling around, one can still move from any island to any other island, and (provided that there are differences in the populations!) the chain is not periodic (because one sometimes stays put), so the chain is ergodic, hence it has a unique stationary distribution. The problem is that the stationary distribution is wrong, i.e., detailed balance is satisfied but not with the intended distribution.

2. Metropolis-Hastings

Your idea is correct. In Metropolis-Hastings, one modifies the acceptance probability to take into account the asymmetry in proposal distributions: the proposal is accepted with probability $\min(r,1)$, where $$r=\frac{p(x')}{p(x)}\,\frac{g(x'\rightarrow x)}{g(x \rightarrow x')}.$$ where $p$ is the target distribution, $x$ is the old state (island), $x'$ is the proposed state and $g(y\rightarrow z)$ is the probability (density in the continuous case) of proposing $z$ when the chain is at $y$. In your cyclic chain $g$ is always $1/2$, thus the latter factor vanishes and we are left with the plain Metropolis algorithm.

In the new chain, $g$ is 1 when moving away from either boundary island. In all other cases $g=0.5$. Thus, the latter factor in the acceptance probability must be taken into account whenever either the old island or the new island is a boundary island (the factor is $0.5$ or $2$, respectively).

3. Implementation

Only the relevant modified parts (proposal and acceptance) attached. This computes directly the ratio of $g$s for the special cases it differs from $1$.

    proposal_ratio = 1 #The ratio of proposal probabilities, change for special cases
# proposal distribution is (-1, +1) cordinates, takng into account the list boundaries
if cur_pos<len(A)-1 and cur_pos>0:
proposed_positions = [cur_pos-1, cur_pos+1]
elif cur_pos == len(A)-1:
proposed_positions =[cur_pos-1]
proposal_ratio = 0.5
elif cur_pos == 0 :
proposed_positions = [cur_pos+1]
proposal_ratio = 0.5

proposed_pos = random.choice(proposed_positions)

if proposed_pos == 0 or proposed_pos == len(A)-1:
proposal_ratio = 2

# decide whether you will move
r = (float(A[proposed_pos])/A[cur_pos]) * proposal_ratio
if r>=1: # definitely move
cur_pos = proposed_pos
else: # move with prob r
u = random.random() # sample from uniform distribution from 0 to 1
if r > u: # move to proposed
cur_pos = proposed_pos
else: # stay at cur_pos
pass

-
+1 Thank you, this is (probably) the most complete answer (adressing eveything I asked) I received on this website. –  Zhubarb Jul 17 at 12:07
Hmm, I ran the code with your modification to get this: [n/float(sum(posterior_A)) for n in posterior_A] = [0.060, 0.12, 0.05, 0.198, 0.312, 0.079, 0.178]. The issue with the last element still persists (it is supposed to be approx. 0.3). –  Zhubarb Jul 17 at 12:26
I tested it and got approx 0.3 with my modification. Check that you haven't made an error when pasting the modification. Maybe upload your modified code to, e.g., pastebin? –  Juho Kokkala Jul 17 at 12:39
sorry, my mistake. i mixed cur_pos and proposed_pos in one of the conditions. It works fine. –  Zhubarb Jul 17 at 12:47
I corrected the error in text (it seems you deleted the comment, but indeed I had mixed up the 0.5 vs. 2 cases in the text - code was correct). –  Juho Kokkala Jul 17 at 13:04