Breaking the multimodalities

[vlas-sokolov]

One issue that often pops up in the sampled spectra is that the posterior is symmetrically multimodal, indicating that the velocity some components are "mirrored" on each other.

This is expected, and here's a reason of why this happens. As we don't have the prior information on where exactly the split in velocity occurs, the easiest way to set a velocity prior is to pass the same prior for centroid positions of all the components. When the live points are generated, they do so randomly within the constrained prior volume, without the knowledge of what spectral components come first or last. So when the posterior is sampled, we would often get a mess where different components are mixed together in

Why it's not such a big of a concern? The "best-fit" parameters don't care whether or not we broke the multimodalities - a point estimate of maximum likelihood would not change. On top of that, the evidences computed integrate over the whole parameter space anyway, and should also be ignorant of the component separation. This is reflected in the smoothness of the evidence and Bayes factor maps.

Why it is crucial to fix anyway? Uncertainty analyses become one hell of a problem when your posterior is multimodal. Also, MLE component maps would look way more coherent then.

[vlas-sokolov]

A neat solution for this would be able to reparametrize velocity centroids of ModelContainer class. The new set of parameters could be, for example, the average velocity of all the components, and velocity jump to the next component, then the jump after it, and so on. This way the velocities in the posterior are bound to be ordered.

This is actually quite a neat trick - we can set minimal priors on the velocity separations! Ongoing work in the ordered-vlsr branch.

[vlas-sokolov]

The time for coordinate transformation should be negligible given that the transformation matrix was precomputed:

In [42]: T = get_xoff_transform(5)
In [43]: Tinv = np.linalg.inv(T)
In [44]: %timeit np.dot(Tinv, [44, 0.2, 0.3, 1, 0.1])
3.08 µs ± 77.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

[vlas-sokolov]

Reproduced the issue on a synthetic Gaussian spectrum.

Summary of the toy model setup:

1. Generate a mix of two bell curve-ish spectral lines. Profiles overlap significantly, and have different line widths and amplitudes. Sprinkle white noise to taste (peak SNR ~ 15).
   

2. Then, consider a set of nested models with 1-, 2-, and 3 line-of-sight components. Due to our ignorance of the parameter sets, the priors on the free parameter triples {amplitude, centroid, line width} x N_components are identical.

3. Sample the npeaks=1 model posterior and get the Bayesian evidence (Z1) for it.

4. Similarly, sample the npeaks=2 model posterior and obtain the Bayesian evidence Z2.
   

And there we have it, the bimodal behavior resulting from components not knowing of each other's existence - local minima are sampled in parallel and thus mixed together in the posterior.

[vlas-sokolov]

Great Scott! The very first try of coordinate transformation yields clear, singly-peaked posteriors:

The wrapped marginals couldn't look any better too!

[vlas-sokolov]

Bayesian evidence for a two-component model doesn't differ much between the two , although more testing wouldn't hurt.

Line centroid specification	Bayesian evidence
Normal coordinates	ln(Z₂) = -508.8 ± 0.7
Transformed coordinates	ln(Z₂) = -507.1 ± 0.7

(both computed for n_live_points = 4000 and sampling_efficiency = 0.3)

Marginalized posterior distributions for both cases:

[vlas-sokolov]

An important note - when we use the reparametrized xoff_transform decorator, it operates on a view of the cube variable passed from within MultiNest itself. As xoff_transform converts {x_mean, d_x1, d_x2, ...} into {x1, x1, x3, ...} coordinates, it modifies cube elements in place so that the MultiNest output is actually in the "normal" coordinate system!

[vlas-sokolov]

Here's a comparison of the two posteriors for the real VLA ammonia data, for n_live_points = 4000 and sampling_efficiency = 0.3. This time the contours are 1σ, 2σ, and 3σ:

[vlas-sokolov]

Despite the feature branch working better than I've expected, a few things still need to be done:

• fix the relative imports of the tests subpackage in the Python 3 version
• add a function that would wrap the priors in the xoff space

[vlas-sokolov]

For now the priors on centroid separations can be input by hand, time to merge the branch.