[audit 2026-07-09] overly Symmetrised proposals in MCMC

[odunbar]

Issue

Not all proposals are symmetric, for example the pCN proposal here:

CalibrateEmulateSample.jl/src/MarkovChainMonteCarlo.jl

Lines 343 to 354 in 26d7455

	function AdvancedMH.propose(
	rng::Random.AbstractRNG,
	sampler::pCNMetropolisHastings,
	model::AdvancedMH.DensityModel,
	current_state::MCMCState;
	stepsize::FT = 1.0,
	) where {FT <: AbstractFloat}
	# Use prescription in Beskos et al (2017) "Geometric MCMC for infinite-dimensional
	# inverse problems." for relating ρ to Euler stepsize:
	ρ = (1 - stepsize / 4) / (1 + stepsize / 4)
	return ρ * current_state.params .+ sqrt(1 - ρ^2) * rand(rng, sampler.proposal)
	end

has a transition density that is not-symmetric i.e. q(u,u') != q(u',u). However the logratio_proposal_density that should do "q(u,u') - q(u',u)" directly when non-symmetric. Instead, all methods currently default to a _get_proposal call that was implemented for RWM, and a transition that assumes a symmetric update.

CalibrateEmulateSample.jl/src/MarkovChainMonteCarlo.jl

Lines 164 to 176 in 26d7455

	AdvancedMH.logratio_proposal_density(
	sampler::pCNMetropolisHastings,
	transition_prev::AdvancedMH.AbstractTransition,
	candidate,
	) = AdvancedMH.logratio_proposal_density(sampler.proposal, transition_prev.params, candidate)
	function MetropolisHastingsSampler(
	::pCNMHSampling{T},
	encoded_prior::ParameterDistribution,
	encoder_schedule::VV,
	) where {T <: AutodiffProtocol, VV <: AbstractVector}
	proposal = _get_proposal(encoded_prior, encoder_schedule)
	return pCNMetropolisHastings{typeof(proposal), T}(proposal)
	end

This applies to pCN and Barker and the effect e.g. for pCN is that there is a doubled prior contribution.

Possible solution:

I think the bug was from some copied boilerplate (due to a misconception of the messy AdvancedMH interface we had to use). It is a relatively straightforward change. One simply writes the transition density for each code, (not based on the same _get_proposal) and then fully write out the ratio as "reverse - forward". In the symmetric case this cancels but otherwise it will not.

# Each sampler implements this
function log_transition_density(sampler::MHS, θ_from, θ_to) where {MHS <: AdvancedMH.MHSampler}      
    ...
end   

# Forward and reverse called by all samplers
function AdvancedMH.logratio_proposal_density(sampler::MHS, transition_prev, candidate) where {MHS   <: AdvancedMH.MHSampler}
       log_q_forward = log_transition_density(sampler, transition_prev.params, candidate)  #    q(cand| prev)       
       log_q_reverse = log_transition_density(sampler, candidate, transition_prev.params)  # q(prev | cand)       
       return log_q_reverse - log_q_forward  
end

I think furthermore, for many parametric types, we could implement the Markov Kernel directly, and then the transition becomes the logpdf(kernel) while the propose() method just becomes a "rand!(kernel)", reducing the multiplicity. However this might not apply to more generic kernels such as Barker...