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 =  * 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]
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.
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?