Posterior Distance Distributions (ligo.skymap.distance)

Distance distribution functions from [1], [2], [3].

References

1

Singer, Chen, & Holz, 2016. “Going the Distance: Mapping Host Galaxies of LIGO and Virgo Sources in Three Dimensions Using Local Cosmography and Targeted Follow-up.” ApJL, 829, L15 <https://doi.org/10.3847/2041-8205/829/1/L15>.

2

Singer, Chen, & Holz, 2016. “Supplement: ‘Going the Distance: Mapping Host Galaxies of LIGO and Virgo Sources in Three Dimensions Using Local Cosmography and Targeted Follow-up’ (2016, ApJL, 829, L15).” ApJS, 226, 10 <https://doi.org/10.3847/0067-0049/226/1/10>.

3

https://asd.gsfc.nasa.gov/Leo.Singer/going-the-distance

Conditional Distributions

ligo.skymap.distance.conditional_pdf(x1, x2, x3, x4, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

Conditional distance probability density function (ansatz).

Parameters
rnumpy.ndarray

Distance (Mpc)

distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

distnormnumpy.ndarray

Distance normalization factor (Mpc^-2)

Returns
pdfnumpy.ndarray

Conditional probability density according to ansatz.

ligo.skymap.distance.conditional_cdf(x1, x2, x3, x4, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

Cumulative conditional distribution of distance (ansatz).

Parameters
rnumpy.ndarray

Distance (Mpc)

distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

distnormnumpy.ndarray

Distance normalization factor (Mpc^-2)

Returns
pdfnumpy.ndarray

Conditional probability density according to ansatz.

Examples

Test against numerical integral of pdf.

>>> import scipy.integrate
>>> distmu = 10.0
>>> distsigma = 5.0
>>> distnorm = 1.0
>>> r = 8.0
>>> expected, _ = scipy.integrate.quad(
...     conditional_pdf, 0, r,
...     (distmu, distsigma, distnorm))
>>> result = conditional_cdf(
...     r, distmu, distsigma, distnorm)
>>> np.testing.assert_almost_equal(expected, result)
ligo.skymap.distance.conditional_ppf(x1, x2, x3, x4, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

Point percent function (inverse cdf) of distribution of distance (ansatz).

Parameters
pnumpy.ndarray

The cumulative distribution function

distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

distnormnumpy.ndarray

Distance normalization factor (Mpc^-2)

Returns
rnumpy.ndarray

Distance at which the cdf is equal to p.

Examples

Test against numerical estimate.

>>> import scipy.optimize
>>> distmu = 10.0
>>> distsigma = 5.0
>>> distnorm = 1.0
>>> p = 0.16  # "one-sigma" lower limit
>>> expected_r16 = scipy.optimize.brentq(
... lambda r: conditional_cdf(r, distmu, distsigma, distnorm) - p, 0.0, 100.0)
>>> r16 = conditional_ppf(p, distmu, distsigma, distnorm)
>>> np.testing.assert_almost_equal(expected_r16, r16)

Moments

ligo.skymap.distance.moments_to_parameters(x1, x2, [out1, out2, out3, ]/, [out=(None, None, None), ]*, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

Convert ansatz moments to parameters. This function is the inverse of parameters_to_moments.

Parameters
distmeannumpy.ndarray

Conditional mean of distance (Mpc)

diststdnumpy.ndarray

Conditional standard deviation of distance (Mpc)

Returns
distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

distnormnumpy.ndarray

Distance normalization factor (Mpc^-2)

ligo.skymap.distance.parameters_to_moments(x1, x2, [out1, out2, out3, ]/, [out=(None, None, None), ]*, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

Convert ansatz parameters to moments. This function is the inverse of moments_to_parameters.

Parameters
distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

Returns
distmeannumpy.ndarray

Conditional mean of distance (Mpc)

diststdnumpy.ndarray

Conditional standard deviation of distance (Mpc)

distnormnumpy.ndarray

Distance normalization factor (Mpc^-2)

Examples

For mu=0, sigma=1, the ansatz is a chi distribution with 3 degrees of freedom, and the moments have simple expressions.

>>> mean, std, norm = parameters_to_moments(0, 1)
>>> expected_mean = 2 * np.sqrt(2 / np.pi)
>>> expected_std = np.sqrt(3 - expected_mean**2)
>>> expected_norm = 2.0
>>> np.testing.assert_allclose(mean, expected_mean)
>>> np.testing.assert_allclose(std, expected_std)
>>> np.testing.assert_allclose(norm, expected_norm)

Check that the moments scale as expected when we vary sigma.

>>> sigma = np.logspace(-8, 8)
>>> mean, std, norm = parameters_to_moments(0, sigma)
>>> np.testing.assert_allclose(mean, expected_mean * sigma)
>>> np.testing.assert_allclose(std, expected_std * sigma)
>>> np.testing.assert_allclose(norm, expected_norm / sigma**2)

Check some more arbitrary values using numerical quadrature:

>>> import scipy.integrate
>>> sigma = 1.0
>>> for mu in np.linspace(-10, 10):
...     mean, std, norm = parameters_to_moments(mu, sigma)
...     moments = np.empty(3)
...     for k in range(3):
...         moments[k], _ = scipy.integrate.quad(
...             lambda r: r**k * conditional_pdf(r, mu, sigma, 1.0),
...             0, np.inf)
...     expected_norm = 1 / moments[0]
...     expected_mean, r2 = moments[1:] * expected_norm
...     expected_std = np.sqrt(r2 - np.square(expected_mean))
...     np.testing.assert_approx_equal(mean, expected_mean, 5)
...     np.testing.assert_approx_equal(std, expected_std, 5)
...     np.testing.assert_approx_equal(norm, expected_norm, 5)
ligo.skymap.distance.parameters_to_marginal_moments(prob, distmu, distsigma)[source]

Calculate the marginal (integrated all-sky) mean and standard deviation of distance from the ansatz parameters.

Parameters
probnumpy.ndarray

Marginal probability (pix^-2)

distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

Returns
distmeanfloat

Mean distance (Mpc)

diststdfloat

Std. deviation of distance (Mpc)

Marginal Distributions

ligo.skymap.distance.marginal_pdf(x1, x2, x3, x4, x5, /, out=None, *, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

Calculate all-sky marginal pdf (ansatz).

Parameters
rnumpy.ndarray

Distance (Mpc)

probnumpy.ndarray

Marginal probability (pix^-2)

distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

distnormnumpy.ndarray

Distance normalization factor (Mpc^-2)

Returns
pdfnumpy.ndarray

Marginal probability density according to ansatz.

Examples

>>> npix = 12
>>> prob, distmu, distsigma, distnorm = np.random.uniform(size=(4, 12))
>>> r = np.linspace(0, 1)
>>> pdf_expected = np.dot(
...     conditional_pdf(r[:, np.newaxis], distmu, distsigma, distnorm), prob)
>>> pdf = marginal_pdf(r, prob, distmu, distsigma, distnorm)
>>> np.testing.assert_allclose(pdf, pdf_expected, rtol=1e-4)
ligo.skymap.distance.marginal_cdf(x1, x2, x3, x4, x5, /, out=None, *, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

Calculate all-sky marginal cdf (ansatz).

Parameters
rnumpy.ndarray

Distance (Mpc)

probnumpy.ndarray

Marginal probability (pix^-2)

distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

distnormnumpy.ndarray

Distance normalization factor (Mpc^-2)

Returns
cdfnumpy.ndarray

Marginal cumulative probability according to ansatz.

Examples

>>> npix = 12
>>> prob, distmu, distsigma, distnorm = np.random.uniform(size=(4, 12))
>>> r = np.linspace(0, 1)
>>> cdf_expected = np.dot(
...     conditional_cdf(r[:, np.newaxis], distmu, distsigma, distnorm), prob)
>>> cdf = marginal_cdf(r, prob, distmu, distsigma, distnorm)
>>> np.testing.assert_allclose(cdf, cdf_expected, rtol=1e-4)
ligo.skymap.distance.marginal_ppf(x1, x2, x3, x4, x5, /, out=None, *, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

Point percent function (inverse cdf) of marginal distribution of distance (ansatz).

Parameters
pnumpy.ndarray

The cumulative distribution function

probnumpy.ndarray

Marginal probability (pix^-2)

distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

distnormnumpy.ndarray

Distance normalization factor (Mpc^-2)

Returns
rnumpy.ndarray

Distance at which the cdf is equal to p.

Examples

>>> from astropy.utils.misc import NumpyRNGContext
>>> npix = 12
>>> with NumpyRNGContext(0):
...     prob, distmu, distsigma, distnorm = np.random.uniform(size=(4, 12))
>>> r_expected = np.linspace(0.4, 0.7)
>>> cdf = marginal_cdf(r_expected, prob, distmu, distsigma, distnorm)
>>> r = marginal_ppf(cdf, prob, distmu, distsigma, distnorm)
>>> np.testing.assert_allclose(r, r_expected, rtol=1e-4)

Kernel Density Estimation

ligo.skymap.distance.conditional_kde(n, datasets, inverse_covariances, weights)[source]
ligo.skymap.distance.cartesian_kde_to_moments(n, datasets, inverse_covariances, weights)[source]

Calculate the marginal probability, conditional mean, and conditional standard deviation of a mixture of three-dimensional kernel density estimators (KDEs), in a given direction specified by a unit vector.

Parameters
nnumpy.ndarray

A unit vector; an array of length 3.

datasetslist of numpy.ndarray

A list 2D Numpy arrays specifying the sample points of the KDEs. The first dimension of each array is 3.

inverse_covariances: list of `numpy.ndarray`

An array of 3x3 matrices specifying the inverses of the covariance matrices of the KDEs. The list has the same length as the datasets parameter.

weightslist

A list of floating-point weights.

Returns
probfloat

The marginal probability in direction n, integrated over all distances.

meanfloat

The conditional mean in direction n.

stdfloat

The conditional standard deviation in direction n.

Examples

>>> # Some imports
>>> import scipy.stats
>>> import scipy.integrate
>>> # Construct random dataset for KDE
>>> np.random.seed(0)
>>> nclusters = 5
>>> ndata = np.random.randint(0, 1000, nclusters)
>>> covs = [np.random.uniform(0, 1, size=(3, 3)) for _ in range(nclusters)]
>>> covs = [_ + _.T + 3 * np.eye(3) for _ in covs]
>>> means = np.random.uniform(-1, 1, size=(nclusters, 3))
>>> datasets = [np.random.multivariate_normal(m, c, n).T
...     for m, c, n in zip(means, covs, ndata)]
>>> weights = ndata / float(np.sum(ndata))
>>>
>>> # Construct set of KDEs
>>> kdes = [scipy.stats.gaussian_kde(_) for _ in datasets]
>>>
>>> # Random unit vector n
>>> n = np.random.normal(size=3)
>>> n /= np.sqrt(np.sum(np.square(n)))
>>>
>>> # Analytically evaluate conditional mean and std. dev. in direction n
>>> datasets = [_.dataset for _ in kdes]
>>> inverse_covariances = [_.inv_cov for _ in kdes]
>>> result_prob, result_mean, result_std = cartesian_kde_to_moments(
...     n, datasets, inverse_covariances, weights)
>>>
>>> # Numerically integrate conditional distance moments
>>> def rkbar(k):
...     def integrand(r):
...         return r ** k * np.sum([kde(r * n) * weight
...             for kde, weight in zip(kdes, weights)])
...     integral, err = scipy.integrate.quad(integrand, 0, np.inf)
...     return integral
...
>>> r0bar = rkbar(2)
>>> r1bar = rkbar(3)
>>> r2bar = rkbar(4)
>>>
>>> # Extract conditional mean and std. dev.
>>> r1bar /= r0bar
>>> r2bar /= r0bar
>>> expected_prob = r0bar
>>> expected_mean = r1bar
>>> expected_std = np.sqrt(r2bar - np.square(r1bar))
>>>
>>> # Check that the two methods give almost the same result
>>> np.testing.assert_almost_equal(result_prob, expected_prob)
>>> np.testing.assert_almost_equal(result_mean, expected_mean)
>>> np.testing.assert_almost_equal(result_std, expected_std)
>>>
>>> # Check that KDE is normalized over unit sphere.
>>> nside = 32
>>> npix = hp.nside2npix(nside)
>>> prob, _, _ = np.transpose([cartesian_kde_to_moments(
...     np.asarray(hp.pix2vec(nside, ipix)),
...     datasets, inverse_covariances, weights)
...     for ipix in range(npix)])
>>> result_integral = prob.sum() * hp.nside2pixarea(nside)
>>> np.testing.assert_almost_equal(result_integral, 1.0, decimal=4)

Miscellaneous

ligo.skymap.distance.ud_grade(prob, distmu, distsigma, *args, **kwargs)[source]

Upsample or downsample a distance-resolved sky map.

Parameters
probnumpy.ndarray

Marginal probability (pix^-2)

distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

*args, **kwargs :

Additional arguments to healpy.ud_grade (e.g., nside, order_in, order_out).

Returns
probnumpy.ndarray

Resampled marginal probability (pix^-2)

distmunumpy.ndarray

Resampled distance location parameter (Mpc)

distsigmanumpy.ndarray

Resampled distance scale parameter (Mpc)

distnormnumpy.ndarray

Resampled distance normalization factor (Mpc^-2)

ligo.skymap.distance.volume_render(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, /, out=None, *, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

Perform volumetric rendering of a 3D sky map.

Parameters
xnumpy.ndarray

X-coordinate in rendered image

ynumpy.ndarray

Y-coordinate in rendered image

max_distancefloat

Limit of integration from -max_distance to +max_distance

axis0int

Index of axis to assign to x-coordinate

axis1int

Index of axis to assign to y-coordinate

Rnumpy.ndarray

Rotation matrix as provided by principal_axes

nestbool

HEALPix ordering scheme

probnumpy.ndarray

Marginal probability (pix^-2)

distmunumpy.ndarray

Distance location parameter (Mpc)

distsigmanumpy.ndarray

Distance scale parameter (Mpc)

distnormnumpy.ndarray

Distance normalization factor (Mpc^-2)

Returns
imagenumpy.ndarray

Rendered image

Examples

Test volume rendering of a normal unit sphere… First, set up the 3D sky map.

>>> nside = 32
>>> npix = hp.nside2npix(nside)
>>> prob = np.ones(npix) / npix
>>> distmu = np.zeros(npix)
>>> distsigma = np.ones(npix)
>>> distnorm = np.ones(npix) * 2.0

The conditional distance distribution should be a chi distribution with 3 degrees of freedom.

>>> from scipy.stats import norm, chi
>>> r = np.linspace(0, 10.0)
>>> actual = conditional_pdf(r, distmu[0], distsigma[0], distnorm[0])
>>> expected = chi(3).pdf(r)
>>> np.testing.assert_almost_equal(expected, actual)

Next, run the volume renderer.

>>> dmax = 4.0
>>> n = 64
>>> s = np.logspace(-dmax, dmax, n)
>>> x, y = np.meshgrid(s, s)
>>> R = np.eye(3)
>>> P = volume_render(x, y, dmax, 0, 1, R, False,
...                   prob, distmu, distsigma, distnorm)

Next, integrate analytically.

>>> P_expected = norm.pdf(x) * norm.pdf(y) * (norm.cdf(dmax) - norm.cdf(-dmax))

Compare the two.

>>> np.testing.assert_almost_equal(P, P_expected, decimal=4)

Check that we get the same answer if the input is in ring ordering. FIXME: this is a very weak test, because the input sky map is isotropic!

>>> P = volume_render(x, y, dmax, 0, 1, R, True,
...                   prob, distmu, distsigma, distnorm)
>>> np.testing.assert_almost_equal(P, P_expected, decimal=4)

Last, check that we don’t have a coordinate singularity at the origin.

>>> x = np.concatenate(([0], np.logspace(1 - n, 0, n) * dmax))
>>> y = 0.0
>>> P = volume_render(x, y, dmax, 0, 1, R, False,
...                   prob, distmu, distsigma, distnorm)
>>> P_expected = norm.pdf(x) * norm.pdf(y) * (norm.cdf(dmax) - norm.cdf(-dmax))
>>> np.testing.assert_allclose(P, P_expected, rtol=1e-4)
ligo.skymap.distance.principal_axes(prob, distmu, distsigma, nest=False)[source]