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

Constraining parameters and model comparison tutorials #81

Open
wants to merge 17 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all 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
3 changes: 2 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
[![Documentation Status](https://readthedocs.org/projects/imagine-code/badge/?version=latest)](https://imagine-code.readthedocs.io/en/latest/?badge=latest)
[![Build Status](https://travis-ci.com/IMAGINE-Consortium/imagine.svg?branch=master)](https://travis-ci.com/IMAGINE-Consortium/imagine)
[![Docker](https://github.com/IMAGINE-Consortium/imagine/workflows/Docker/badge.svg)](https://github.com/IMAGINE-Consortium/imagine/actions?query=workflow%3ADocker)
[![Release Version](https://img.shields.io/github/release/IMAGINE-Consortium/imagine.svg)](https://github.com/IMAGINE-Consortium/imagine/releases)

# IMAGINE

Expand All @@ -10,7 +11,7 @@ analysing the Galactic components with observables.
It is a modular framework for doing inference on generic (non)parametric models of the Galaxy.
The analysing pipeline works with MPI support and MultiNest sampler.

The Galactic magnetic field (GMF) has a huge impact on the evolution of the Milky Way,
The Galactic magnetic field (GMF) has significant impact on the evolution of the Milky Way,
yet currently there exists no standard model for it as its structure is not fully understood.
In the past many parametric GMF models of varying complexity have been developed that
all have been fitted to an individual set of observational data complicating comparability.
Expand Down
7 changes: 0 additions & 7 deletions doc/source/evidence.rst

This file was deleted.

4 changes: 2 additions & 2 deletions doc/source/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -42,8 +42,8 @@ please refer to the
installation
design
components
posterior
evidence
tutorial_fitting
tutorial_comparison
parallel

.. toctree::
Expand Down
10 changes: 0 additions & 10 deletions doc/source/posterior.rst

This file was deleted.

3 changes: 3 additions & 0 deletions doc/source/tutorial_comparison.nblink
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
{
"path": "../../tutorials/tutorial_comparison.ipynb"
}
3 changes: 3 additions & 0 deletions doc/source/tutorial_fitting.nblink
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
{
"path": "../../tutorials/tutorial_fitting.ipynb"
}
16 changes: 10 additions & 6 deletions imagine/pipelines/dynesty_pipeline.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ class DynestyPipeline(Pipeline):

This pipeline may use
:py:class:`DynamicNestedSampler <dynesty.DynamicNestedSampler>` if the
sampling parameter 'dynamic' is set to `True` or
sampling parameter 'dynamic' is set to `True` (default) or
:py:class:`NestedSampler <dynesty.NestedSampler>`
if 'dynamic` is False (default).

Expand All @@ -37,9 +37,9 @@ class DynestyPipeline(Pipeline):
Sampling controllers
--------------------
dynamic : bool
If `True`, use
If `True` (default), use
:py:class:`dynesty.DynamicNestedSampler` otherwise uses
:py:class:`dynesty.NestedSampler`
:py:class:`dynesty.NestedSampler`.
dlogz : float
Iteration will stop, in the `dynamic==False` case,
when the estimated contribution of the
Expand Down Expand Up @@ -291,8 +291,12 @@ def call(self, **kwargs):

self.results = self.sampler.results

self._samples_array = self.results['samples']
self._evidence = self.results['logz']
self._evidence_err = self.results['logzerr']
# Creates equal weighted samples from the weighted samples
samples = self.results.samples
weights = np.exp(self.results.logwt - self.results.logz[-1])

self._samples_array = dynesty.utils.resample_equal(samples, weights)
self._evidence = self.results['logz'][-1]
self._evidence_err = self.results['logzerr'][-1]

return self.results
4 changes: 3 additions & 1 deletion imagine/pipelines/emcee_pipeline.py
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,7 @@ def call(self, **kwargs):

old_tau = np.inf
nsteps = 0
nsteps_0 = self.sampler.iteration

# Iterates trying to reach convergence
while nsteps < params['max_nsteps']:
Expand All @@ -145,9 +146,10 @@ def call(self, **kwargs):
self.tau * params['convergence_factor'] < self.sampler.iteration)
self.converged &= np.all(np.abs(old_tau - self.tau) / self.tau < 0.01)
if self.converged:
print('pipeline converged')
break
old_tau = self.tau
nsteps = self.sampler.iteration
nsteps = self.sampler.iteration -nsteps_0

burnin = int(params['burnin_factor'] * np.max(self.tau))
thin = int(params['thin_factor'] * np.min(self.tau))
Expand Down
47 changes: 47 additions & 0 deletions imagine/pipelines/pipeline.py
Original file line number Diff line number Diff line change
Expand Up @@ -142,6 +142,7 @@ def __init__(self, *, simulator, factory_list, likelihood, ensemble_size=1,
self._MAP_simulation = None
self._MAP_model = None
self._MAP = None
self._BIC = None

# Report settings
self.show_summary_reports = show_summary_reports
Expand Down Expand Up @@ -558,6 +559,52 @@ def MAP_simulation(self):
self._MAP_simulation = self.simulator(self.MAP_model)
return self._MAP_simulation

@property
def BIC(self):
if self._BIC is None:
self.get_BIC()
return self._BIC

def get_BIC(self):
r"""
Computes the Bayesian Information Criterion (BIC), this is done
using the simple expression

.. math::

BIC = k \ln n - 2 \ln L(\theta_{\rm MAP})

where :math:`k` is the number of parameters, :math:`n` is the total
number of data points, and :math:`\hat{L}` is the likelihood function
at a reference point :math:`\theta_{\rm MAP}.

Traditionally, this information criterion uses maximum likelihood as
a reference point (i.e. :math:`\theta_{\rm MLE}`).
By default, however, this method uses the likelihod
at the MAP as the reference. motivated by the heuristic that, if the
choice of prior is a sensible one, the MAP a better representation of
the model performance than the MLE.
"""
k = len(self._active_parameters)

# Counts the total number of datapoints/pixels
n = 0
for key in self.likelihood.measurement_dict:
obs = self.likelihood.measurement_dict[key]
n += obs.shape[1]

# Gets MAP
if self._MAP is None:
self.get_MAP()
# Evaluates the log likelihood at that value
# (the units need to be removed)
MAP = [param.value for param in self._MAP]
logL = self._likelihood_function(MAP)

self._BIC = k*np.log(n) - 2*logL

return self._BIC

@property
def samples(self):
"""
Expand Down
6 changes: 5 additions & 1 deletion imagine/pipelines/ultranest_pipeline.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,6 @@
# All declaration
__all__ = ['UltranestPipeline']


# %% CLASS DEFINITIONS
class UltranestPipeline(Pipeline):
"""
Expand Down Expand Up @@ -151,6 +150,10 @@ def call(self, **kwargs):
# Creates directory, if needed
os.makedirs(ultranest_dir, exist_ok=True)

# Adjusts UltraNest's logger
ultranest_logger = log.getLogger("ultranest")
ultranest_logger.setLevel(log.INFO)

# Runs UltraNest
sampler = ultranest.ReactiveNestedSampler(
param_names=list(self.active_parameters),
Expand All @@ -161,6 +164,7 @@ def call(self, **kwargs):
wrapped_params=self.wrapped_parameters,
**init_params)


self.results = sampler.run(viz_callback=ultranest.viz.nicelogger,
**run_params)

Expand Down
2 changes: 1 addition & 1 deletion imagine/simulators/test_simulator.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ class TestSimulator(Simulator):
performing the integration, i.e. computes:

.. math ::
$t(x,y,z) = B_y\,n_e\,$
t(x,y,z) = B_y\,n_e\,

"""

Expand Down
74 changes: 74 additions & 0 deletions tutorials/tutorial_comparison.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,74 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Model comparison\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Model evidence\n",
"\n",
"In the context of Bayesian statistics, the most rigorous way of comparing a set of models is looking at the Bayes *evidence* associated with each model.\n",
"Given a dataset $\\mathcal D$, model $\\mathcal M$ and a set of parameter $\\theta$, the evidence or\n",
"[marginal likelihood](https://en.wikipedia.org/wiki/Marginal_likelihood) is given by\n",
"\n",
"$$\n",
" \\mathcal{Z} = P(\\mathcal D|\\mathcal M) = \\int_\\theta P(\\mathcal D | \\theta, \\mathcal M) P(\\theta | \\mathcal M) \\mathrm{d}\\theta \n",
" = \\int_\\theta \\mathcal{L}_{\\mathcal M}(\\theta) \\pi(\\theta) \\mathrm{d}\\theta\n",
"$$\n",
"\n",
"To compare different models, one wants to compute the [Bayes factor](https://en.wikipedia.org/wiki/Bayes_factor), which is given by\n",
"\n",
"$$\n",
" R = \\frac{P(\\mathcal{M}_1 | d)}{P(\\mathcal{M}_2 | d)} \n",
" = \\frac{P(\\mathcal{D} | \\mathcal{M}_1) P(\\mathcal{M}_1)}{P(\\mathcal{D} | \\mathcal{M}_2) P(\\mathcal{M}_2)}\n",
" = \\frac{\\mathcal{Z}_1}{\\mathcal{Z}_2} \\times \\frac{P(\\mathcal{M}_1)}{P(\\mathcal{M}_2)} . \n",
"$$\n",
"\n",
"If there is no a priori reason for preferring one model over the other (as it is often the case), the prior ratio $P(\\mathcal{M}_1)/P(\\mathcal{M}_2)$\n",
"becomes unity and the Bayes factor becomes the ratio of evidences. \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Information criteria"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python (imagine)",
"language": "python",
"name": "imagine"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.8"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Loading