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 *
from threeML.plugins.XYLike import XYLike

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

21:34:38 WARNING   The naima package is not available. Models that depend on it will not be         functions.py:47
                  available

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

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

21:34:40 INFO      Using Gaussian statistic (equivalent to chi^2) with the provided errors.            XYLike.py:93

21:34:41 INFO      Using Gaussian statistic (equivalent to chi^2) with the provided errors.            XYLike.py:93

../_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:202

21:34:44 INFO      Mean acceptance fraction: 0.7112                                            emcee_sampler.py:157

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

21:34:45 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.002 -0.023 +0.024	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(9.968 -0.025 +0.026) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-7.711225
total	-7.711225

Values of statistical measures:

	statistical measures
AIC	20.128332
BIC	21.413915
DIC	19.355913
PDIC	1.966289

[5]:

../_images/notebooks_sampler_docs_7_12.png

../_images/notebooks_sampler_docs_7_13.png

../_images/notebooks_sampler_docs_7_14.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                                                bayesian_analysis.py:202

 *****************************************************
 MultiNest v3.10
 Copyright Farhan Feroz & Mike Hobson
 Release Jul 2015

 no. of live points =  400
 dimensionality =    2
 *****************************************************
  analysing data from chains/fit-.txt ln(ev)=  -17.273050735622686      +/-  0.14614904848402743
 Total Likelihood Evaluations:         6354
 Sampling finished. Exiting MultiNest

21:34:47 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.001 -0.021 +0.023	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(9.967 -0.028 +0.030) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-7.712639
total	-7.712639

Values of statistical measures:

	statistical measures
AIC	20.131160
BIC	21.416743
DIC	19.482007
PDIC	2.029519
log(Z)	-7.501591

         INFO      deleting the chain directory chains                                     multinest_sampler.py:255

WARNING:root:Too few points to create valid contours

[6]:

../_images/notebooks_sampler_docs_9_12.png

../_images/notebooks_sampler_docs_9_13.png

../_images/notebooks_sampler_docs_9_14.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                                           bayesian_analysis.py:202

4269it [00:04, 888.07it/s, +400 | bound: 9 | nc: 1 | ncall: 19188 | eff(%): 24.851 | loglstar:   -inf < -7.710 <    inf | logz: -17.336 +/-  0.147 | dlogz:  0.001 >  0.409]

21:34:52 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.002 -0.022 +0.023	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(9.967 -0.026 +0.027) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-7.711237
total	-7.711237

Values of statistical measures:

	statistical measures
AIC	20.128356
BIC	21.413938
DIC	19.387834
PDIC	1.982629
log(Z)	-7.528871

[7]:

../_images/notebooks_sampler_docs_11_10.png

../_images/notebooks_sampler_docs_11_11.png

../_images/notebooks_sampler_docs_11_12.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()

21:34:53 INFO      sampler set to dynesty_dynamic                                          bayesian_analysis.py:202

7684it [00:08, 1659.86it/s, batch: 0 | bound: 13 | nc: 1 | ncall: 27499 | eff(%): 27.923 | loglstar:   -inf < -7.710 <    inf | logz: -17.492 +/-  0.133 | dlogz:  0.000 >  0.010]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

8665it [00:09, 1237.50it/s, batch: 1 | bound: 3 | nc: 1 | ncall: 28870 | eff(%): 29.863 | loglstar: -9.384 < -7.847 < -8.199 | logz: -17.492 +/-  0.136 | stop:  1.417]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

9426it [00:10, 1000.75it/s, batch: 2 | bound: 2 | nc: 1 | ncall: 29719 | eff(%): 31.696 | loglstar: -9.912 < -7.780 < -9.381 | logz: -17.498 +/-  0.113 | stop:  1.113]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

9702it [00:11, 599.64it/s, batch: 3 | bound: 2 | nc: 1 | ncall: 30020 | eff(%): 31.951 | loglstar: -10.367 < -9.192 < -9.911 | logz: -17.500 +/-  0.105 | stop:  1.073]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

10038it [00:12, 809.09it/s, batch: 3 | bound: 2 | nc: 1 | ncall: 30365 | eff(%): 33.058 | loglstar: -10.367 < -7.714 < -9.911 | logz: -17.500 +/-  0.105 | stop:  0.869]

21:35:06 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.002 -0.022 +0.023	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(9.968 -0.027 +0.026) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-7.711254
total	-7.711254

Values of statistical measures:

	statistical measures
AIC	20.128389
BIC	21.413972
DIC	19.400741
PDIC	1.989119
log(Z)	-7.600641

[8]:

../_images/notebooks_sampler_docs_12_10.png

../_images/notebooks_sampler_docs_12_11.png

../_images/notebooks_sampler_docs_12_12.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()

21:35:07 INFO      sampler set to zeus                                                     bayesian_analysis.py:202

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:11<00:00, 55.80it/s]

21:35:19 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Summary
-------
Number of Generations: 625
Number of Parameters: 2
Number of Walkers: 20
Number of Tuning Generations: 44
Scale Factor: 1.34706
Mean Integrated Autocorrelation Time: 3.33
Effective Sample Size: 3751.67
Number of Log Probability Evaluations: 64814
Effective Samples per Log Probability Evaluation: 0.057884

None

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.002 -0.023 +0.021	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(9.968 -0.028 +0.026) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-7.711303
total	-7.711303

Values of statistical measures:

	statistical measures
AIC	20.128488
BIC	21.414070
DIC	19.285606
PDIC	1.931198

[9]:

../_images/notebooks_sampler_docs_14_12.png

../_images/notebooks_sampler_docs_14_13.png

../_images/notebooks_sampler_docs_14_14.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()

21:35:20 INFO      sampler set to ultranest                                                bayesian_analysis.py:202

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

[ultranest] Explored until L=-8
[ultranest] Likelihood function evaluations: 10194
[ultranest]   logZ = -17.01 +- 0.09287
[ultranest] Effective samples strategy satisfied (ESS = 965.7, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.06 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.09 tail:0.41 total:0.42 required:<0.50
[ultranest] done iterating.

21:35:29 INFO      fit restored to maximum of posterior                                         sampler_base.py:178

         INFO      fit restored to maximum of posterior                                         sampler_base.py:178

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.002 -0.023 +0.024	1 / (cm2 keV s)
demo.spectrum.main.Sin.f	(9.966 -0.026 +0.027) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-7.712662
total	-7.712662

Values of statistical measures:

	statistical measures
AIC	20.131206
BIC	21.416788
DIC	19.402409
PDIC	1.988829
log(Z)	-7.396792

[10]:

../_images/notebooks_sampler_docs_16_12.png

../_images/notebooks_sampler_docs_16_13.png

../_images/notebooks_sampler_docs_16_14.png