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.6483
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.53 -0.09 +0.5) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.5 +0.4 +0.5) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -3.147496
total -3.147496

Values of statistical measures:

statistical measures
AIC 11.000875
BIC 12.286457
DIC -9.178552
PDIC -61.230253
[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.76 -0.25 +0.26) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.95 -0.05 +0.04) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -3.147351
total -3.147351

Values of statistical measures:

statistical measures
AIC 11.000585
BIC 12.286167
DIC 10.486654
PDIC 2.095988
log(Z) -5.259971
[INFO    ] deleting the chain directory chains
WARNING:root:Too few points to create valid contours
[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)=  -12.111531093828289      +/-  0.14041035407109745
 Total Likelihood Evaluations:         5169
 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()
[INFO    ] sampler set to dynesty_nested
4014it [00:10, 372.82it/s, +400 | bound: 8 | nc: 1 | ncall: 19641 | eff(%): 22.941 | loglstar:   -inf < -3.171 <    inf | logz: -12.168 +/-  0.142 | dlogz:  0.001 >  0.409]
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.76 -0.26 +0.27) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.95 +/- 0.05) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -3.147457
total -3.147457

Values of statistical measures:

statistical measures
AIC 11.000796
BIC 12.286378
DIC 10.513147
PDIC 2.109102
log(Z) -5.284329
[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
7210it [00:16, 1536.62it/s, batch: 0 | bound: 12 | nc: 1 | ncall: 26463 | eff(%): 27.225 | loglstar:   -inf < -3.172 <    inf | logz: -12.008 +/-  0.125 | dlogz:  0.000 >  0.010]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

8262it [00:19, 935.83it/s, batch: 1 | bound: 3 | nc: 1 | ncall: 27969 | eff(%): 29.341 | loglstar: -5.150 < -3.339 < -3.632 | logz: -12.007 +/-  0.129 | stop:  1.472]
WARNING DeprecationWarning: This an old stopping function that will be removed in future releases

8441it [00:19, 426.72it/s, batch: 1 | bound: 3 | nc: 1 | ncall: 28159 | eff(%): 29.976 | loglstar: -5.150 < -3.172 < -3.632 | logz: -12.007 +/-  0.129 | stop:  0.987]
Maximum a posteriori probability (MAP) point:

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

Values of -log(posterior) at the minimum:

-log(posterior)
demo -3.147337
total -3.147337

Values of statistical measures:

statistical measures
AIC 11.000557
BIC 12.286139
DIC 10.230251
PDIC 1.967197
log(Z) -5.206736
[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.59it/s]
Summary
-------
Number of Generations: 625
Number of Parameters: 2
Number of Walkers: 20
Number of Tuning Generations: 37
Scale Factor: 1.403902
Mean Integrated Autocorrelation Time: 2.81
Effective Sample Size: 4448.72
Number of Log Probability Evaluations: 64789
Effective Samples per Log Probability Evaluation: 0.068665
None
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.75 +/- 0.25) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.95 -0.05 +0.04) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -3.14734
total -3.14734

Values of statistical measures:

statistical measures
AIC 11.000563
BIC 12.286145
DIC 10.238191
PDIC 1.971085
[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=-3
[ultranest] Likelihood function evaluations: 9806
[ultranest]   logZ = -11.92 +- 0.1091
[ultranest] Effective samples strategy satisfied (ESS = 970.7, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.08 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.11 tail:0.40 total:0.42 required:<0.50
[ultranest] done iterating.
Maximum a posteriori probability (MAP) point:

result unit
parameter
demo.spectrum.main.Sin.K (9.76 -0.26 +0.25) x 10^-1 1 / (cm2 keV s)
demo.spectrum.main.Sin.f (9.95 -0.04 +0.05) x 10^-2 rad / keV

Values of -log(posterior) at the minimum:

-log(posterior)
demo -3.147647
total -3.147647

Values of statistical measures:

statistical measures
AIC 11.001176
BIC 12.286758
DIC 10.168432
PDIC 1.936168
log(Z) -5.163283
[10]:
../_images/notebooks_sampler_docs_16_8.png
../_images/notebooks_sampler_docs_16_9.png
../_images/notebooks_sampler_docs_16_10.png