Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Normalised stats #396

Merged
merged 6 commits into from
Sep 27, 2024
Merged
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
anesthetic: nested sampling post-processing
===========================================
:Authors: Will Handley and Lukas Hergt
:Version: 2.8.14
:Version: 2.9.0
:Homepage: https://github.com/handley-lab/anesthetic
:Documentation: http://anesthetic.readthedocs.io/

Expand Down
2 changes: 1 addition & 1 deletion anesthetic/_version.py
Original file line number Diff line number Diff line change
@@ -1 +1 @@
__version__ = '2.8.14'
__version__ = '2.9.0'
37 changes: 32 additions & 5 deletions anesthetic/samples.py
Original file line number Diff line number Diff line change
Expand Up @@ -761,7 +761,7 @@ def ns_output(self, *args, **kwargs):
" as well as average loglikelihoods: help(samples.logL_P)"
)

def stats(self, nsamples=None, beta=None):
def stats(self, nsamples=None, beta=None, norm=None):
r"""Compute Nested Sampling statistics.

Using nested sampling we can compute:
Expand Down Expand Up @@ -822,20 +822,27 @@ def stats(self, nsamples=None, beta=None):
beta : float, array-like, optional
inverse temperature(s) beta=1/kT. Default self.beta

norm : Series, :class:`Samples`, optional
:meth:`NestedSamples.stats` output used for normalisation.
Can be either a Series of mean values or Samples produced with
matching `nsamples` and `beta`. In addition to the columns
['logZ', 'D_KL', 'logL_P', 'd_G'], this adds the normalised
versions ['Delta_logZ', 'Delta_D_KL', 'Delta_logL_P', 'Delta_d_G'].

lukashergt marked this conversation as resolved.
Show resolved Hide resolved
Returns
-------
if beta is scalar and nsamples is None:
Series, index ['logZ', 'd_G', 'DK_L', 'logL_P']
Series, index ['logZ', 'd_G', 'D_KL', 'logL_P']
elif beta is scalar and nsamples is int:
:class:`Samples`, index range(nsamples),
columns ['logZ', 'd_G', 'DK_L', 'logL_P']
columns ['logZ', 'd_G', 'D_KL', 'logL_P']
elif beta is array-like and nsamples is None:
:class:`Samples`, index beta,
columns ['logZ', 'd_G', 'DK_L', 'logL_P']
columns ['logZ', 'd_G', 'D_KL', 'logL_P']
elif beta is array-like and nsamples is int:
:class:`Samples`, index :class:`pandas.MultiIndex` the product of
beta and range(nsamples)
columns ['logZ', 'd_G', 'DK_L', 'logL_P']
columns ['logZ', 'd_G', 'D_KL', 'logL_P']
"""
logw = self.logw(nsamples, beta)
if nsamples is None and beta is None:
Expand All @@ -861,6 +868,26 @@ def stats(self, nsamples=None, beta=None):
samples.set_label('d_G', r'$d_\mathrm{G}$')

samples.label = self.label

if norm is not None:
lukashergt marked this conversation as resolved.
Show resolved Hide resolved
samples['Delta_logZ'] = samples['logZ'] - norm['logZ']
samples.set_label('Delta_logZ',
r"$\Delta\ln\mathcal{Z}$")

samples['Delta_D_KL'] = samples['D_KL'] - norm['D_KL']
samples.set_label('Delta_D_KL',
r"$\Delta\mathcal{D}_\mathrm{KL}$")

samples['Delta_logL_P'] = samples['logL_P'] - norm['logL_P']
samples.set_label(
'Delta_logL_P',
r"$\Delta\langle\ln\mathcal{L}\rangle_\mathcal{P}$"
)

samples['Delta_d_G'] = samples['d_G'] - norm['d_G']
samples.set_label('Delta_d_G',
r"$\Delta d_\mathrm{G}$")

return samples

def logX(self, nsamples=None):
Expand Down
54 changes: 54 additions & 0 deletions tests/test_samples.py
Original file line number Diff line number Diff line change
Expand Up @@ -946,39 +946,93 @@ def test_stats():
beta = [0., 0.5, 1.]

vals = ['logZ', 'D_KL', 'logL_P', 'd_G']
delta_vals = ['Delta_logZ', 'Delta_D_KL', 'Delta_logL_P', 'Delta_d_G']

labels = [r'$\ln\mathcal{Z}$',
r'$\mathcal{D}_\mathrm{KL}$',
r'$\langle\ln\mathcal{L}\rangle_\mathcal{P}$',
r'$d_\mathrm{G}$']
delta_labels = [r'$\Delta\ln\mathcal{Z}$',
r'$\Delta\mathcal{D}_\mathrm{KL}$',
r'$\Delta\langle\ln\mathcal{L}\rangle_\mathcal{P}$',
r'$\Delta d_\mathrm{G}$']

stats = pc.stats()
assert isinstance(stats, WeightedLabelledSeries)
assert_array_equal(stats.drop_labels().index, vals)
assert_array_equal(stats.get_labels(), labels)

stats = pc.stats(norm=pc.stats())
assert isinstance(stats, WeightedLabelledSeries)
assert_array_equal(stats.drop_labels().index, vals + delta_vals)
assert_array_equal(stats.get_labels(), labels + delta_labels)

stats = pc.stats(nsamples=nsamples)
assert isinstance(stats, WeightedLabelledDataFrame)
assert_array_equal(stats.drop_labels().columns, vals)
assert_array_equal(stats.get_labels(), labels)
assert stats.index.name == 'samples'
assert_array_equal(stats.index, range(nsamples))

stats = pc.stats(nsamples=nsamples, norm=pc.stats())
assert isinstance(stats, WeightedLabelledDataFrame)
assert_array_equal(stats.drop_labels().columns, vals + delta_vals)
assert_array_equal(stats.get_labels(), labels + delta_labels)
assert stats.index.name == 'samples'
assert_array_equal(stats.index, range(nsamples))

stats = pc.stats(nsamples=nsamples, norm=pc.stats(nsamples=nsamples))
assert isinstance(stats, WeightedLabelledDataFrame)
assert_array_equal(stats.drop_labels().columns, vals + delta_vals)
assert_array_equal(stats.get_labels(), labels + delta_labels)
assert stats.index.name == 'samples'
assert_array_equal(stats.index, range(nsamples))

stats = pc.stats(beta=beta)
assert isinstance(stats, WeightedLabelledDataFrame)
assert_array_equal(stats.drop_labels().columns, vals)
assert_array_equal(stats.get_labels(), labels)
assert stats.index.name == 'beta'
assert_array_equal(stats.index, beta)

stats = pc.stats(beta=beta, norm=pc.stats())
assert isinstance(stats, WeightedLabelledDataFrame)
assert_array_equal(stats.drop_labels().columns, vals + delta_vals)
assert_array_equal(stats.get_labels(), labels + delta_labels)
assert stats.index.name == 'beta'
assert_array_equal(stats.index, beta)

stats = pc.stats(beta=beta, norm=pc.stats(beta=beta))
assert isinstance(stats, WeightedLabelledDataFrame)
assert_array_equal(stats.drop_labels().columns, vals + delta_vals)
assert_array_equal(stats.get_labels(), labels + delta_labels)
assert stats.index.name == 'beta'
assert_array_equal(stats.index, beta)

stats = pc.stats(nsamples=nsamples, beta=beta)
assert isinstance(stats, WeightedLabelledDataFrame)
assert_array_equal(stats.drop_labels().columns, vals)
assert_array_equal(stats.get_labels(), labels)
assert stats.index.names == ['beta', 'samples']
assert stats.index.levshape == (len(beta), nsamples)

stats = pc.stats(nsamples=nsamples, beta=beta, norm=pc.stats())
assert isinstance(stats, WeightedLabelledDataFrame)
assert_array_equal(stats.drop_labels().columns, vals + delta_vals)
assert_array_equal(stats.get_labels(), labels + delta_labels)
assert stats.index.names == ['beta', 'samples']
assert stats.index.levshape == (len(beta), nsamples)

stats = pc.stats(nsamples=nsamples, beta=beta,
norm=pc.stats(nsamples=nsamples, beta=beta))
assert isinstance(stats, WeightedLabelledDataFrame)
assert_array_equal(stats.drop_labels().columns, vals + delta_vals)
assert_array_equal(stats.get_labels(), labels + delta_labels)
assert stats.index.names == ['beta', 'samples']
assert stats.index.levshape == (len(beta), nsamples)

for beta in [1., 0., 0.5]:
np.random.seed(42)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do you need a seed here?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I managed to get failing tests in what follows, because of changes to the seed from the tests I added above... Nothing dramatic, just bad chance in combination with the tests checking against 3*std or similar. I could have fiddled around with tolerances, but just trying a new seed there was enough, so I did not dig deeper...

Copy link
Collaborator

@AdamOrmondroyd AdamOrmondroyd Sep 27, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok, I'll reapprove when you've done the version number ofc this is a minor bump

pc.beta = beta
n = 1000
PC = pc.stats(n, beta)
Expand Down
Loading