# Implementing a Hit-And-Run MCMC Sampler in Python

I'm trying to implement the Metropolisized Hit-And-Run Sampler as described in Chapter 7.2 of Chen and Schmeiser, Performance of the Gibbs, Hit-and-Run, and Metropolis Samplers https://www.jstor.org/stable/pdf/1390645.pdf

I chose the function $$g$$ to be a piecewise constant distribution function based on evaluations of the log probability along the direction of movement.

It seems that the chain does not converge to the posterior distribution. I particular both mean and standard deviation of the produced samples are smaller than the true values.

Here's my python implementation. Am I missing something? I haven't found a working version of this sampler anywhere.

import numpy as np

# test log probability
def logprob(x):
return - 0.5*np.sum((x-5.0)**2.0)

x0 = np.array([-1.0,-1.0]) # Initial point
nsteps = 50000 # Number of steps
ndim = 2 # Number of dimensions

# Piecewise constant density g
def gfunc(logprob,x0,d,lmax=1.0,npieces=2):

dl = 2. * lmax / npieces  # Compute width of piecewise constant regions
# Compute centres of those regions
centres = np.linspace(x0-(npieces-1.0)/2.0*dl*d,x0+(npieces-1.0)/2.0*dl*d,npieces)

logprobvals = list(map(logprob,centres)) # Evaluate the logprob at those points
probvals = np.exp(logprobvals) # Get the probability
gvals = probvals / np.sum(probvals) # Normalise the probability to sum to 1

piece = np.random.choice(np.arange(-int(npieces/2),int(npieces/2)),p=gvals) # Choose a random region with probability goals

gprime = gvals[int(npieces/2)+piece] # Compute g(x0+l*d)

l = np.random.uniform(piece*dl, (piece+1)*dl, size=ndim) # Uniformly sample l from that region

return l, gprime

x = np.empty((nsteps,ndim)) # Initialise array to save samples
g = np.empty((nsteps)) # Initialise array to save g evaluations
g[-1] = 0.1 # This is arbitrary
x[-1] = x0
for i in range(nsteps):
d = np.random.normal(0,1,size=ndim) # Compute unit direction vector
d /= np.linalg.norm(d)

l, gp = gfunc(logprob,x[i-1],d,lmax=5.0,npieces=100) # Sample l and get g(x0+l*d)
ratio = g[i-1]/gp # Compute g/g' ratio to be used in M-H acceptance criterion

Z = x[i-1] + l * d # Compute new point

# Check Metropolis-Hastings acceptance criterion
if np.random.uniform(0.0,1.0) <= min(1.0, ratio * np.exp(logprob(Z)-logprob(x[i-1]))): # Accept
x[i] = Z
g[i] = gp
else: # Reject
x[i] = x[i-1]
g[i] = g[i-1]


Update

Here are the produced Markov chains and samples respectively:

• Can you explain a little better what you did? For example, what's your target density? Oct 21, 2019 at 1:03
• Also, a thought experiment (and a unit test to implement): what is the integral of your $g$ on the whole real line, and what should it be? Oct 21, 2019 at 1:06
• Is your gfunc drawing a sample, or returning a representation of the whole function, or something else? Oct 21, 2019 at 1:10
• Just out of curiosity, would you post a plot of the samples it's generating? Hopefully colored by $i$? Oct 21, 2019 at 1:17
• I just added the plots. The target density is a 2D Normal distribution of mean (5,5) and unit variance (i.e. see logprob function). Oct 21, 2019 at 1:28

The problem may be that your $$g$$ depends on $$x_i$$. Your source says "We now require that $$g$$ depend on $$x_i$$ only through the line segment $$S_i$$." But your $$g$$ is basically completely flexible; if you were to plot $$\lambda$$ versus $$g(\lambda)$$, it could be different with every new point sampled.
• This make sense! I tried using different definitions of $g$ that depend on $x_{i}$ only through the line segment $S_{i}$, but I fail to find a proper working one. In particular, the main difficulty is constructing a general proposal $g$ that returns consistent values for $g(x_{i})$ in the current point using the previous line segment (i.e. $g_{i-1}(x_{i-1}+\lambda_{i-1}\vec{d}_{i-1})$) and the proposed one (i.e. $g_{i}(x_{i})$). Since the current point is the intersection of the previous line $g$ and the current $g$ one, the two of them should agree, right? Oct 22, 2019 at 1:59
• I'm not sure that's the right way to say $g$ should depend on $S_i$ but not $x_i$. I would say you want $g_{x_{i-1}, d}(y) = g_{x_i, d}(y)$ for all $y$ and $d$. (Note I keep the direction the same; it's $d$ and not $d_i$.) Oct 22, 2019 at 14:19
• Can you just somehow choose your grid points without reference to $x_i$? Like a permanently fixed grid? Oct 22, 2019 at 14:35