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