Algorithm 6.1 MCMC for the Posterior Distribution

import numpy as np
import matplotlib.pyplot as plt
def propose(x):
  u = delta*np.random.uniform(-1,1)
  y = x+u
 # keep it in [0,1]
  if( y<0 ):
    y = -y
  elif( y > 1 ):
    y = 2 - y
  return y
def lny(t,m):
  return m*(t-1)+0.01
def sumSquares(m):
  acc = 0
  for i in range(len(tdat)):
    acc += (ydat[i] - lny(tdat[i],m))**2
  return acc
def prior(m):
  y = (1/((2*np.pi)**0.5*sigma_m))*np.exp(-0.5*((m-mu_m)/sigma_m)**2)
  return y
def prob_obsGiven_m(m):
  y = (1/((2*np.pi)**0.5*sigma_l)**5)*np.exp(-sumSquares(m)/(2*sigma_l**2))
  return y
def f(m):
  y = prob_obsGiven_m(m)
  y *= prior(m)
  return y
delta = 0.05
tdat = [1,5,10,15,20] #in days
ydat = [.01,0.1,0.9,4,9] #in ln(increase in gm)
## normal prior with a known sigma
mu_m = 0.20
sigma_m = 0.04
## normal likelihood, fixed sigma
sigma_l = 0.6
trace = []
x = np.random.normal(mu_m,sigma_m)
fx = f(x)
burnin = 1000
for i in range(100000):
  y = propose(x)
  fy = f(y)
  h = fy/fx
  u = np.random.uniform()
  if( u < h ):
    x = y
    fx = fy
  if( i > burnin):
    trace.append(x)
meanValue = np.mean(trace)
variance = np.var(trace)
print("mean",meanValue,"var",variance)
plt.hist(trace,bins=60,density=True)
plt.show()

Algorithm 6.2 Pymc3 for the Coin Tossing Problem

import numpy as np
import pymc3 as pm
import scipy.stats as stats
from scipy.stats import beta
import matplotlib.pyplot as plt

data=np.array([1,1,1,1,1,1,1,0,0,0]) #our observation
#print(sum(data)/float(len(data)))
with pm.Model() as model:
    # set a Beta prior on the Bernoulli parameter
    theta=pm.Beta('theta', alpha=3, beta=3)
    # set a Bernoulli likelihood for the data
    y = pm.Bernoulli("obs", p = theta, observed = data)
    # inference step
    step=pm.Metropolis()
    trace=pm.sample(18000, step=step,random_seed=123)
    burned_trace=trace[1000:]

plt.hist(burned_trace["theta"], bins=40, histtype='stepfilled',
density=True,color="darkgrey")
x=np.arange(0,1.04,0.04)
y = beta.pdf(x,10,6)  #plot beta(10,6)
plt.plot(x, y, color='#FF7F00')
xmax = 0.642
plt.vlines(xmax,0.3, 3.5, linestyle='--', label=r"MAP estimate", color='red')
plt.text(xmax+0.025,0.40, "MAP", color='red', fontsize=12)
plt.text(xmax-0.04,0.08,str(xmax), color='red', fontsize=12)
plt.text(0.82, 2.1, "beta(10,6)", color='black')
plt.hlines(2.15,0.755,0.80)
plt.title(
"Pymc3 MCMC histogram for heads probability theta",pad=10,fontsize=14)
plt.xlabel("theta")
#plt.savefig("pymcSeventenfig.png")

# Plot a simple line chart without any feature
plt.plot(x, y)
plt.show()

Algorithm 6.3 Pymc3 for the Wheat Growth Problem

import numpy as np
import pymc3 as pm
import arviz as az
import matplotlib.pyplot as plt
#import scipy.stats as stats
#from scipy.stats import norm

mu_m = 0.20; sig_m = 0.04;
tdat = np.array([1.,5.,10.,15.,20.]) #in days
ydat = np.array([0.01,0.1,0.9,4.,9.]) #in ln(increase in gm)
sigma_l = 0.6 # normal likelihood, fixed sigma

with pm.Model() as model:
    # set a normal prior on the parameter m
    m = pm.Normal('m',mu_m,sig_m)
    # set the likelihood for the data
    lny = 0.01+m*(tdat-1)
    y = pm.Normal('y',mu=lny,sd=sigma_l,observed=ydat)
    # inference step
    step=pm.Metropolis()
    trace=pm.sample(20000, step=step,random_seed=123)
with model:
    az.plot_posterior(trace,var_names=['m'])
plt.title(
"Pymc3 MCMC histogram for Growth Rate m",pad=10,fontsize=14)
plt.savefig("pymcGrowthfig.png")
plt.show()





Algorithm 6.4 Gibbs Sampling

import numpy as np
import matplotlib.pyplot as plt

def fXY(y):
  return np.random.normal(rho*y,var)
def fYX(x):
  return np.random.normal(rho*x,var)
np.random.seed(227)
rho = 0.5
var = 1-rho**2
x=[]; y=[]
x.append(1); y.append(1)
for i in range(20):
  x.append(fXY(y[i])); #sample x given current y
  y.append(fYX(x[i+1])); #sample y given new x

xa = np.linspace(-3,3,100)
ya = np.linspace(-3,3,100)
X,Y = np.meshgrid(xa,ya)
Z = (1/(2*np.pi*np.sqrt(var)))*np.exp((-1/(2*var))*(X**2-2*rho*X*Y+Y**2))
#do axes
plt.hlines(0,-3,3, color='k'); #x-axis
plt.vlines(0,-3,3, color='k'); #y-axis
plt.contour(X, Y, Z, colors='black')
#plot a few points
plt.scatter(x,y,color='red')
#plot two steps
plt.scatter(x[0],y[0],color='orange')
plt.scatter(x[1],y[0],color='orange');
plt.hlines(y[0],x[0],x[1],color='blue');
plt.scatter(x[1],y[1],color='orange');
plt.vlines(x[1],y[0],y[1],color='blue')
plt.suptitle("Density plot generated by Gibbs sampling")
#plt.savefig("gibbsfig")
plt.show()





Algorithm 6.5 Sample Mean Estimator

import numpy as np
import matplotlib.pyplot as plt
def f1(x): # on [0,infty)
  return 1/((x+1)**2)
def f2(x):
  return -2+(x-2)**2 + 1/np.sqrt(x)
def exponDen(x,lam):
  return lam*np.exp(-lam*x)
def exponSample(lam):
  U = np.random.uniform()
  return (-1/lam)*np.log(1-U)
nTrials = 1000
lam = 0.1
acc = 0; # accumulator
for i in range(nTrials):
  U = exponSample(lam)
  acc += f1(U)/exponDen(U,lam)
print("est=",acc/float(nTrials))
acc = 0
for i in range(nTrials):
  U = np.random.uniform(0,4)
  acc += 4*f2(U); # b-a=4
print("est=",acc/float(nTrials))
fig = plt.figure(figsize=(6.5,2.5))
x = np.linspace(0,40,100)
y = f1(x)
ax0=fig.add_subplot(121)
ax0.plot(x,y)
y = exponDen(x,lam)
ax0.plot(x,y)
plt.title("f(x)={\scriptsize\$}1/(x+1)\^{}2{\scriptsize\$}")
x = np.linspace(0,4,100)
y = f2(x)
ax1=fig.add_subplot(122)
ax1.plot(x,y)
y = [0.25]*len(x)
ax1.plot(x,y)
plt.title("f(x)={\scriptsize\$}-2+(x-2)\^{}2+\
(1/\\sqrt{\scriptsize\{}x{\scriptsize\}}){\scriptsize\$}")
plt.show()






Algorithm 6.6 Estimating a Double Integral

import numpy as np
import matplotlib.pyplot as plt

def fcn(x,y):
  return np.sqrt(1-(x*x+y*y))

def q(x,y):
  return 4/np.pi

def sample1(): #uniform over [0,1]x[0,1]
  x = np.random.uniform()
  y = np.random.uniform()
  return x,y

def sample2(): #uniform over x^2+y^2<1 with rejection
  while( True):
    x,y = sample1()
    if( x*x+y*y < 1):
      break
  return x,y

def sample3(): #uniform over x^2+y^2<1 directly
  u1 = np.random.uniform()
  u2 = np.random.uniform()
  r = np.sqrt(u1); th = np.pi*u2/2
  x = r*np.cos(th); y = r*np.sin(th)
  return x,y

def integral(nTrials):
  acc = 0
  for i in range(nTrials):
    x,y = sample3()
    if( x*x + y*y > 1): #only used for sample2()
      nTrials -= 1
      continue
    acc += fcn(x,y)/q(x,y)
  est = acc/float(nTrials)
  error = abs(est - np.pi/6)
  return est,error

n = [1000,3000,5000,7000,9000,110000]
err = [0]*len(n)
for i in range(len(n)):
  navg = 20
  acc = 0
  for j in range(navg):
    est,error = integral(n[i])
    acc += error
  err[i] = np.log(acc/navg)
  print(".",end="",flush=True)
  #pcError = 100*err[i]/(np.pi/6)
  #print("est=",est,"cf 1, pcError=",pcError)
ln = np.log(n)
plt.plot(ln,err)
plt.xlabel("log(n)")
plt.ylabel("log(error)")
plt.title("Error chart for $f(x)=\sqrt{1-(x*x+y*y)}$")
a,b = np.polyfit(ln,err,1)
print("a=",a)
ye = []
y5 = []
for i in range(len(n)):
  ye.append(a*ln[i]+b)
  y5.append(-0.5*ln[i]+b)
plt.plot(ln,ye)
plt.plot(ln,y5)
plt.show()

#pcError = 100*error/(np.pi/6)
#print("est=",est,"cf pi/6=",np.pi/6,"pcError=",pcError)





Algorithm 6.7 Estimating the Triple Integral for the 4-Sphere

import numpy as np
import matplotlib.pyplot as plt

def rho(x,y,z):
  return 6/np.pi

def sample3():
  r = (np.random.uniform())**0.333
  phi = np.arccos(np.random.uniform())
  theta = np.random.uniform(0,np.pi/2.)
  x = r*np.sin(phi)*np.cos(theta)
  y = r* np.sin(phi)*np.sin(theta)
  z = r*np.cos(phi)
  return x,y,z


def integral(nTrials):
  acc = 0
  for i in range(nTrials):
    x,y,z = sample3()
    acc += fcn2(x,y,z)/rho(x,y,z)
  est = acc/nTrials
  error = abs(est - exact)
  return est,error

exact = 0.308552;
n = [1000,3000,5000,7000,9000,110000]
err = [0]*len(n)
for i in range(len(n)):
  navg = 20; #number to average error over
  acc = 0
  for j in range(navg):
    est,error = integral(n[i])
    acc += error
  err[i] = np.log(acc/navg)
  print(".",end="",flush=True)
  pcError = 100*error/exact
  print("n=",n[i],"est=",est,"cf ",exact,"pcError=",pcError)
ln = np.log(n)
plt.plot(ln,err)
plt.xlabel("log(n)")
plt.ylabel("log(error)")
plt.title("Error chart for $f(x)=\sqrt{1-(x*x+y*y+z*z)}$")
a,b = np.polyfit(ln,err,1)
print("a=",a)
ye = []
y5 = []
for i in range(len(n)):
  ye.append(a*ln[i]+b)
  y5.append(-0.5*ln[i]+b)
plt.plot(ln,ye)
plt.plot(ln,y5)
plt.show()





Algorithm 6.8 Control Vaiate Method

##
# use cauchy on [-2,2] to approximate normal [-2,2]
##

import numpy as np
import matplotlib.pyplot as plt

def h1(x): #uniform on [-2,2]
  return 0.25

def samp1():
  u = np.random.uniform(-2,2)
  return u

def normal(x,sigma):
  a = 1/(sigma*np.sqrt(2*np.pi))
  b = np.exp(-0.5*(x/sigma)**2)
  return a*b

def cauchy(x):
  return 1/(np.pi*(1+x*x))

sigma = 1

nTrials = 4000
acc = 0
for i in range(nTrials):
  u = samp1()
  acc += normal(u,sigma)/h1(u)
acc /= float(nTrials)
est = acc;
error = abs(est - 0.9545)
percentError = 100*(error/0.9545)
#print("straight est=",est,"cf. ",0.9545,"error=",error)
print("straight est=",est,"cf. ",0.9545,"percenterror=",percentError)

G = (2/np.pi)*np.arctan(2); #0.70483
acc = 0
for i in range(nTrials):
  u = samp1()
  acc += (normal(u,sigma) - cauchy(u))/h1(u)
acc /= float(nTrials)
est = G + acc;
error = abs(est - 0.9545)
percentError = 100*(error/0.9545)
#print("est=",est,"cf. ",0.9545,"error=",error)
print("est=",est,"cf. ",0.9545,"percenterror=",percentError)

"""
x = np.linspace(-2,2,100)
y = normal(x,sigma)
plt.plot(x,y,label="normal density")
y = cauchy(x)
plt.plot(x,y,label="Cauchy density")
y = [h1(x)]*len(x)
plt.plot(x,y,label="uniform density")
plt.legend()
plt.suptitle("The Cauchy density as control variate for the normal");
#plt.savefig("cauchyasCVfig.png")
plt.show()
"""





Algorithm 6.9 Importance Sampling Method in 1d

##
# use normal to estimate integral of cauchy on [-2,2]
##

import numpy as np
import matplotlib.pyplot as plt

def qu(x): #uniform on [-2,2]
  return 0.25

def sampU(): #uniform on [-2,2=]
  u = np.random.uniform(-2,2)
  return u

def qn(x): #normal density on [-2,2]
  a = 1/(sigma*np.sqrt(2*np.pi))
  b = np.exp(-0.5*(x/sigma)**2)
  c = 0.9545; # integral normal(x) [-2,2]
  return a*b/c

def sampN(): #from normal on [-2,2]
  while( True ):
    x = np.random.normal()
    if( x<2 and x > -2 ):
      break
  return x

def normal(x,sigma):
  a = 1/(sigma*np.sqrt(2*np.pi))
  b = np.exp(-0.5*(x/sigma)**2)
  return a*b

def cauchy(x):
  return 1/(np.pi*(1+x*x))

delta = 0.1
def mcmcNormal(x): #mcmc normal samples from [-2,2]
  y = x-delta + 2*delta*np.random.uniform()
  if( y>2 ):
    y = 4-y
  elif( y < -2 ):
    y = -4-y
  h = normal(y,sigma)/normal(x,sigma)
  u = np.random.uniform()
  if( u < h):
    x = y
  return x

sigma = 1
exact = (2/np.pi)*np.arctan(2)

nTrials = 40000
Q = 0
for i in range(nTrials):
  x = np.random.uniform(-2,2)
  Q += normal(x,sigma)
Q = 4*Q/float(nTrials)
print("Q=",Q,"cf. 0.9545")


nTrials = 400
acc = 0
for i in range(nTrials):
  u = sampU()
  acc += cauchy(u)/qu(u)
est = acc/float(nTrials)
error = abs(est - exact)
pcError = 100*error/exact
print(nTrials,"trials uniform sampling:",'%.6f'%est,"cf. ",
round(exact,6),"percenterror=",round(pcError,6))

acc = 0
for i in range(nTrials):
  u = sampN()
  acc += cauchy(u)/qn(u)
est = acc/float(nTrials)
error = abs(est - exact)
pcError = 100*error/exact
print(nTrials,"trials normal sampling:",'%.6f'%est,"cf. ",
round(exact,6),"percenterror=",round(pcError,6))

nTrials = 4000
acc = 0
x = 0
for i in range(nTrials):
  x = mcmcNormal(x)
  acc += cauchy(x)/qn(x)
est = acc/float(nTrials)
error = abs(est - exact)
pcError = 100*error/exact
print(nTrials,"trials mcmc w/exact Q:",'%.6f'%est,"cf. ",
round(exact,6),"percenterror=",round(pcError,6))

acc = 0; Q = 0
x = 0
for i in range(nTrials):
  x = mcmcNormal(x)
  acc += cauchy(x)/normal(x,sigma)
  u = np.random.uniform(-2,2)
  Q += normal(u,sigma)
Q = 4*Q/float(nTrials); #b-a = 4
print("Q=",Q,"cf. 0.9545")
est = Q*acc/float(nTrials)
error = abs(est - exact)
pcError = 100*error/exact
print(nTrials,"trials mcmc w/estimated Q:",'%.6f'%est,"cf. ",
round(exact,6),"percenterror=",round(pcError,6))





Algorithm 6.10 Importance Sampling Method in 2D

#
# with 30000 iterates: slant est=0.7644 cf 0.7627, pcError= 0.2261
#
import numpy as np
import matplotlib.pyplot as plt

def fcnR(x,y): # defined on triangle
  return -1+1/np.sqrt(x*x+y*y)

def hUniform(x,y): #(1/2)*base*height = 1
  return 1

def sampleUniform():
  u1 = np.random.uniform()
  u2 = np.random.uniform()
  X = np.sqrt(2)*(1 - np.sqrt(u1))
  Y = (np.sqrt(2)-X)*u2
  return X,Y

def hSlant(x,y):
  return np.sqrt(2)-x-y

def qSlant(x,y):
  c = 6/2**1.5
  return c*hSlant(x,y)

def sampleSlant():
  u1 = np.random.uniform()
  u2 = np.random.uniform()
  X = np.sqrt(2)*(1 - u1**0.334)
  Y = (np.sqrt(2)-X)*(1-u2**0.5)
  return X,Y

delta = 0.1
def sampM(x,y):
  while( True ):
    ux = np.random.uniform()
    uy = np.random.uniform()
    xn = x-delta + 2*delta*ux
    yn = y-delta + 2*delta*uy
    if( xn < 0):
      xn = -xn
    if( yn < 0):
      yn = -yn
    if( xn+yn > np.sqrt(2)):
     continue
    break; #new point in D
  h = hSlant(xn,yn)/hSlant(x,y)
  u = np.random.uniform()
  if( u < h):
    x = xn; y = yn;
  return x,y

## mmm ##
n = 30000;

x = [0.01,0,np.sqrt(2),0]
y = [0,np.sqrt(2),0,0.01]
#plt.plot(x,y)

uni = 0;
slant = 1;
mcmc = 0;
if( uni+slant+mcmc != 1 ):
  print("exactly one of uni,slant,mcmc = 1 please")
  quit()

acc=0; Q = 0;
x = 0.7; y = 0.7;
for i in range(n):
  if( uni == 1 ):
    x,y = sampleUniform()
    acc += fcnR(x,y)/hUniform(x,y)
  elif( slant== 1 ):
    x,y = sampleSlant()
    acc += fcnR(x,y)/qSlant(x,y)
  else:
    x,y = sampM(x,y)
    acc += fcnR(x,y)/hSlant(x,y)
    u,v = sampleUniform()
    Q += hSlant(u,v); #area of triange = 1
  #plt.scatter(x,y,color='black',s=8)
  if( i%1000 == 0 ):
    print(".",end="",flush=True)
if( mcmc == 1 ):
  Q /= float(n)
  print("Q=",Q,"cf ",2**1.5/6)
else:
  Q = 1
est = Q*acc/float(n)
exact = 0.7627
error = abs(est - exact)
pcError = 100*error/exact
if( uni == 1 ):
  print("unif est=",end="")
if( slant == 1 ):
  print("slant est=",end="")
else:
  print("mcmc est=",end="")
print(round(est,4),"cf 0.7627, pcError=",round(pcError,4))
#plt.show()





Algorithm 6.11 Estimating Against a Weight Function

import numpy as np
def f(x):
  return x*x
def expWeight(x):
  return np.exp(-0.5*x)
delta = 1
def sampM(x): # Metropolis sampling
  u = np.random.uniform()
  y = x-delta + 2*delta*u
  if( y < 0 ):
    y = -y
  h = expWeight(y)/expWeight(x)
  u = np.random.uniform()
  if( u < h):
    x = y
  return x
nTrials = 100000
acc = 0; W = 0;
x = 1; # MCMC starter
for i in range(nTrials):
  x = sampM(x)
  acc += f(x)
  u = np.random.uniform(0,10)
  W += expWeight(u)
W = 10*W/float(nTrials); #b-a=10
print("W=",W,"cf. 2")
est = W*acc/float(nTrials)
exact = 16
error = abs(est - exact)
pcError = 100*error/exact
print("est=",est,"cf ",
exact," pcError",pcError)





Algorithm 6.12 Using an Arbitrary Non-negative Function

import numpy as np
def fcnR(x,y): #defined on triangle
  return -1+1/np.sqrt(x*x+y*y)
def hUniform(x,y):
  return 1; #(1/2)*base*height = 1
def sampleUniform():
  u1 = np.random.uniform()
  u2 = np.random.uniform()
  X = np.sqrt(2)*(1 - np.sqrt(u1))
  Y = (np.sqrt(2)-X)*u2
  return X,Y
def hSlant(x,y):
  return np.sqrt(2)-x-y
delta = 0.1
def sampleM(x,y): #MCMC sample
  while( True ):
    ux = np.random.uniform()
    uy = np.random.uniform()
    xn = x-delta + 2*delta*ux
    yn = y-delta + 2*delta*uy
    if( xn < 0):
      xn = -xn
    if( yn < 0):
      yn = -yn
    if( xn+yn > np.sqrt(2)):
      continue; # outside D
    break; #new point in D
  h = hSlant(xn,yn)/hSlant(x,y)
  u = np.random.uniform()
  if( u < h): #accept new
    x = xn; y = yn;
  return x,y
nTrials = 200000;
acc=0; Q = 0;
x = 0.7; y = 0.7; #starter points
for i in range(nTrials):
  x,y = sampleM(x,y)
  acc += fcnR(x,y)/hSlant(x,y)
  u,v = sampleUniform()
  Q += hSlant(u,v)/hUniform(u,v)
Q /= float(nTrials)
print("Q=",Q,"cf ",2**1.5/6)
est = Q*acc/float(nTrials)
exact = 0.7627
pcError = 100*abs(est - exact)/exact
print("mcmc est=",end="")
print(round(est,4),"cf 0.7627,pcError=",round(pcError,4))