metropolis_island_hopping.py

import random

# Codility (Prefix sum) inspired - MCMC
def get_posterior_from_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 values 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

        proposal_ratio = 1 #The ratio of proposal probabilities, change when moving to or from 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
            proposal_ratio = 0.5 # moving away from a boundary
            #proposed_positions = [0,cur_pos-1] This makes the islands cyclic, you can jump from the last island to the first.
        elif cur_pos == 0 :
            proposed_positions =[cur_pos+1]
            proposal_ratio = 0.5 # moving away from a boundary
            #proposed_positions = [len(A)-1, cur_pos+1] # if proposed_positions =[cur_pos+1], the posterior does not look like A

        proposed_pos = random.choice(proposed_positions)

        #  if moving to a boundary (0th or the nth element in the list)
        if proposed_pos == 0 or proposed_pos == len(A)-1:
            proposal_ratio = 2 # prob of moving away from boundary to proposed_pos (1) is twice the prob of moving into cur_pos from proposed_pos (1/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 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 r > u: # move to proposed
                 cur_pos = proposed_pos 
  
        # if past burn_in point, increment the location by 1 in the posterior distribution.
        if i > burn_in: 
            Posterior[cur_pos] +=  1

    return Posterior

A = [2,3,1,5,8,2,9] # distribution to be predicted.
posterior_A = get_posterior_from_Metropolis(A, 4000, 3, 1000)

# Validate the sampling works:
print [n/float(sum(A)) for n in A]
print [n/float(sum(posterior_A)) for n in posterior_A]