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
import dynesty
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
[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.7153999999999999
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.56 -0.20 +0.21) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.96 +/- 0.06) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -7.582203
total -7.582203

Values of statistical measures:

statistical measures
AIC 19.870289
BIC 21.155871
DIC 19.362047
PDIC 2.098194
[5]:
../_images/notebooks_sampler_docs_7_9.png
../_images/notebooks_sampler_docs_7_10.png
../_images/notebooks_sampler_docs_7_11.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.57 +/- 0.19) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.96 -0.07 +0.06) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -7.582715
total -7.582715

Values of statistical measures:

statistical measures
AIC 19.871313
BIC 21.156895
DIC 19.394755
PDIC 2.111313
log(Z) -7.133150
[INFO    ] deleting the chain directory chains
[6]:
../_images/notebooks_sampler_docs_9_7.png
 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  400
 dimensionality =    2
 *****************************************************
 ln(ev)=  -16.424684070481497      +/-  0.13936420317415385
 Total Likelihood Evaluations:         5858
 Sampling finished. Exiting MultiNest
../_images/notebooks_sampler_docs_9_9.png
../_images/notebooks_sampler_docs_9_10.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
4014it [00:10, 375.65it/s, +400 | bound: 8 | nc: 1 | ncall: 19280 | eff(%): 23.379 | loglstar:   -inf < -7.632 <    inf | logz: -16.618 +/-  0.141 | dlogz:  0.001 >  0.409]
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.56 -0.19 +0.20) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.96 +/- 0.06) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -7.584833
total -7.584833

Values of statistical measures:

statistical measures
AIC 19.875549
BIC 21.161131
DIC 19.249353
PDIC 2.042485
log(Z) -7.216942
[7]:
../_images/notebooks_sampler_docs_11_8.png
../_images/notebooks_sampler_docs_11_9.png
../_images/notebooks_sampler_docs_11_10.png
[8]:
bayes_analysis.set_sampler("dynesty_dynamic")
bayes_analysis.sampler.setup(
    stop_function=dynesty.utils.old_stopping_function, n_effective=None
)
bayes_analysis.sample()

xyl.plot()
bayes_analysis.results.corner_plot()
[INFO    ] sampler set to dynesty_dynamic
7033it [00:16, 972.65it/s, batch: 0 | bound: 12 | nc: 1 | ncall: 25827 | eff(%): 27.030 | loglstar:   -inf < -7.631 <    inf | logz: -16.457 +/-  0.125 | dlogz:  0.004 >  0.010]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

8198it [00:18, 1048.09it/s, batch: 1 | bound: 3 | nc: 1 | ncall: 27470 | eff(%): 29.596 | loglstar: -9.635 < -7.841 < -8.108 | logz: -16.453 +/-  0.128 | stop:  1.326]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

8935it [00:20, 660.06it/s, batch: 2 | bound: 2 | nc: 1 | ncall: 28253 | eff(%): 31.520 | loglstar: -10.193 < -8.008 < -9.633 | logz: -16.439 +/-  0.105 | stop:  1.143]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

9531it [00:22, 524.82it/s, batch: 3 | bound: 2 | nc: 1 | ncall: 28885 | eff(%): 32.894 | loglstar: -10.593 < -8.042 < -10.189 | logz: -16.413 +/-  0.098 | stop:  1.011]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

9619it [00:22, 420.66it/s, batch: 3 | bound: 2 | nc: 1 | ncall: 28975 | eff(%): 33.198 | loglstar: -10.593 < -7.628 < -10.189 | logz: -16.413 +/-  0.098 | stop:  0.820]
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.55 -0.20 +0.19) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.96 +/- 0.06) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -7.582349
total -7.582349

Values of statistical measures:

statistical measures
AIC 19.870580
BIC 21.156162
DIC 19.221893
PDIC 2.028862
log(Z) -7.127890
[8]:
../_images/notebooks_sampler_docs_12_8.png
../_images/notebooks_sampler_docs_12_9.png
../_images/notebooks_sampler_docs_12_10.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.30it/s]
Summary
-------
Number of Generations: 625
Number of Parameters: 2
Number of Walkers: 20
Number of Tuning Generations: 25
Scale Factor: 1.105713
Mean Integrated Autocorrelation Time: 3.11
Effective Sample Size: 4020.95
Number of Log Probability Evaluations: 65977
Effective Samples per Log Probability Evaluation: 0.060945
None
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.56 +/- 0.19) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.96 +/- 0.06) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -7.582138
total -7.582138

Values of statistical measures:

statistical measures
AIC 19.870158
BIC 21.155740
DIC 19.133807
PDIC 1.984780
[9]:
../_images/notebooks_sampler_docs_14_8.png
../_images/notebooks_sampler_docs_14_9.png
../_images/notebooks_sampler_docs_14_10.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=-8
[ultranest] Likelihood function evaluations: 6854
[ultranest]   logZ = -16.53 +- 0.1018
[ultranest] Effective samples strategy satisfied (ESS = 993.9, 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.14 bs:0.10 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.56 +/- 0.20) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.96 +/- 0.06) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -7.583919
total -7.583919

Values of statistical measures:

statistical measures
AIC 19.873721
BIC 21.159303
DIC 19.434120
PDIC 2.133980
log(Z) -7.174210
[10]:
../_images/notebooks_sampler_docs_16_8.png
../_images/notebooks_sampler_docs_16_9.png
../_images/notebooks_sampler_docs_16_10.png