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. doi: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. doi:10.3847/0067-0049/226/1/10

Conditional Distributions

conditional_pdf(x1, x2, x3, x4, /[, out, ...])

Conditional distance probability density function (ansatz).

conditional_cdf(x1, x2, x3, x4, /[, out, ...])

Cumulative conditional distribution of distance (ansatz).

conditional_ppf(x1, x2, x3, x4, /[, out, ...])

Point percent function (inverse cdf) of distribution of distance (ansatz).
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(result, expected)
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(r16, expected_r16)

Moments

moments_to_parameters(x1, x2[, out1, out2, ...)

Convert ansatz moments to parameters.

parameters_to_moments(x1, x2[, out1, out2, ...)

Convert ansatz parameters to moments.

parameters_to_marginal_moments(prob, distmu, ...)

Calculate the marginal (integrated all-sky) mean and standard deviation of distance from the ansatz parameters.
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] [edit on github]

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

marginal_pdf(x1, x2, x3, x4, x5, /[, out, ...])

Calculate all-sky marginal pdf (ansatz).

marginal_cdf(x1, x2, x3, x4, x5, /[, out, ...])

Calculate all-sky marginal cdf (ansatz).

marginal_ppf(x1, x2, x3, x4, x5, /[, out, ...])

Point percent function (inverse cdf) of marginal distribution of distance (ansatz).
ligo.skymap.distance.marginal_pdf(x1, x2, x3, x4, x5, /, out=None, *, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj, axes, axis])

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, axes, axis])

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, axes, axis])

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

conditional_kde(n, datasets, ...)

cartesian_kde_to_moments(n, datasets, ...)

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.
ligo.skymap.distance.conditional_kde(n, datasets, inverse_covariances, weights)[source] [edit on github]
ligo.skymap.distance.cartesian_kde_to_moments(n, datasets, inverse_covariances, weights)[source] [edit on github]

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 = ah.nside_to_npix(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

ud_grade(prob, distmu, distsigma, *args, ...)

Upsample or downsample a distance-resolved sky map.

volume_render(x1, x2, x3, x4, x5, x6, x7, ...)

Perform volumetric rendering of a 3D sky map.

principal_axes(prob, distmu, distsigma[, nest])

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

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, axes, axis])

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 = ah.nside_to_npix(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(actual, expected)

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] [edit on github]