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
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
WARNING RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject
WARNING RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject
WARNING RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject
[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.
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.7145
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.013 +/- 0.018 | 1 / (cm2 keV s) |
| demo.spectrum.main.Sin.f | (9.97 +/- 0.04) x 10^-2 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -6.147317 |
| total | -6.147317 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 17.000517 |
| BIC | 18.286099 |
| DIC | 16.236335 |
| PDIC | 1.969801 |
WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
[5]:
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 | 1.013 -0.017 +0.018 | 1 / (cm2 keV s) |
| demo.spectrum.main.Sin.f | (9.97 +/- 0.04) x 10^-2 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -6.150649 |
| total | -6.150649 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 17.007180 |
| BIC | 18.292762 |
| DIC | 16.289595 |
| PDIC | 1.995901 |
| log(Z) | -6.585631 |
[INFO ] deleting the chain directory chains
WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
[6]:
*****************************************************
MultiNest v3.10
Copyright Farhan Feroz & Mike Hobson
Release Jul 2015
no. of live points = 400
dimensionality = 2
*****************************************************
ln(ev)= -15.163976518284725 +/- 0.14182676163695318
Total Likelihood Evaluations: 6988
Sampling finished. Exiting MultiNest
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
4297it [00:05, 723.31it/s, +400 | bound: 10 | nc: 1 | ncall: 20057 | eff(%): 23.418 | loglstar: -inf < -6.142 < inf | logz: -15.833 +/- 0.208 | dlogz: 0.001 > 0.409]
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.014 -0.017 +0.018 | 1 / (cm2 keV s) |
| demo.spectrum.main.Sin.f | (9.97 +/- 0.04) x 10^-2 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -6.147713 |
| total | -6.147713 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 17.001308 |
| BIC | 18.286891 |
| DIC | 16.277723 |
| PDIC | 1.991272 |
| log(Z) | -6.876261 |
WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
[7]:
[8]:
bayes_analysis.set_sampler("dynesty_dynamic")
bayes_analysis.sampler.setup(n_live_points=400)
bayes_analysis.sample()
xyl.plot()
bayes_analysis.results.corner_plot()
[INFO ] sampler set to dynesty_dynamic
24043it [03:48, 105.35it/s, batch: 12 | bound: 62 | nc: 1 | ncall: 57842 | eff(%): 41.567 | loglstar: -12.907 < -6.142 < -6.481 | logz: -15.482 +/- 0.184 | stop: 0.941]
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.013 +/- 0.017 | 1 / (cm2 keV s) |
| demo.spectrum.main.Sin.f | (9.97 +/- 0.04) x 10^-2 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -6.147332 |
| total | -6.147332 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 17.000547 |
| BIC | 18.286129 |
| DIC | 16.192629 |
| PDIC | 1.948349 |
| log(Z) | -6.728869 |
WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
[8]:
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:17<00:00, 36.69it/s]
Summary
-------
Number of Generations: 625
Number of Parameters: 2
Number of Walkers: 20
Number of Tuning Generations: 30
Scale Factor: 1.205591
Mean Integrated Autocorrelation Time: 2.69
Effective Sample Size: 4648.02
Number of Log Probability Evaluations: 65223
Effective Samples per Log Probability Evaluation: 0.071263
None
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.014 +/- 0.018 | 1 / (cm2 keV s) |
| demo.spectrum.main.Sin.f | (9.97 +/- 0.04) x 10^-2 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -6.14742 |
| total | -6.14742 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 17.000722 |
| BIC | 18.286304 |
| DIC | 16.297779 |
| PDIC | 2.000503 |
WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
[9]:
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=-6
[ultranest] Likelihood function evaluations: 7238
[ultranest] logZ = -15.44 +- 0.1177
[ultranest] Effective samples strategy satisfied (ESS = 976.3, 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.31, need <0.5)
[ultranest] logZ error budget: single: 0.14 bs:0.12 tail:0.41 total:0.42 required:<0.50
[ultranest] done iterating.
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.013 +/- 0.016 | 1 / (cm2 keV s) |
| demo.spectrum.main.Sin.f | (9.97 +/- 0.04) x 10^-2 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -6.147556 |
| total | -6.147556 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 17.000994 |
| BIC | 18.286576 |
| DIC | 16.138559 |
| PDIC | 1.921339 |
| log(Z) | -6.711466 |
WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("viridis").copy()
[10]:
