Suppose you can calculate a continuous distribution’s weight at any point but cannot sample from it directly. This happens with a posterior: multiply a prior’s weight by the likelihood of the observed evidence, and you obtain a weight function without necessarily knowing its normalizing constant. Rejection sampling may require details of the posterior’s shape that you do not yet know. The Metropolis algorithm instead uses relative weights to generate a sequence of samples. section 16.3 section 16.4
Keep one sample as the current state and generate the next as follows:
- For the first sample, draw from an initial distribution and return that value. section 16.4
- Otherwise, draw a proposed value from a distribution centered on the current value. Divide the proposed value’s target weight by the current value’s target weight. section 16.4
- Draw a standard uniform value. If it is no greater than the ratio, return the proposal as the new current value; otherwise, return the current value again. section 16.4
A proposal with greater weight is always accepted. One with one-tenth the current weight is accepted one-tenth of the time. Because the decision uses a ratio, a common normalization factor cancels. A rejected proposal does not trigger another attempt: the current value is emitted again, so each generated sample takes a constant amount of work rather than an unbounded rejection loop. The resulting sequence is a Markov process, with each transition depending on its current value. section 16.4 section 16.6
Listing 16.7 The Metropolis algorithm
public sealed class Metropolis: IContinuousDistribution
{
private readonly IEnumerator<double> enumerator;
private readonly Func<double, double> target;
public static Metropolis Distribution(
Func<double, double> target,
IDistribution<double> initial,
Func<double, IDistribution<double>> proposal)
{
var markov = Markov<double>.Distribution(initial, transition);
return new Metropolis(target, markov.Sample().GetEnumerator());
IDistribution<double> transition(double current)
{
double proposed = proposal(current).Sample();
return new Flip<double>(proposed, current,
target(proposed) / target(current));
}
}
private Metropolis(
Func<double, double> target, IEnumerator<double> enumerator)
{
this.enumerator = enumerator;
this.target = target;
}
public double Sample()
{
enumerator.MoveNext();
return enumerator.Current;
}
public double Weight(double d) => target(d);
}
The C# implementation makes that transition explicit: sample `proposal(current)`, then use a weighted coin flip to choose between `proposed` and `current`. It supplies the transition and initial distribution to `Markov<double>.Distribution`. section 16.4
The initial and proposal distributions still need sensible choices. For posterior sampling, the prior can supply the initial value because it is expected to share the posterior’s support. A normal proposal centered on the current value can use a standard deviation chosen for the process’s expected range. The demonstrated implementation is one-dimensional, although the algorithm is particularly suited to multidimensional distributions. section 16.5 section 16.4
Listing 16.8 Completing the implementation of continuous posterior sampling
static class Extensions
{
public static Func<B, IContinuousDistribution> Posterior<B>(
this IContinuousDistribution prior,
Func<double, IWeightedDistribution<B>> likelihood,
double stddev = 1.0)
{
return (B evidence) =>
{
Func<double, double> posteriorPDF = a =>
prior.Weight(a) * likelihood(a).Weight(evidence);
Func<double, IContinuousDistribution> proposal =
current => Normal.Distribution(current, stddev);
return Metropolis.Distribution(posteriorPDF, prior, proposal);
};
}
}
The `Posterior` implementation combines `prior.Weight(a) * likelihood(a).Weight(evidence)` into a target function, uses the prior for initialization, and passes a current-centered normal proposal to `Metropolis.Distribution`. section 16.5
Updating a coin-fairness model
Model a factory coin’s probability of heads with a normal prior centered at 0.5, clamped to the interval from 0 to 1. A binomial likelihood models the number of heads in ten flips. After observing nine heads, sample the posterior and compare its histogram with the prior’s: the observed flips shift the distribution toward higher probabilities of heads. section 16.5
Check how likelihood weights are defined before combining them with a prior. If likelihood distributions for different candidate values have different total weights for reasons unrelated to the evidence, multiplying their weights can distort the posterior. Metropolis can sample from a non-normalized target; it does not repair an incorrectly constructed target. section 16.3
