|
|
Here we review the changes made for pesummary v0.6.0.
|
|
|
|
|
|
[[_TOC_]]
|
|
|
|
|
|
## Comparison to `cbcBayesPostProc`
|
|
|
|
|
|
The result of running the `summaryreview` script can be found here:
|
|
|
|
|
|
#### Review comments
|
|
|
|
|
|
## Jensen-Shannon divergence calculation
|
|
|
|
|
|
The `pesummary.utils.utils.jensen_shannon_divergence` allows for both any kdes to be used. We will compare `gaussian_kde` and `Bounded_1d_kde` of these to the [`scipy.spatial.distance.jensenshannon`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.jensenshannon.html) to confirm the values. We also check the values against the code used in the [GWTC-1 paper](https://git.ligo.org/publications/O2/cbc-catalog/-/blob/master/event_plots/event_plots_utility.py#L231-260)
|
|
|
|
|
|
### Gaussian kde
|
|
|
|
|
|
```python
|
|
|
from pesummary.utils import utils
|
|
|
from scipy.spatial.distance import jensenshannon
|
|
|
from scipy import stats
|
|
|
import numpy as np
|
|
|
|
|
|
np.random.seed(123456789)
|
|
|
|
|
|
mus = np.random.uniform(10, 100, 100)
|
|
|
sigmas = np.random.uniform(0.1, 0.7, 100)
|
|
|
|
|
|
diff = []
|
|
|
for mu, sigma in zip(mus, sigmas):
|
|
|
samples = [
|
|
|
np.random.uniform(mu, sigma, 100),
|
|
|
np.random.uniform(mu, sigma, 100)
|
|
|
]
|
|
|
x = np.linspace(np.min(samples), np.max(samples), 100)
|
|
|
kde = [stats.gaussian_kde(i)(x) for i in samples]
|
|
|
_scipy = jensenshannon(*kde)**2
|
|
|
_pesummary = utils.jensen_shannon_divergence(samples, decimal=9)
|
|
|
np.testing.assert_almost_equal(_scipy, _pesummary)
|
|
|
diff.append(np.abs(_scipy - _pesummary))
|
|
|
|
|
|
print(np.max(diff))
|
|
|
```
|
|
|
|
|
|
with output:
|
|
|
|
|
|
```bash
|
|
|
4.988973179888279e-10
|
|
|
```
|
|
|
|
|
|
### Bounded 1d kde
|
|
|
|
|
|
```python
|
|
|
from pesummary.utils import utils
|
|
|
from scipy.spatial.distance import jensenshannon
|
|
|
from pesummary.core.plots.bounded_1d_kde import Bounded_1d_kde
|
|
|
from scipy import stats
|
|
|
import numpy as np
|
|
|
|
|
|
np.random.seed(123456789)
|
|
|
|
|
|
mus = np.random.uniform(0.75, 0.5, 100)
|
|
|
sigmas = np.random.uniform(0.1, 0.4, 100)
|
|
|
|
|
|
diff = []
|
|
|
for mu, sigma in zip(mus, sigmas):
|
|
|
samples = [
|
|
|
np.random.uniform(mu, sigma, 100),
|
|
|
np.random.uniform(mu, sigma, 100)
|
|
|
]
|
|
|
x = np.linspace(np.min(samples), np.max(samples), 100)
|
|
|
kde = [Bounded_1d_kde(i, xlow=0.0, xhigh=1.0)(x) for i in samples]
|
|
|
_scipy = jensenshannon(*kde)**2
|
|
|
_pesummary = utils.jensen_shannon_divergence(samples, decimal=9, kde=Bounded_1d_kde, xlow=0.0, xhigh=1.0)
|
|
|
np.testing.assert_almost_equal(_scipy, _pesummary)
|
|
|
diff.append(np.abs(_scipy - _pesummary))
|
|
|
|
|
|
print(np.max(diff))
|
|
|
```
|
|
|
|
|
|
with output:
|
|
|
|
|
|
```bash
|
|
|
4.988971141588194e-10
|
|
|
```
|
|
|
|
|
|
### Comparison to GWTC-1
|
|
|
|
|
|
```python
|
|
|
from pesummary.core.plots.bounded_1d_kde import Bounded_1d_kde
|
|
|
from pesummary.utils import utils
|
|
|
from scipy.spatial.distance import jensenshannon
|
|
|
import scipy as sp
|
|
|
from scipy.stats import entropy
|
|
|
import numpy as np
|
|
|
|
|
|
np.random.seed(123456789)
|
|
|
|
|
|
def makeKde(data, Nbins, xMin, xMax, bounded=False, bw_method=None):
|
|
|
"""Compute KDE of `data` normalised across interval [`xMin`, `xMax`]
|
|
|
using Nbins bins.
|
|
|
If `bounded` is true, use a bounded KDE method with bounds `xMin`, `xMax`.
|
|
|
"""
|
|
|
xPoints = np.linspace(xMin, xMax, Nbins)
|
|
|
if bounded:
|
|
|
kernel = Bounded_1d_kde(np.atleast_2d(data).T,
|
|
|
xlow=xMin, xhigh=xMax, bw_method=bw_method)
|
|
|
return kernel(np.atleast_2d(xPoints).T)
|
|
|
else:
|
|
|
kernel = sp.stats.gaussian_kde(data)
|
|
|
return kernel.evaluate(xPoints) / kernel.integrate_box_1d(xMin, xMax)
|
|
|
|
|
|
def JSD(pos1, pos2, Q_min=None, Q_max=None, binNum=1000, bounded=False,
|
|
|
bw_method=None, nSamples=None, seed=1234, base=2):
|
|
|
"""Compute Jensen-Shannon divergence between marginal
|
|
|
posteriors `pos1` and `pos2`.
|
|
|
Use `binNum` bins for evaluating the KDEs.
|
|
|
Can specify bounds and whether to use a bounded KDE.
|
|
|
Use `base` for the logarithm in the entropy.
|
|
|
`base`=2 will return the JSD in bits.
|
|
|
"""
|
|
|
if not bounded:
|
|
|
if not Q_min:
|
|
|
Q_min = min(np.min(pos1), np.min(pos2))
|
|
|
if not Q_max:
|
|
|
Q_max = max(np.max(pos1), np.max(pos2))
|
|
|
if nSamples:
|
|
|
np.random.seed(seed)
|
|
|
pos1_draws = np.random.choice(pos1, nSamples, replace=True)
|
|
|
pos2_draws = np.random.choice(pos2, nSamples, replace=True)
|
|
|
else:
|
|
|
pos1_draws = pos1
|
|
|
pos2_draws = pos2
|
|
|
kde1 = makeKde(pos1_draws, binNum, Q_min, Q_max, bounded=bounded,
|
|
|
bw_method=bw_method)
|
|
|
kde2 = makeKde(pos2_draws, binNum, Q_min, Q_max, bounded=bounded,
|
|
|
bw_method=bw_method)
|
|
|
kde_m = 0.5*(kde1 + kde2)
|
|
|
return 0.5 * (entropy(kde1, kde_m, base=base) + entropy(kde2, kde_m, base=base))
|
|
|
|
|
|
mus = np.random.uniform(10, 100, 100)
|
|
|
sigmas = np.random.uniform(0.1, 0.7, 100)
|
|
|
|
|
|
diff = []
|
|
|
for mu, sigma in zip(mus, sigmas):
|
|
|
samples = [
|
|
|
np.random.uniform(mu, sigma, 100),
|
|
|
np.random.uniform(mu, sigma, 100)
|
|
|
]
|
|
|
_GWTC1 = JSD(samples[0], samples[1], base=2, binNum=100)
|
|
|
_pesummary = utils.jensen_shannon_divergence(samples, decimal=9, base=2)
|
|
|
np.testing.assert_almost_equal(_GWTC1, _pesummary, 7)
|
|
|
diff.append(np.abs(_GWTC1 - _pesummary))
|
|
|
|
|
|
print(np.max(diff))
|
|
|
```
|
|
|
|
|
|
with output:
|
|
|
|
|
|
```bash
|
|
|
1.9302468302548337e-08
|
|
|
```
|
|
|
|
|
|
#### Review comments
|
|
|
|
|
|
Phil - Looks good. |