Island Hopping with Metropolis Algorithm

Question

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 burn_in, 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?

Answer 1

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 \begin{equation} r=\frac{p(x')}{p(x)}\,\frac{g(x'\rightarrow x)}{g(x \rightarrow x')}. \end{equation} 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