Bayesian Sampler Examples

Examples of running each sampler avaiable in 3ML.

Before, that, let’s discuss setting up configuration default sampler with default parameters. We can set in our configuration a default algorithm and default setup parameters for the samplers. This can ease fitting when we are doing exploratory data analysis.

With any of the samplers, you can pass keywords to access their setups. Read each pacakges documentation for more details.

[1]:
from threeML import *

import numpy as np

from jupyterthemes import jtplot

%matplotlib inline
jtplot.style(context="talk", fscale=1, ticks=True, grid=False)
silence_warnings()
set_threeML_style()
Welcome to JupyROOT 6.22/02
[WARNING ] The naima package is not available. Models that depend on it will not be available
[WARNING ] The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it will not be available.
[WARNING ] The ebltable package is not available. Models that depend on it will not be available

WARNING RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject


WARNING RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject


WARNING RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject

[2]:
threeML_config.bayesian.default_sampler
[2]:
<Sampler.emcee: 'emcee'>
[3]:
threeML_config.bayesian.emcee_setup
[3]:
{'n_burnin': None, 'n_iterations': 500, 'n_walkers': 50, 'seed': 5123}

If you simply run bayes_analysis.sample() the default sampler and its default parameters will be used.

Let’s make some data to fit.

[4]:
sin = Sin(K=1, f=0.1)
sin.phi.fix = True
sin.K.prior = Log_uniform_prior(lower_bound=0.5, upper_bound=1.5)
sin.f.prior = Uniform_prior(lower_bound=0, upper_bound=0.5)

model = Model(PointSource("demo", 0, 0, spectral_shape=sin))

x = np.linspace(-2 * np.pi, 4 * np.pi, 20)
yerr = np.random.uniform(0.01, 0.2, 20)


xyl = XYLike.from_function("demo", sin, x, yerr)
xyl.plot()

bayes_analysis = BayesianAnalysis(model, DataList(xyl))
[INFO    ] Using Gaussian statistic (equivalent to chi^2) with the provided errors.
[INFO    ] Using Gaussian statistic (equivalent to chi^2) with the provided errors.
../_images/notebooks_sampler_docs_5_1.png

emcee

[5]:
bayes_analysis.set_sampler("emcee")
bayes_analysis.sampler.setup(n_walkers=20, n_iterations=500)
bayes_analysis.sample()

xyl.plot()
bayes_analysis.results.corner_plot()
[INFO    ] sampler set to emcee
[INFO    ] Mean acceptance fraction: 0.7101999999999998
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K 1.000 +/- 0.008 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (1.0029 -0.0018 +0.0019) x 10^-1 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -10.912001
total -10.912001

Values of statistical measures:

statistical measures
AIC 26.529884
BIC 27.815466
DIC 25.986943
PDIC 2.078197

WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()

[5]:
../_images/notebooks_sampler_docs_7_10.png
../_images/notebooks_sampler_docs_7_11.png
../_images/notebooks_sampler_docs_7_12.png

multinest

[6]:
bayes_analysis.set_sampler("multinest")
bayes_analysis.sampler.setup(n_live_points=400, resume=False, auto_clean=True)
bayes_analysis.sample()

xyl.plot()
bayes_analysis.results.corner_plot()
[INFO    ] sampler set to multinest
  analysing data from chains/fit-.txt
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K 1.001 +/- 0.008 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (1.0029 -0.0020 +0.0019) x 10^-1 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -10.911873
total -10.911873

Values of statistical measures:

statistical measures
AIC 26.529628
BIC 27.815211
DIC 25.880588
PDIC 2.020462
log(Z) -9.499729
[INFO    ] deleting the chain directory chains

WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()

WARNING:root:Too few points to create valid contours
[6]:
../_images/notebooks_sampler_docs_9_8.png
 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  400
 dimensionality =    2
 *****************************************************
 ln(ev)=  -21.873934532674973      +/-  0.15761873368553861
 Total Likelihood Evaluations:         5927
 Sampling finished. Exiting MultiNest
../_images/notebooks_sampler_docs_9_10.png
../_images/notebooks_sampler_docs_9_11.png

dynesty

[7]:
bayes_analysis.set_sampler("dynesty_nested")
bayes_analysis.sampler.setup(n_live_points=400)
bayes_analysis.sample()

xyl.plot()
bayes_analysis.results.corner_plot()
[INFO    ] sampler set to dynesty_nested
4876it [00:06, 808.38it/s, +400 | bound: 12 | nc: 1 | ncall: 21180 | eff(%): 24.910 | loglstar:   -inf < -10.912 <    inf | logz: -22.048 +/-  0.225 | dlogz:  0.001 >  0.409]
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (10.00 +/- 0.08) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (1.0028 -0.0018 +0.0019) x 10^-1 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -10.912204
total -10.912204

Values of statistical measures:

statistical measures
AIC 26.530291
BIC 27.815873
DIC 25.688763
PDIC 1.932182
log(Z) -9.575204

WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()

[7]:
../_images/notebooks_sampler_docs_11_9.png
../_images/notebooks_sampler_docs_11_10.png
../_images/notebooks_sampler_docs_11_11.png
[8]:
bayes_analysis.set_sampler("dynesty_dynamic")
bayes_analysis.sampler.setup(n_live_points=400)
bayes_analysis.sample()

xyl.plot()
bayes_analysis.results.corner_plot()
[INFO    ] sampler set to dynesty_dynamic
19151it [02:37, 121.24it/s, batch: 9 | bound: 48 | nc: 1 | ncall: 48909 | eff(%): 39.156 | loglstar: -16.126 < -10.912 < -11.215 | logz: -22.127 +/-  0.202 | stop:  0.905]
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (10.00 +/- 0.08) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (1.0028 +/- 0.0019) x 10^-1 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -10.911905
total -10.911905

Values of statistical measures:

statistical measures
AIC 26.529693
BIC 27.815275
DIC 25.891711
PDIC 2.033994
log(Z) -9.613190

WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()

[8]:
../_images/notebooks_sampler_docs_12_9.png
../_images/notebooks_sampler_docs_12_10.png
../_images/notebooks_sampler_docs_12_11.png

zeus

[9]:
bayes_analysis.set_sampler("zeus")
bayes_analysis.sampler.setup(n_walkers=20, n_iterations=500)
bayes_analysis.sample()

xyl.plot()
bayes_analysis.results.corner_plot()
[INFO    ] sampler set to zeus
WARNING:root:The sampler class has been deprecated. Please use the new EnsembleSampler class.
The run method has been deprecated and it will be removed. Please use the new run_mcmc method.
Initialising ensemble of 20 walkers...
Sampling progress : 100%|██████████| 625/625 [00:15<00:00, 40.26it/s]
Summary
-------
Number of Generations: 625
Number of Parameters: 2
Number of Walkers: 20
Number of Tuning Generations: 26
Scale Factor: 1.275024
Mean Integrated Autocorrelation Time: 3.0
Effective Sample Size: 4165.48
Number of Log Probability Evaluations: 64857
Effective Samples per Log Probability Evaluation: 0.064226
None
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (10.00 +/- 0.07) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (1.0028 -0.0019 +0.0018) x 10^-1 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -10.912057
total -10.912057

Values of statistical measures:

statistical measures
AIC 26.529997
BIC 27.815579
DIC 25.705008
PDIC 1.940607

WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()

[9]:
../_images/notebooks_sampler_docs_14_9.png
../_images/notebooks_sampler_docs_14_10.png
../_images/notebooks_sampler_docs_14_11.png

ultranest

[10]:
bayes_analysis.set_sampler("ultranest")
bayes_analysis.sampler.setup(
    min_num_live_points=400, frac_remain=0.5, use_mlfriends=False
)
bayes_analysis.sample()

xyl.plot()
bayes_analysis.results.corner_plot()
[INFO    ] sampler set to ultranest
[ultranest] Sampling 400 live points from prior ...
[ultranest] Explored until L=-1e+01
[ultranest] Likelihood function evaluations: 16956
[ultranest]   logZ = -21.92 +- 0.1188
[ultranest] Effective samples strategy satisfied (ESS = 985.6, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.05 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.42, need <0.5)
[ultranest]   logZ error budget: single: 0.16 bs:0.12 tail:0.41 total:0.42 required:<0.50
[ultranest] done iterating.
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K 1.000 +/- 0.008 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (1.0028 -0.0019 +0.0020) x 10^-1 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -10.920104
total -10.920104

Values of statistical measures:

statistical measures
AIC 26.546090
BIC 27.831673
DIC 26.023635
PDIC 2.099922
log(Z) -9.531704

WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()

[10]:
../_images/notebooks_sampler_docs_16_9.png
../_images/notebooks_sampler_docs_16_10.png
../_images/notebooks_sampler_docs_16_11.png