Algorithm 3.1 MCMC for Dice Rolls

import numpy as np
import random
import matplotlib.pyplot as plt

# target density
f = [0,0,1,2,3,4,5,6,5,4,3,2,1]
nTrials = 5000; #20000
d = [0]*nTrials
x = 5 # starting state
for i in range(nTrials):
  # do the proposal
  u = random.random()
  if( x == 2):
    if( u < 0.5):
      y = 3
    else:
      y = 2
  elif( x == 12):
    if( u < 0.5):
      y=11
    else:
      y = 12
  else:
    if( u < 0.5):
      y = x+1
    else:
      y = x-1
  # now do the acceptance
  h = min(1,f[y]/f[x])
  u = random.random()
  if( u < h ):
    x=y #move to y, else stay in x
  d[i] = x

plt.hist(d,bins = 11)
plt.savefig("tentfig.png")
plt.show()

Algorithm 3.2 Annealing the Ising Model

import numpy as np
import matplotlib.pyplot as plt

def flipk(x,k):
  y = [0]*len(x)
  for i in range(len(x)):
    y[i] = x[i]
  y[k] = -x[k]
  return y

nTrials=2000; burnin=3000;
d = [0]*nTrials
beta=.01; M=0;
#beta=1; M=0;
x = 2*np.random.randint(0,2,4)-1 #random 4-vector of +/-1
Hx = -(x[0]*x[1]+x[2]*x[3]+x[0]*x[2]+x[1]*x[3]);
for t in range(burnin+nTrials):
  k = np.random.randint(0,4) # index to flip
  y = flipk(x,k)
  Hy = -(y[0]*y[1]+y[2]*y[3]+y[0]*y[2]+y[1]*y[3]);
  h = min(1,np.exp(-beta*(Hy-Hx)));
  #print(x)
  #print(y)
  #print(
#"Hx "+str(Hx)+" Hy "+str(Hy)+", beta*deltH "+str(beta*(Hy-Hx))+" h="+str(h))
  if( np.random.rand() < h ):
    x = y
    Hx = Hy
  if( t > burnin):
    Mx = x[0]+x[1]+x[2]+x[3];
    #print("Mx= "+str(Mx))
    M = M + Mx; # to obtain avg mag moment
    d[t-burnin] = Mx;
 # compute avg moment per spin
M = (M/float(nTrials))/4.0;
print(M)

def countElements(seq) -> dict:
  hist = {}
  for i in seq:
    hist[i] = hist.get(i, 0) + 1
  return hist
counts = countElements(d)

plt.bar(counts.keys(),counts.values())
plt.show()

Algorithm 3.3 Hastings Modification of MCMC

import numpy as np
import random
import matplotlib.pyplot as plt

def g(x,y):
  if( x == 2 or x == 12):
    return 0.5
  if( y == x+1 ):
    return 0.67
  else:
    return 0.33

def run(hastings):
  d = [0]*nTrials
  x = 5 # starting state
  nAccept = 0
  for i in range(nTrials):
    # do the proposal
    u = random.random()
    if( x == 2):
      if( u < 0.5):
        y = 3
      else:
        y = 2
    elif( x == 12):
      if( u < 0.5):
        y=11
      else:
        y = 12
    else:
      if( u < 0.67):
        y = x+1
      else:
        y = x-1
    gxy = g(x,y)
    gyx = g(y,x)
    # now do the acceptance
    if( hastings ):
      h = min(1,f[y]*gyx/(f[x]*gxy))
    else:
      h = min(1,f[y]/f[x])
                                 
    u = random.random()
    if( u < h ):
      x=y #move to y, else stay in x
      nAccept += 1
    d[i] = x
  return d,nAccept

# target density
f = [0,0,1,2,3,4,5,6,5,4,3,2,1]
nTrials = 20000

"""
d,nAccept = run(True)  #with hastings
print("acceptanceRate= "+str(nAccept/float(nTrials)))

plt.hist(d,bins=11,density=True)
plt.title('Metropolis-Hastings Dice Roll Calculation')
plt.savefig('mhDiceRollfig.png')
"""

#fig,axs = plt.subplots(nrows=1, ncols=2)
#axs[0].hist(d,bins = 11)

d,nAccept = run(False)  #w/o hastings
print("acceptanceRate= "+str(nAccept/float(nTrials)))
#axs[1].hist(d,bins = 11)

plt.hist(d,bins=11,density=True)
plt.title('Metropolis Dice Roll Calculation')
plt.savefig('metroDiceRollfig.png')

plt.show()

Algorithm 3.4 MCMC for Graph Matching

import numpy as np
import random
import matplotlib.pyplot as plt

# initialize the first two rows of E to the given graph, a 3x24 matrix
# 1st item in a column is one end of the edge, 2nd is the other end
# initialize third row to 0, start with an empty match set
E = [[0]*14, [0]*14, [0]*14]
E[0][0]=0; E[1][0]=6; E[0][1]=0; E[1][1]=8;
E[0][2]=0; E[1][2]=10; E[0][3]=1; E[1][3]=7;
E[0][4]=1; E[1][4]=10; E[0][5]=2; E[1][5]=9;
E[0][6]=2; E[1][6]=11; E[0][7]=3; E[1][7]=6;
E[0][8]=3; E[1][8]=8; E[0][9]=3; E[1][9]=11;
E[0][10]=4; E[1][10]=10; E[0][11]=4; E[1][11]=11;
E[0][12]=5; E[1][12]=7; E[0][13]=5; E[1][13]=9;

s = [0]*14; # array to count no. matches of size k
k=0; # match set size, initially empty
nTrials = 990000
for t in range(nTrials):
  i=np.random.randint(14) # select edge unif. at random
  if( E[2][i] == 1): #check if i is in the current match set
    E[2][i]=0; # remove it
    k=k-1; # decrement matching size
  else:
    l=E[0][i]; r=E[1][i]; # left, right ends of i
    kl=-1; kr=-1; #search current match set for l and r
    for j in range(14): # loop over all edges
      if( E[2][j] == 1): # this edge is in the match set
        if( l == E[0][j] or l == E[1][j]): # l is covered
          kl=j;
        if( r == E[0][j] or r == E[1][j]): # r is covered
          kr = j;
      #end if j in match set}
    #end loop over edges}
   if( kl == -1 and kr == -1): # neither l or r covered
      E[2][i]=1; #add edge i to the match set
      k=k+1; # increment match set size
    elif( kl >= 0 and kr == -1): # only left end covered
      E[2][kl]=0; E[2][i]=1; # swap edges
    elif( kl == -1 and kr >= 0): # only right end covered
      E[2][kr]=0; E[2][i]=1; # swap edges
    #end check covered
  #end check i in current match set}
  #increment the count for the current size}
  s[k] += 1;
#end loop over trials
print("number of matchings of size k = 0,1,...: "+str(s));