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

from packaging.version import Version
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()

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

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

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.014 -0.015 +0.017	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.98 +/- 0.04) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.245753
total	-10.245753

Values of statistical measures:

	statistical measures
AIC	25.197389
BIC	26.482971
DIC	24.505069
PDIC	2.002502

[5]:

../_images/notebooks_sampler_docs_7_8.png

../_images/notebooks_sampler_docs_7_9.png

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

 *****************************************************
 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)=  -19.092733372449974      +/-  0.14082216261646543
 Total Likelihood Evaluations:         6358
 Sampling finished. Exiting MultiNest

Maximum a posteriori probability (MAP) point:

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

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.247942
total	-10.247942

Values of statistical measures:

	statistical measures
AIC	25.201766
BIC	26.487349
DIC	24.289048
PDIC	1.896654
log(Z)	-8.291869

WARNING:root:Too few points to create valid contours

[6]:

../_images/notebooks_sampler_docs_9_8.png

../_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(nlive=400)
bayes_analysis.sample()

xyl.plot()
bayes_analysis.results.corner_plot()

100%|█████████▉| 4207/4208 [00:04<00:00, 1003.69it/s, +400 | bound: 11 | nc: 1 | ncall: 19217 | eff(%): 24.483 | loglstar: -10.232 | logz: -19.693 +/-  0.145 | dlogz:  0.001 >  0.409]

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.014 -0.015 +0.016	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.97 +/- 0.04) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.246335
total	-10.246335

Values of statistical measures:

	statistical measures
AIC	25.198552
BIC	26.484134
DIC	24.257026
PDIC	1.882282
log(Z)	-8.552653

[7]:

../_images/notebooks_sampler_docs_11_7.png

../_images/notebooks_sampler_docs_11_8.png

../_images/notebooks_sampler_docs_11_9.png

[8]:

bayes_analysis.set_sampler("dynesty_dynamic")
bayes_analysis.sampler.setup()

if Version(dynesty.__version__) >= Version("3.0.0"):
    bayes_analysis.sample(n_effective=None)
else:
    bayes_analysis.sample(
        stop_function=dynesty.utils.old_stopping_function, n_effective=None
    )

xyl.plot()
bayes_analysis.results.corner_plot()

100%|█████████▉| 17337/17339 [00:16<00:00, 1057.70it/s, batch: 9 | bound: 5 | nc: 1 | ncall: 38841 | eff(%): 44.608 | loglstar: -15.275 < -10.231 < -10.474 | logz: -19.925 +/-  0.074 | stop:  0.844]

Maximum a posteriori probability (MAP) point:

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

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.245761
total	-10.245761

Values of statistical measures:

	statistical measures
AIC	25.197404
BIC	26.482987
DIC	24.451251
PDIC	1.979743
log(Z)	-8.652203

[8]:

../_images/notebooks_sampler_docs_12_7.png

../_images/notebooks_sampler_docs_12_8.png

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

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:03<00:00, 170.28it/s]
fit restored to maximum of posterior
fit restored to maximum of posterior
Summary
-------
Number of Generations: 625
Number of Parameters: 2
Number of Walkers: 20
Number of Tuning Generations: 35
Scale Factor: 1.212627
Mean Integrated Autocorrelation Time: 2.73
Effective Sample Size: 4579.37
Number of Log Probability Evaluations: 65777
Effective Samples per Log Probability Evaluation: 0.06962

None

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.014 -0.017 +0.016	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.98 +/- 0.04) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.245759
total	-10.245759

Values of statistical measures:

	statistical measures
AIC	25.197401
BIC	26.482983
DIC	24.608373
PDIC	2.057013

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

sampler set to [blue]ultranest[/blue]

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

[ultranest] Explored until L=-1e+01
[ultranest] Likelihood function evaluations: 7954
[ultranest]   logZ = -19.74 +- 0.1116
[ultranest] Effective samples strategy satisfied (ESS = 990.3, 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.15 bs:0.11 tail:0.41 total:0.42 required:<0.50
[ultranest] done iterating.

fit restored to maximum of posterior
fit restored to maximum of posterior

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.014 -0.016 +0.017	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.98 +/- 0.04) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.24732
total	-10.24732

Values of statistical measures:

	statistical measures
AIC	25.200521
BIC	26.486104
DIC	24.497682
PDIC	2.000689
log(Z)	-8.575851

[10]:

../_images/notebooks_sampler_docs_16_11.png

../_images/notebooks_sampler_docs_16_12.png

../_images/notebooks_sampler_docs_16_13.png