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?