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.7188
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.87 +/- 0.17) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.98 +/- 0.05) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -6.019459
total -6.019459

Values of statistical measures:

statistical measures
AIC 16.744801
BIC 18.030383
DIC 15.897667
PDIC 1.921855

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 (9.89 -0.17 +0.18) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.99 +/- 0.05) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -6.020122
total -6.020122

Values of statistical measures:

statistical measures
AIC 16.746126
BIC 18.031708
DIC 16.045475
PDIC 2.002528
log(Z) -6.663083
[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()

[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)=  -15.342315101991641      +/-  0.14426542337390316
 Total Likelihood Evaluations:         5633
 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
4114it [00:05, 735.12it/s, +400 | bound: 10 | nc: 1 | ncall: 20180 | eff(%): 22.369 | loglstar:   -inf < -6.025 <    inf | logz: -15.260 +/-  0.202 | dlogz:  0.001 >  0.409]
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.88 -0.16 +0.17) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.98 +/- 0.05) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -6.020453
total -6.020453

Values of statistical measures:

statistical measures
AIC 16.746788
BIC 18.032371
DIC 15.870814
PDIC 1.915349
log(Z) -6.627171

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
24306it [03:52, 104.49it/s, batch: 12 | bound: 60 | nc: 1 | ncall: 58073 | eff(%): 41.854 | loglstar: -12.979 < -6.025 < -6.381 | logz: -15.357 +/-  0.182 | stop:  0.966]
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.88 +/- 0.17) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.98 +/- 0.05) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -6.019564
total -6.019564

Values of statistical measures:

statistical measures
AIC 16.745010
BIC 18.030592
DIC 16.049929
PDIC 2.004820
log(Z) -6.670842

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:16<00:00, 37.51it/s]
Summary
-------
Number of Generations: 625
Number of Parameters: 2
Number of Walkers: 20
Number of Tuning Generations: 21
Scale Factor: 1.113809
Mean Integrated Autocorrelation Time: 3.0
Effective Sample Size: 4166.6
Number of Log Probability Evaluations: 66072
Effective Samples per Log Probability Evaluation: 0.063062
None
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.88 +/- 0.17) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.99 +/- 0.05) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -6.01981
total -6.01981

Values of statistical measures:

statistical measures
AIC 16.745502
BIC 18.031084
DIC 16.071081
PDIC 2.015786

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=-6
[ultranest] Likelihood function evaluations: 8866
[ultranest]   logZ = -15.41 +- 0.119
[ultranest] Effective samples strategy satisfied (ESS = 975.7, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.45+-0.07 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.42, need <0.5)
[ultranest]   logZ error budget: single: 0.14 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 (9.89 +/- 0.17) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.99 +/- 0.05) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -6.021
total -6.021

Values of statistical measures:

statistical measures
AIC 16.747882
BIC 18.033464
DIC 16.053768
PDIC 2.007373
log(Z) -6.702627

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