-
Notifications
You must be signed in to change notification settings - Fork 11
Adaptive MCMC sampling in PROfit
This is a quick breakdown of PROfit framework's implementation to Adaptive Markov Chain Monte Carlo (MCMC) sampling. This page provides a brief guide to the adaptive proposal mechanism used in PROfit, the motivation behind it, and technical details to help users and developers understand, debug, and hopefully extend the codebase.
For the most part, Heikki Haario, Eero Saksman, Johanna Tamminen (2001). An adaptive Metropolis algorithm is a good reference.
We are actually not (yet) using MCMC for the actual fit itself, but rather as a robust way to estimate error bars post-fit. We have two types of systematics, covariance and spline (pull-term/nuisance) and pre-fit the errors due to covariance are simply the sqrt of diagonal and we generally treat the splines as uncorrelated gaussians so can easily show their individual and combined effects on the spectra. Post-fit, however, is different as now many of the splines may be correlated and the data itself constrains many of the splines stronger than the priors.
This is where MCMC comes into play. Its a technique for generating correctly generated samples from complicated, high-dimensional probability distributions. In our case we want to estimate the distribution of model parameters given observed data, i.e we want to have the correlated distributions of all splines and physics params after the fit has taken place.
MCMC solves this by:
- Constructing a chain of parameter samples, where each one depends on the prior sample
- Ensuring that, after an initial burn-in period, these samples are distributed according to the target (posterior) distribution.
- Allowing us to estimate means, correlations, and the direct effect on final observable using these samples.
Why adaptive? MCMC
Standard MCMC (Metropolis-Hastings) requires a user-tuned proposal distribution, which is often inefficient or needs time-consuming tuning. Adaptive MCMC learns the optimal proposal width, shape, and orientation during sampling, making it easier to use and more efficient. We do note that we don't use MCMC for the final fits, results or chi2 itself just for drawing visual error bars.
PROfit uses an adaptive multivariate Gaussian proposal, where the proposal covariance and global scale are updated on the fly.
- Start with a diagonal identity covariance .
-
At each accepted proposal:
- Update the empirical mean and covariance using running statistics. (Welford's online algorithm)
- Covariance matrix "shrinks" if the chain is too exploratory, or "grows" if it is stuck.
-
Every N steps (the adaptation window, usually 1000):
- Calculate the acceptance rate for recent proposals.
- Adjust the global scale to target a desired acceptance rate (23.4% is optimal for high-dimensional problems
- If acceptance rate > 23.4%, increase the scale
- If acceptance rate < 23.4%, decrease the scale
To draw a random sample from a multivariate normal with covariance A, you need to transform standard normal vectors (z) into the target space. The transformation is:
sample = mean + L * z
where L satisfies A = LL^T (the matrix square root).
Any way to generate an L such that the above is true works, and there are two main methods supported in PROfit. Both using Eigen inbuilt tools.
-
A = L D Lᵀis the LDLT decomposition. - We construct
LasPᵗ * Lp * D^{1/2}(wherePis permutation,Lpis lower triangular,Dis diagonal). - Fast and numerically stable if the matrix is positive-definite.
-
A = U S Uᵀis the Singular Value Decomposition. - Transform with
L = U * sqrt(S). - More robust to singularities or small eigenvalues, but slower
Even with adaptation, MCMC samples are correlated, each sample depends on the last. Autocorrelation measures how many steps you must skip before samples become effectively independent. We can study this by comparing events to events in the chain lagging behind teh current event
- Short autocorrelation: Efficient exploration, more "effective" samples per run.
- Long autocorrelation: Chain is "sticky", this is overall bad. you need many more total samples to get the same statistical precision.
It's common to plot an autocorrelation factor for each parameter. PROfit can output info for these, but the plots are made in the jupyter-notebook contained in the other folder. Here is an example showing good short autocorrelation that is pretty much gone from all samples by 50 lag.
You ideally want autocorrelation times (the number of steps to near-independence) to be as small as possible—this validates that the adaptive proposal is tuned well.