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

22:49:03 WARNING   The naima package is not available. Models that depend on it will not be         functions.py:43
                  available

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

22:49:04 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))

22:49:05 INFO      Using Gaussian statistic (equivalent to chi^2) with the provided errors.            XYLike.py:90

22:49:06 INFO      Using Gaussian statistic (equivalent to chi^2) with the provided errors.            XYLike.py:90

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

22:49:07 INFO      Mean acceptance fraction: 0.7133                                            emcee_sampler.py:145

22:49:08 INFO      fit restored to maximum of posterior                                         sampler_base.py:168

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

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.56 -0.13 +0.11) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.989 -0.030 +0.04) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.447057
total	-4.447057

Values of statistical measures:

	statistical measures
AIC	13.599996
BIC	14.885579
DIC	13.098420
PDIC	2.094977

[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:186

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

 no. of live points =  400
 dimensionality =    2
 *****************************************************

 MultiNest Warning!
 Parameter            1  of mode            6  is converging towards the edge of the prior.

 MultiNest Warning!
 Parameter            1  of mode            6  is converging towards the edge of the prior.
  analysing data from chains/fit-.txt ln(ev)=  -14.505665349415551      +/-  0.15044379120258144
 Total Likelihood Evaluations:         6163
 Sampling finished. Exiting MultiNest

22:49:09 INFO      fit restored to maximum of posterior                                         sampler_base.py:168

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

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.55 -0.10 +0.13) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.989 -0.030 +0.029) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.4515
total	-4.4515

Values of statistical measures:

	statistical measures
AIC	13.608882
BIC	14.894465
DIC	12.546291
PDIC	1.825818
log(Z)	-6.299730

         INFO      deleting the chain directory chains                                     multinest_sampler.py:234

[6]:

../_images/notebooks_sampler_docs_9_11.png

../_images/notebooks_sampler_docs_9_12.png

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

         INFO      sampler set to dynesty_nested                                           bayesian_analysis.py:186

4382it [00:04, 1054.08it/s, +400 | bound: 12 | nc: 1 | ncall: 19081 | eff(%): 25.598 | loglstar:   -inf < -4.492 <    inf | logz: -14.404 +/-  0.149 | dlogz:  0.001 >  0.409]

22:49:13 INFO      fit restored to maximum of posterior                                         sampler_base.py:168

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

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.57 +/- 0.12) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.990 -0.04 +0.031) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.447886
total	-4.447886

Values of statistical measures:

	statistical measures
AIC	13.601655
BIC	14.887237
DIC	13.008095
PDIC	2.056081
log(Z)	-6.255759

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

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

22:49:14 INFO      sampler set to dynesty_dynamic                                          bayesian_analysis.py:186

16715it [00:15, 1092.92it/s, batch: 8 | bound: 5 | nc: 1 | ncall: 38699 | eff(%): 43.135 | loglstar: -9.326 < -4.492 < -4.781 | logz: -14.463 +/-  0.077 | stop:  0.895]

22:49:29 INFO      fit restored to maximum of posterior                                         sampler_base.py:168

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

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.56 +/- 0.12) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.989 -0.031 +0.032) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.447062
total	-4.447062

Values of statistical measures:

	statistical measures
AIC	13.600006
BIC	14.885588
DIC	12.849871
PDIC	1.977755
log(Z)	-6.282101

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

22:49:30 INFO      sampler set to zeus                                                     bayesian_analysis.py:186

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:04<00:00, 146.62it/s]

22:49:35 INFO      fit restored to maximum of posterior                                         sampler_base.py:168

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

Summary
-------
Number of Generations: 625
Number of Parameters: 2
Number of Walkers: 20
Number of Tuning Generations: 23
Scale Factor: 1.653765
Mean Integrated Autocorrelation Time: 2.67
Effective Sample Size: 4676.96
Number of Log Probability Evaluations: 63969
Effective Samples per Log Probability Evaluation: 0.073113

None

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.56 -0.11 +0.13) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.990 -0.032 +0.034) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.44732
total	-4.44732

Values of statistical measures:

	statistical measures
AIC	13.600523
BIC	14.886105
DIC	12.965771
PDIC	2.034743

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

         INFO      sampler set to ultranest                                                bayesian_analysis.py:186

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

[ultranest] Explored until L=-4
[ultranest] Likelihood function evaluations: 8197
[ultranest]   logZ = -14.48 +- 0.1049
[ultranest] Effective samples strategy satisfied (ESS = 975.6, 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.42, need <0.5)
[ultranest]   logZ error budget: single: 0.15 bs:0.10 tail:0.41 total:0.42 required:<0.50
[ultranest] done iterating.

22:49:39 INFO      fit restored to maximum of posterior                                         sampler_base.py:168

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

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.56 -0.12 +0.13) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.990 -0.031 +0.033) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-4.448147
total	-4.448147

Values of statistical measures:

	statistical measures
AIC	13.602176
BIC	14.887758
DIC	12.971236
PDIC	2.038159
log(Z)	-6.265935

[10]:

../_images/notebooks_sampler_docs_16_12.png

../_images/notebooks_sampler_docs_16_13.png

../_images/notebooks_sampler_docs_16_14.png