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

15:09:34 WARNING   The naima package is not available. Models that depend on it will not be         functions.py:51
                  available

         WARNING   The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it  functions.py:72
                  will not be available.

15:09:35 WARNING   The ebltable package is not available. Models that depend on it will not be     absorption.py:37
                  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))

15:09:38 INFO      Using Gaussian statistic (equivalent to chi^2) with the provided errors.            XYLike.py:62

15:09:39 INFO      Using Gaussian statistic (equivalent to chi^2) with the provided errors.            XYLike.py:62

../_images/notebooks_sampler_docs_5_2.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                                                    bayesian_analysis.py:233

15:09:44 INFO      Mean acceptance fraction: 0.7127                                            emcee_sampler.py:157

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	(1.007 +/- 0.005) x 10^-1	rad / keV


Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-5.789283
total	-5.789283


Values of statistical measures:

	statistical measures
AIC	16.284448
BIC	17.570030
DIC	15.662018
PDIC	2.039614

[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()

15:09:46 INFO      sampler set to multinest                                                bayesian_analysis.py:233

  analysing data from chains/fit-.txt
Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.89 +/- 0.16) x 10^-1	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(1.007 +/- 0.005) x 10^-1	rad / keV


Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-5.790174
total	-5.790174


Values of statistical measures:

	statistical measures
AIC	16.286230
BIC	17.571813
DIC	15.583248
PDIC	2.001242
log(Z)	-6.495891

15:09:47 INFO      deleting the chain directory chains                                     multinest_sampler.py:261

[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)=  -14.957341088739314      +/-  0.14286640001680764
 Total Likelihood Evaluations:         5786
 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()

15:09:48 INFO      sampler set to dynesty_nested                                           bayesian_analysis.py:233

4154it [00:10, 388.18it/s, +400 | bound: 8 | nc: 1 | ncall: 20060 | eff(%): 23.164 | loglstar:   -inf < -5.800 <    inf | logz: -15.134 +/-  0.144 | dlogz:  0.001 >  0.409]

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.89 -0.17 +0.16) x 10^-1	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(1.007 +/- 0.005) x 10^-1	rad / keV


Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-5.789555
total	-5.789555


Values of statistical measures:

	statistical measures
AIC	16.284992
BIC	17.570574
DIC	15.541130
PDIC	1.980991
log(Z)	-6.572678

[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()

15:10:01 INFO      sampler set to dynesty_dynamic                                          bayesian_analysis.py:233

7172it [00:15, 1165.46it/s, batch: 0 | bound: 12 | nc: 1 | ncall: 26585 | eff(%): 26.705 | loglstar:   -inf < -5.805 <    inf | logz: -15.066 +/-  0.128 | dlogz:  0.005 >  0.010]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

8348it [00:18, 913.47it/s, batch: 1 | bound: 3 | nc: 2 | ncall: 28145 | eff(%): 29.341 | loglstar: -7.663 < -6.050 < -6.294 | logz: -15.061 +/-  0.132 | stop:  1.747]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

9015it [00:20, 560.16it/s, batch: 2 | bound: 2 | nc: 1 | ncall: 28899 | eff(%): 30.917 | loglstar: -8.179 < -6.643 < -7.662 | logz: -15.069 +/-  0.108 | stop:  1.105]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

9816it [00:21, 624.33it/s, batch: 3 | bound: 2 | nc: 2 | ncall: 29759 | eff(%): 32.967 | loglstar: -8.556 < -5.883 < -8.172 | logz: -15.066 +/-  0.101 | stop:  1.053]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

9832it [00:22, 436.12it/s, batch: 3 | bound: 2 | nc: 1 | ncall: 29775 | eff(%): 33.021 | loglstar: -8.556 < -5.804 < -8.172 | logz: -15.066 +/-  0.101 | stop:  0.916]

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.90 -0.16 +0.17) x 10^-1	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(1.007 +/- 0.005) x 10^-1	rad / keV


Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-5.789308
total	-5.789308


Values of statistical measures:

	statistical measures
AIC	16.284498
BIC	17.570081
DIC	15.660321
PDIC	2.040791
log(Z)	-6.551229

[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()

15:10:25 INFO      sampler set to zeus                                                     bayesian_analysis.py:233

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, 38.37it/s]
Summary
-------
Number of Generations: 625
Number of Parameters: 2
Number of Walkers: 20
Number of Tuning Generations: 27
Scale Factor: 1.390873
Mean Integrated Autocorrelation Time: 2.92
Effective Sample Size: 4275.54
Number of Log Probability Evaluations: 64522
Effective Samples per Log Probability Evaluation: 0.066265

None
Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.90 +/- 0.17) x 10^-1	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(1.007 +/- 0.005) x 10^-1	rad / keV


Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-5.789189
total	-5.789189


Values of statistical measures:

	statistical measures
AIC	16.284261
BIC	17.569843
DIC	15.749748
PDIC	2.085597

[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()

15:10:43 INFO      sampler set to ultranest                                                bayesian_analysis.py:233

[ultranest] Sampling 400 live points from prior ...

[ultranest] Explored until L=-6
[ultranest] Likelihood function evaluations: 6069
[ultranest]   logZ = -15.11 +- 0.1316
[ultranest] Effective samples strategy satisfied (ESS = 976.9, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.07 nat, need <0.50 nat)
[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.43, need <0.5)
[ultranest]   logZ error budget: single: 0.14 bs:0.13 tail:0.41 total:0.43 required:<0.50
[ultranest] done iterating.
Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.90 -0.17 +0.16) x 10^-1	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(1.007 +/- 0.005) x 10^-1	rad / keV


Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-5.793097
total	-5.793097


Values of statistical measures:

	statistical measures
AIC	16.292076
BIC	17.577658
DIC	15.542795
PDIC	1.981204
log(Z)	-6.584871

[10]:

../_images/notebooks_sampler_docs_16_9.png

../_images/notebooks_sampler_docs_16_10.png

../_images/notebooks_sampler_docs_16_11.png