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.053 -0.025 +0.027	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.90 +/- 0.05) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.259644
total	-4.259644

Values of statistical measures:

	statistical measures
AIC	13.225171
BIC	14.510753
DIC	12.614149
PDIC	2.047563

[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)=  -13.716250509596829      +/-  0.14497824299534093
 Total Likelihood Evaluations:         4866
 Sampling finished. Exiting MultiNest

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.055 -0.029 +0.027	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.90 -0.05 +0.06) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.263485
total	-4.263485

Values of statistical measures:

	statistical measures
AIC	13.232853
BIC	14.518435
DIC	12.720367
PDIC	2.099633
log(Z)	-5.956892

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%|█████████▉| 4108/4110 [00:04<00:00, 1019.29it/s, +400 | bound: 12 | nc: 1 | ncall: 18822 | eff(%): 24.471 | loglstar: -4.207 | logz: -13.427 +/-  0.143 | dlogz:  0.001 >  0.409]

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.054 -0.025 +0.023	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.89 -0.05 +0.06) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.259473
total	-4.259473

Values of statistical measures:

	statistical measures
AIC	13.224829
BIC	14.510411
DIC	12.582453
PDIC	2.030366
log(Z)	-5.831368

[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%|█████████▉| 16497/16499 [00:14<00:00, 1178.08it/s, batch: 8 | bound: 5 | nc: 1 | ncall: 38705 | eff(%): 42.597 | loglstar: -9.127 < -4.207 < -4.493 | logz: -13.398 +/-  0.073 | stop:  0.884]

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.054 -0.026 +0.025	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.89 -0.05 +0.06) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.259557
total	-4.259557

Values of statistical measures:

	statistical measures
AIC	13.224996
BIC	14.510578
DIC	12.519402
PDIC	1.997839
log(Z)	-5.816595

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

Initialising ensemble of 20 walkers...
Sampling progress : 100%|██████████| 625/625 [00:03<00:00, 182.01it/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: 23
Scale Factor: 1.016353
Mean Integrated Autocorrelation Time: 2.87
Effective Sample Size: 4348.68
Number of Log Probability Evaluations: 66530
Effective Samples per Log Probability Evaluation: 0.065364

None

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	1.053 -0.026 +0.025	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.90 +/- 0.05) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.259346
total	-4.259346

Values of statistical measures:

	statistical measures
AIC	13.224574
BIC	14.510157
DIC	12.522953
PDIC	2.000754

[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=-4
[ultranest] Likelihood function evaluations: 9301
[ultranest]   logZ = -13.31 +- 0.1168
[ultranest] Effective samples strategy satisfied (ESS = 984.9, 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.12 tail:0.40 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.055 -0.028 +0.023	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.89 -0.05 +0.06) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.26052
total	-4.26052

Values of statistical measures:

	statistical measures
AIC	13.226922
BIC	14.512504
DIC	12.758363
PDIC	2.114010
log(Z)	-5.785510

[10]:

../_images/notebooks_sampler_docs_16_11.png

../_images/notebooks_sampler_docs_16_12.png

../_images/notebooks_sampler_docs_16_13.png