Skip to main content

Come analyze HEASARC, IRSA, and MAST data in the cloud! The Fornax Initiative is now welcoming all interested beta users.

Xspec Home Page

next up previous contents
Next: Samplers Up: Priors Previous: Practical advice   Contents

Joint priors: the smoothness regularizer

All the priors above act on one parameter at a time. Some inference problems are instead under-determined in a way that only a prior coupling several parameters can fix. The canonical X-ray example is spatially-resolved spectroscopy with the projct model: the de-projection that recovers the emissivity of each 3-D shell from the projected annular spectra is an ill-conditioned inverse, so the per-shell normalizations come out noisy and strongly anti-correlated, and finer radial binning makes it worse rather than better. No per-parameter prior helps, because the problem is not that any single shell is unconstrained — it is that neighboring shells trade emission back and forth.

The bayes smooth command supplies the standard remedy, a Tikhonov smoothness prior. It couples a chosen model parameter (the normalization, say, or the temperature) across the datagroups and penalizes shell-to-shell roughness,

$\displaystyle \ln p(\bm{\theta}) = -\tfrac{1}{2}\,\lambda\,\lVert L\,\bm{q}\rVert^{2}, \qquad q_k = p_k \ $   or$\displaystyle \ \ln p_k,$ (3.2)

where $L$ is a spacing-aware first- or second-difference operator built from the shell radii (penalizing the gradient or the curvature of the profile respectively), and $\lambda$ sets the strength. This is exactly a zero-mean Gaussian prior on the shell parameters with precision $\lambda L^{\mathsf T}L$ — a correlated prior, the first whose density does not factorize over the parameters. The radii are taken automatically from each datagroup's XFLT major keyword, so for a projct fit the whole prior is set with a single line: 
bayes smooth 2 order 2 space log lambda 5
The strength can be fixed, or selected automatically at the corner of the misfit-versus-roughness L-curve (lcurve) or by the discrepancy principle (discrep); the latter raises the smoothing until the data misfit equals the number of degrees of freedom.

Because the prior enters through the same $\ln p(\bm{\theta})$ that the per-parameter priors use, it regularizes everything downstream with no further work: the Levenberg–Marquardt fit and its error bars, and all four samplers below. Under nest the prior of Eq. 3.2 needs one extra step, because certain profile shapes leave the penalty unchanged: adding the same constant to every shell does not alter the differences a first-difference $L$ measures, and a constant or a linear-in-radius ramp does not alter the curvatures a second-difference $L$ measures. For those components $\lVert L\,\bm{q}\rVert = 0$ (they span the null space of $L$, and make the precision matrix $\lambda L^{\mathsf T}L$ rank-deficient), so the prior expresses no preference about them at all. The prior is therefore sampled through a spectral transform: a Gaussian draw along each penalized eigen-direction, and a flat draw, within the parameter limits, along the unconstrained directions — the overall level and, for second differences, the overall slope. The posterior then shows the adjacent shells smoothed together while the overall amplitude of the profile stays free.

The smoothness prior is the first concrete case of a more general multi-parameter prior facility; the per-parameter priors of the previous section are simply its one-parameter limit. Settings round-trip through save/@file and are listed by show bayes (PyXspec Fit.bayesPriors).

next up previous contents
Next: Samplers Up: Priors Previous: Practical advice   Contents