1
$\begingroup$

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: Markov Chains Last 2500 out 5000 samples

$\endgroup$
7
  • $\begingroup$ Can you explain a little better what you did? For example, what's your target density? $\endgroup$ Oct 21, 2019 at 1:03
  • $\begingroup$ 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? $\endgroup$ Oct 21, 2019 at 1:06
  • $\begingroup$ Is your gfunc drawing a sample, or returning a representation of the whole function, or something else? $\endgroup$ Oct 21, 2019 at 1:10
  • $\begingroup$ Just out of curiosity, would you post a plot of the samples it's generating? Hopefully colored by $i$? $\endgroup$ Oct 21, 2019 at 1:17
  • $\begingroup$ 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). $\endgroup$ Oct 21, 2019 at 1:28

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$ Oct 22, 2019 at 1:59
  • $\begingroup$ 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$.) $\endgroup$ Oct 22, 2019 at 14:19
  • $\begingroup$ This seems weird and hard to construct though. I may have understood wrong. $\endgroup$ Oct 22, 2019 at 14:23
  • $\begingroup$ Can you just somehow choose your grid points without reference to $x_i$? Like a permanently fixed grid? $\endgroup$ Oct 22, 2019 at 14:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.