MultivariateGaussianDist(): use "standard" multivariate normal distribution internally

A student of mine (Hank Hua) is working on a CW project using CWInPy. We're using a MultivariateGaussianDist() to set priors on the CW frequency and frequency derivative (spindown) parameters. Because the typical scales of these parameters are very different (frequency ~ 100, spindown ~1 / year), the covariance matrix can appear close to singular, which leads to NumPy errors, e.g.

numpy.linalg.LinAlgError: When `allow_singular is False`, the input matrix must be symmetric positive definite.

This MR modifies MultivariateGaussianDist() to internally use a "standard" multivariate normal distribution. In analogy to a one-dimensional standard normal distribution, where

\text{PDF}( x \,|\, \mu=\mu, \sigma=\sigma ) = \frac{1}{\sigma} \text{PDF}[ (x - \mu)/\sigma \,|\, \mu=0, \sigma=1 ]

the multidimensional analogue is

\text{PDF}( \vec x \,|\, \vec\mu=\vec\mu, \mathbf\Sigma=\mathbf\Sigma ) = \frac{1}{\prod_i \sigma_i} \text{PDF}[ (\vec x - \vec\mu)/\vec\sigma \,|\, \vec\mu=\vec 0, \mathbf\Sigma=\mathbf C ]

where

\mathbf C

is the correlation matrix. I've attached a PDF of a Mathematica notebook which confirms the derivation.

The code modifications are:

• Add a new array to MultivariateGaussianDist(): self.logprodsigmas to compute the log-product of self.sigmas for the PDF normalisation
• In add_mode() add to self.mvn a "standard" multivariate normal distribution as above with mus=0 and cov=corrcoefs
• In _ln_prob() call self.mvn[i].logpdf() with normalised sample z = (x - mus) / sigmas and subtract normalisation in self.logprodsigmas
• Add tests in test/core/prior/joint_test.py for _ln_prob() for MultivariateGaussianDist() with a) order-unity-scale parameters b) CW frequency/spindown parameter scalings

@karl-wette and Hank, thanks for this. The change makes sense to me. It'll probably be later in the week before I can check this properly, but I don't see any immediate problems.

Thanks @matthew-pitkin ! No time urgency on our side, just wanted to make sure this change made its way back into Bilby for future users.

@karl-wette I've run the multivariate_gaussian_prior.py example using this MR and get consistent results with the original code (the first plot is the original result and the second is with your patch), so I'm happy that it is correct.

approved this merge request

changed milestone to %1.2.1

approved this merge request

merged

mentioned in commit d53e5c4e