Select covariates with reversible-jump MCMC

Scientific question

Which predictors belong in the regression?

The checked-in example has sixteen observations and three candidate predictors: dose, genotype, and a correlated proxy_marker. The response was generated from the first two, while the proxy creates a realistic competing-locus model that exercises swap moves. A model identifier is a bit mask, so model 3 contains dose and genotype while model 0 contains only the common intercept.

For subset M, the model is:

y[i] ~ Normal(alpha + sum(beta[j] * x[i,j], j in M), 0.75)
alpha ~ Normal(0, 5)
beta[j] ~ Normal(0, 2)
Pr(j in M) = 0.30 independently

Correct target

Keep every model-dependent constant

SubsetSelectionTarget takes a Java callback that receives common parameters, sorted active candidate indices, and active coefficients. Unlike an ordinary within-model target, it must return the complete normalized log joint. Constants that differ between models directly affect posterior model probabilities.

SubsetSelectionTarget target = new SubsetSelectionTarget(
    new String[] {"alpha"},
    new String[] {"dose", "genotype", "proxy_marker"},
    (common, active, coefficients) -> {
        double result = normalLogDensity(common[0], 0, 5);
        result += active.length * Math.log(0.30)
            + (3 - active.length) * Math.log(0.70);
        for (double beta : coefficients)
            result += normalLogDensity(beta, 0, 2);
        // Add every normalized observation density here.
        return result;
    });

Detailed balance

Add, drop, swap, then update within the model

SubsetSelectionRj schedules three reversible structure moves. Add chooses uniformly among inactive predictors and samples a new coefficient from its declared birth density. Drop chooses uniformly among active predictors and evaluates the matching reverse birth density. Swap chooses one of each and preserves dimension. The sampler also includes the probability of selecting each move among those available at the current and proposed boundaries.

GaussianRjBirthProposal birth =
    new GaussianRjBirthProposal(0, 2);
ReversibleJumpSampler sampler =
    SubsetSelectionRj.sampler(birth, 0.20, 0.30);

The general engine evaluates:

log alpha = logJoint(proposed) - logJoint(current)
          + log q(reverse) - log q(forward)
          + log p(select reverse) - log p(select forward)
          + logAbsJacobian

CoordinateInsertionTransformation supplies the unit-Jacobian birth/death map and DimensionMatchingValidator checks total dimensions, inverse round trips, and reciprocal Jacobians. Custom ReversibleJumpMove implementations can supply nonzero Jacobians and arbitrary Java model structures.

Four independent chains

Warm up proposals, freeze them, and retain ragged draws

ReversibleJumpSamplingOptions options =
    ReversibleJumpSamplingOptions.builder()
        .warmupIterations(1500)
        .sampleIterations(4000)
        .targetJumpAcceptance(0.25)
        .adaptMoveWeights(true)
        .build();

ReversibleJumpResult chain = sampler.sample(
    target,
    target.state(0, new double[] {0}, new double[0]),
    options,
    new MersenneTwister(20260829));

Move weights, per-model random-walk scales, and adaptive birth proposals may change only during discarded warmup. They are frozen before retained sampling. Use a separate sampler instance and explicit random engine for each independent chain.

Compile and run the complete file

The classpath wildcard stays valid across JDistlib releases.

./gradlew compileDocumentationExamples
java -cp "build/libs/*:build/documentation-examples" \
  examples.WorkedReversibleJumpSelectionExample

# Windows
gradlew.bat compileDocumentationExamples
java -cp "build\libs\*;build\documentation-examples" examples.WorkedReversibleJumpSelectionExample

An optional first argument selects the output directory. The default is build/example-output/rjmcmc.

Review the trans-dimensional chain

Model occupancy alone is not enough

ReversibleJumpDiagnostics.analyze combines chains and reports posterior model and candidate inclusion probabilities with ESS/MCSE/R-hat, retained transition counts, round trips between the two most visited models, add/drop/swap acceptance and invalid proposals, and parameter summaries conditional on presence.

ReversibleJumpDiagnosticReport report =
    ReversibleJumpDiagnostics.analyze(target, chains);

double[] inclusion = report.inclusionProbabilities();
double[] inclusionEss = report.inclusionEffectiveSampleSizes();
double[] inclusionRhat = report.inclusionRHats();

Investigate a candidate with low inclusion ESS or R-hat above 1.01, a move direction with near-zero acceptance, no round trips, models seen in only one chain, or active coefficients with very few conditional draws. The example writes rjmcmc-diagnostics.json using schema jdistlib.rjmcmc-diagnostics/1 and rjmcmc-chain-1.csv in tidy ragged form.

Reproducibility

Resume the exact frozen process

The portable checkpoint includes the model identifier, variable-length parameter vector, log joint, completed iteration, cloned random stream, frozen move weights, and every component's adaptation state. Model and option fingerprints reject incompatible restores.

ReversibleJumpCheckpointIO.write(path, chain.checkpoint(),
    "worked-rj-linear-v1", "warmup-1500-samples-4000-v1");

PortableReversibleJumpCheckpoint restored =
    ReversibleJumpCheckpointIO.read(path,
        "worked-rj-linear-v1", "warmup-1500-samples-4000-v1");

ReversibleJumpResult continued = newSampler.resume(
    target, restored.checkpoint(),
    ReversibleJumpSamplingOptions.builder()
        .warmupIterations(0).sampleIterations(250).build());