When a Random Walk Monte Carlo move is proposed by IMP::core::MonteCarlo, Model::update() is called during scoring function evaluation. For variables with ScoreStates, this causes their score states to be enforced before evaluation. When the score state enforces a boundary constraint, this can become problematic, as it modifies the proposal transition density. The proper way to handle this would be with rejection sampling, rejecting all moves that would take the variable outside of its boundary.

Instead, for IMP::isd::Nuisance, and potentially other constrained variables, whenever a move outside the boundary is proposed, it is moved to the boundary for evaluation. This is analogous to placing a magnet at the boundary; the transition density exactly at the boundary is much higher than it ought to be. When the score at the boundary is undefined, as with an IMP::isd::JeffreysRestraint, then this move will be rejected, and we effectively get rejection sampling. When it is defined, then the transition is no longer symmetric, and the resulting Markov chains do not sample the right distribution. In general, they may become frozen at the boundary for long periods of time, biasing the resulting sample.

Below is an MWE adapted from this hmc example. Here we sample a Nuisance with lower bound of 0 using a half-normal distribution (an IMP::isd::Gaussian with mean 0 and standard deviation 1).

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

import IMP
import IMP.core
import IMP.isd

# Seeds
IMP.random_number_generator.seed(42)
np.random.seed(8675309)

# Config params
nsteps = 10000
nwarmup = nsteps
nsample = nsteps

# Set-up model
m = IMP.Model()
p = IMP.isd.Nuisance.setup_particle(IMP.Particle(m))
p.set_nuisance_is_optimized(True)
p.set_lower(0)
r = IMP.isd.GaussianRestraint(p.get_particle(), 0, 1)
sf = IMP.core.RestraintsScoringFunction([r])
p.set_nuisance(abs(np.random.normal()))

# Set-up RWMC
mvr = IMP.core.NormalMover(
    m, p.get_particle_index(), [p.get_nuisance_key()], 2.0
)
mc = IMP.core.MonteCarlo(m)
mc.set_return_best(False)
mc.set_scoring_function(sf)
mc.add_movers([mvr])

# Sample RWMC
mc.optimize(nwarmup)
mc.reset_statistics()
xs_rwmc = []
naccept = 0
for n in range(nsample):
    mc.optimize(1)
    naccept += mc.get_number_of_accepted_steps()
    xs_rwmc.append(p.get_nuisance())

# Sample exact
xs_exact = np.fabs(np.random.normal(size=nsample))

Exploring the results

>>> naccept / nsample
0.6952
>>> hist1 = sns.distplot(xs_exact, kde=False, label="Exact")
>>> hist2 = sns.distplot(xs_rwmc, kde=False, label="RWMC in IMP")
>>> plt.legend()
>>> plt.show()

As expected, IMP's samples using RWMC differ substantially from the true distribution:

Examining the trace closer, we can see the predicted freezing behavior:

>>> plt.plot(xs_rwmc[:300])
>>> plt.show()

There are several ways to handle this. One is to just not support restraints that are defined at constraint boundaries. This is very restrictive. Half-normal, for example, is a very common default prior for standard deviations and is far more sensible than a Jeffreys prior (which has infinite probability mass!).

An alternative would be to find a way to reject such moves instead of moving them to the boundary. This has its own problems, however. As I covered in my group meeting, constraints like these essentially introduce an infinite potential that can make it difficult for the Markov chain to sample the probability mass near the constraint. In practice, it will generally exhibit freezing behavior to a lesser extent, and posterior expectations will have poor convergence, if they converge at all in finite time.

A better solution was proposed in #1007 and consists of transforming the variables to an unconstrained space for sampling, making the necessary correction to the density. Currently, HMCUtilities.jl provides these transformations and corrections for the hmc module, but it would be nice if these were standard in IMP.