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	(9.96 -0.22 +0.21) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(1.002 +/- 0.004) x 10^-1	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-9.550171
total	-9.550171

Values of statistical measures:

	statistical measures
AIC	23.806225
BIC	25.091807
DIC	23.164881
PDIC	2.030415

[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
 *****************************************************
 Starting MultiNest
 generating live points
 live points generated, starting sampling
Acceptance Rate:                        0.986842
Replacements:                                450
Total Samples:                               456
Nested Sampling ln(Z):              -8497.967033
Acceptance Rate:                        0.974659
Replacements:                                500
Total Samples:                               513
Nested Sampling ln(Z):              -6998.517450
Acceptance Rate:                        0.943396
Replacements:                                550
Total Samples:                               583
Nested Sampling ln(Z):              -6081.429384
Acceptance Rate:                        0.924499
Replacements:                                600
Total Samples:                               649
Nested Sampling ln(Z):              -5282.318752
Acceptance Rate:                        0.890411
Replacements:                                650
Total Samples:                               730
Nested Sampling ln(Z):              -4784.625593
Acceptance Rate:                        0.854701
Replacements:                                700
Total Samples:                               819
Nested Sampling ln(Z):              -4344.172901
Acceptance Rate:                        0.847458
Replacements:                                750
Total Samples:                               885
Nested Sampling ln(Z):              -4056.402620
Acceptance Rate:                        0.824742
Replacements:                                800
Total Samples:                               970
Nested Sampling ln(Z):              -3740.949058
Acceptance Rate:                        0.822846
Replacements:                                850
Total Samples:                              1033
Nested Sampling ln(Z):              -3359.761246
Acceptance Rate:                        0.816042
Replacements:                                936
Total Samples:                              1147
Nested Sampling ln(Z):              -2790.526259
Acceptance Rate:                        0.797896
Replacements:                               1062
Total Samples:                              1331
Nested Sampling ln(Z):              -2164.335849
Acceptance Rate:                        0.782946
Replacements:                               1212
Total Samples:                              1548
Nested Sampling ln(Z):              -1174.847150
Acceptance Rate:                        0.778286
Replacements:                               1362
Total Samples:                              1750
Nested Sampling ln(Z):               -717.331437
Acceptance Rate:                        0.778980
Replacements:                               1512
Total Samples:                              1941
Nested Sampling ln(Z):               -500.605805
Acceptance Rate:                        0.776636
Replacements:                               1662
Total Samples:                              2140
Nested Sampling ln(Z):               -379.966696
Acceptance Rate:                        0.775353
Replacements:                               1812
Total Samples:                              2337
Nested Sampling ln(Z):               -273.980076
Acceptance Rate:                        0.773964
Replacements:                               1962
Total Samples:                              2535
Nested Sampling ln(Z):               -209.081201
Acceptance Rate:                        0.773060
Replacements:                               2112
Total Samples:                              2732
Nested Sampling ln(Z):               -158.745077
Acceptance Rate:                        0.767821
Replacements:                               2262
Total Samples:                              2946
Nested Sampling ln(Z):               -123.659912
Acceptance Rate:                        0.766444
Replacements:                               2412
Total Samples:                              3147
Nested Sampling ln(Z):                -94.811776
Acceptance Rate:                        0.769138
Replacements:                               2562
Total Samples:                              3331
Nested Sampling ln(Z):                -72.917087
Acceptance Rate:                        0.768707
Replacements:                               2712
Total Samples:                              3528
Nested Sampling ln(Z):                -60.083382
Acceptance Rate:                        0.768322
Replacements:                               2862
Total Samples:                              3725
Nested Sampling ln(Z):                -50.121220
Acceptance Rate:                        0.771516
Replacements:                               3012
Total Samples:                              3904
Nested Sampling ln(Z):                -41.744533
Acceptance Rate:                        0.771973
Replacements:                               3162
Total Samples:                              4096
Nested Sampling ln(Z):                -36.173174
Acceptance Rate:                        0.772388
Replacements:                               3312
Total Samples:                              4288
Nested Sampling ln(Z):                -31.441577
Acceptance Rate:                        0.771906
Replacements:                               3462
Total Samples:                              4485
Nested Sampling ln(Z):                -28.438677
Acceptance Rate:                        0.769001
Replacements:                               3612
Total Samples:                              4697
Nested Sampling ln(Z):                -25.776219
Acceptance Rate:                        0.763767
Replacements:                               3731
Total Samples:                              4885
Nested Sampling ln(Z):                -23.709446
Acceptance Rate:                        0.759215
Replacements:                               3831
Total Samples:                              5046
Nested Sampling ln(Z):                -22.248807
Acceptance Rate:                        0.755090
Replacements:                               3931
Total Samples:                              5206
Nested Sampling ln(Z):                -21.379384
Acceptance Rate:                        0.751557
Replacements:                               3981
Total Samples:                              5297
Nested Sampling ln(Z):                -20.669919
Acceptance Rate:                        0.750233
Replacements:                               4031
Total Samples:                              5373
Nested Sampling ln(Z):                -20.177505
Acceptance Rate:                        0.747710
Replacements:                               4081
Total Samples:                              5458
Nested Sampling ln(Z):                -19.839960
Acceptance Rate:                        0.746881
Replacements:                               4131
Total Samples:                              5531
Nested Sampling ln(Z):                -19.601989
Acceptance Rate:                        0.747408
Replacements:                               4181
Total Samples:                              5594
Nested Sampling ln(Z):                -19.430100
Acceptance Rate:                        0.748320
Replacements:                               4231
Total Samples:                              5654
Nested Sampling ln(Z):                -19.301938
Acceptance Rate:                        0.747713
Replacements:                               4250
Total Samples:                              5684
Nested Sampling ln(Z):                -19.260323
  analysing data from chains/fit-.txt ln(ev)=  -18.913174767983953      +/-  0.14451907981360401
 Total Likelihood Evaluations:         5684
 Sampling finished. Exiting MultiNest

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.96 -0.24 +0.23) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(1.002 -0.004 +0.005) x 10^-1	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-9.551042
total	-9.551042

Values of statistical measures:

	statistical measures
AIC	23.807967
BIC	25.093549
DIC	23.222957
PDIC	2.057475
log(Z)	-8.213887

[6]:

../_images/notebooks_sampler_docs_9_7.png

../_images/notebooks_sampler_docs_9_8.png

../_images/notebooks_sampler_docs_9_9.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%|█████████▉| 4145/4146 [00:04<00:00, 935.75it/s, +400 | bound: 12 | nc: 1 | ncall: 19050 | eff(%): 24.370 | loglstar: -9.556 | logz: -18.869 +/-  0.144 | dlogz:  0.001 >  0.409]

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.95 -0.22 +0.23) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(1.002 -0.004 +0.005) x 10^-1	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-9.551415
total	-9.551415

Values of statistical measures:

	statistical measures
AIC	23.808712
BIC	25.094294
DIC	23.160226
PDIC	2.023626
log(Z)	-8.194609

[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%|█████████▉| 16258/16260 [00:15<00:00, 1029.00it/s, batch: 8 | bound: 5 | nc: 1 | ncall: 38460 | eff(%): 42.098 | loglstar: -14.310 < -9.555 < -9.904 | logz: -19.242 +/-  0.077 | stop:  0.906]

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.96 +/- 0.22) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(1.002 -0.004 +0.005) x 10^-1	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-9.550273
total	-9.550273

Values of statistical measures:

	statistical measures
AIC	23.806428
BIC	25.092010
DIC	23.138716
PDIC	2.018062
log(Z)	-8.353248

[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:04<00:00, 150.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: 20
Scale Factor: 1.139967
Mean Integrated Autocorrelation Time: 2.77
Effective Sample Size: 4512.94
Number of Log Probability Evaluations: 66250
Effective Samples per Log Probability Evaluation: 0.06812

None

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.96 -0.22 +0.21) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(1.002 -0.004 +0.005) x 10^-1	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-9.550173
total	-9.550173

Values of statistical measures:

	statistical measures
AIC	23.806228
BIC	25.091810
DIC	22.990262
PDIC	1.942127

[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: 10191
[ultranest]   logZ = -19.17 +- 0.1183
[ultranest] Effective samples strategy satisfied (ESS = 974.9, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.47+-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.15 bs:0.12 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	(9.96 -0.21 +0.22) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(1.002 -0.005 +0.004) x 10^-1	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-9.551328
total	-9.551328

Values of statistical measures:

	statistical measures
AIC	23.808539
BIC	25.094121
DIC	23.004872
PDIC	1.952076
log(Z)	-8.322507

[10]:

../_images/notebooks_sampler_docs_16_11.png

../_images/notebooks_sampler_docs_16_12.png

../_images/notebooks_sampler_docs_16_13.png

nautilus

[11]:

bayes_analysis.set_sampler("nautilus")
bayes_analysis.sampler.setup(n_live=3000, verbose=False)
bayes_analysis.sample()

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

sampler set to [blue]nautilus[/blue]

Starting the nautilus sampler...
Please report issues at github.com/johannesulf/nautilus.
Status    | Bounds | Ellipses | Networks | Calls    | f_live | N_eff | log Z
Finished  | 21     | 1        | 4        | 63600    | N/A    | 16358 | -19.10

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	(9.96 +/- 0.22) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(1.002 -0.004 +0.005) x 10^-1	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-9.550175
total	-9.550175

Values of statistical measures:

	statistical measures
AIC	23.806232
BIC	25.091814
DIC	23.111729
PDIC	2.003807
log(Z)	-19.097414

[11]:

../_images/notebooks_sampler_docs_18_9.png

../_images/notebooks_sampler_docs_18_10.png

../_images/notebooks_sampler_docs_18_11.png

autoemcee

[12]:

bayes_analysis.set_sampler("autoemcee")
bayes_analysis.sampler.setup(
    num_global_samples=10000,
    num_chains=4,
    num_walkers=None,
    max_ncalls=1000000,
    max_improvement_loops=4,
    num_initial_steps=100,
    min_autocorr_times=0,
)
bayes_analysis.sample()

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

sampler set to [blue]autoemcee[/blue]

[autoemcee] finding starting points and running initial 100 MCMC steps

finding starting points and running initial 100 MCMC steps
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 278.82it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 278.72it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 277.37it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 278.62it/s]
checking convergence (iteration 1) ...
acceptance rates: [47 48 53 56 56 57 57 57]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [53 57 58 58 59 59 59 60]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [52 53 56 58 58 59 59 60]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [50 54 54 55 56 56 57 58]% (worst few)
autocorrelation length: tau=inf -> 0x lengths

[autoemcee] rhat chain diagnostic: [1.04062013 1.03553537] (<1.010 is good)

rhat chain diagnostic: [1.04062013 1.03553537] (<1.010 is good)

[autoemcee] not converged yet at iteration 1 after 80400 evals

not converged yet at iteration 1 after 80400 evals
expected memory usage: 0.00 GiB

[autoemcee] Running 1000 MCMC steps ...

Running 1000 MCMC steps ...

[autoemcee] Starting points chosen: {np.int64(53)}, L=-12.0

Starting points chosen: {np.int64(53)}, L=-12.0

[autoemcee] Starting at [0.66168803 0.19960164] +- [8.20625154e-05 3.84342328e-06]

Starting at [0.66168803 0.19960164] +- [8.20625154e-05 3.84342328e-06]
100%|██████████| 100/100 [00:00<00:00, 275.44it/s]
100%|██████████| 1000/1000 [00:03<00:00, 280.21it/s]

[autoemcee] Starting points chosen: {np.int64(80)}, L=-12.0


Starting points chosen: {np.int64(80)}, L=-12.0

[autoemcee] Starting at [0.62864264 0.19948287] +- [7.80487281e-05 3.93888139e-06]

Starting at [0.62864264 0.19948287] +- [7.80487281e-05 3.93888139e-06]
100%|██████████| 100/100 [00:00<00:00, 279.79it/s]
100%|██████████| 1000/1000 [00:03<00:00, 280.28it/s]

[autoemcee] Starting points chosen: {np.int64(18)}, L=-12.0


Starting points chosen: {np.int64(18)}, L=-12.0

[autoemcee] Starting at [0.63672612 0.19968858] +- [7.95111708e-05 4.27892080e-06]

Starting at [0.63672612 0.19968858] +- [7.95111708e-05 4.27892080e-06]
100%|██████████| 100/100 [00:00<00:00, 278.75it/s]
100%|██████████| 1000/1000 [00:03<00:00, 278.56it/s]

[autoemcee] Starting points chosen: {np.int64(94)}, L=-12.0


Starting points chosen: {np.int64(94)}, L=-12.0

[autoemcee] Starting at [0.62269512 0.20022483] +- [9.53931288e-05 3.84857564e-06]

Starting at [0.62269512 0.20022483] +- [9.53931288e-05 3.84857564e-06]
100%|██████████| 100/100 [00:00<00:00, 276.19it/s]
100%|██████████| 1000/1000 [00:03<00:00, 278.18it/s]

[autoemcee] Used 440800 calls in last MCMC run


Used 440800 calls in last MCMC run
checking convergence (iteration 2) ...
acceptance rates: [68 68 68 68 68 68 68 69]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [66 67 68 68 68 69 69 69]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [67 68 68 68 68 68 68 68]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [66 67 67 68 68 68 68 68]% (worst few)
autocorrelation length: tau=inf -> 0x lengths

[autoemcee] rhat chain diagnostic: [1.00023328 1.00034563] (<1.010 is good)

rhat chain diagnostic: [1.00023328 1.00034563] (<1.010 is good)

[autoemcee] converged!!!

converged!!!
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	(9.97 -0.23 +0.21) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(1.002 -0.004 +0.005) x 10^-1	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-9.550357
total	-9.550357

Values of statistical measures:

	statistical measures
AIC	23.806596
BIC	25.092179
DIC	23.136753
PDIC	2.016253

[12]:

../_images/notebooks_sampler_docs_20_37.png

../_images/notebooks_sampler_docs_20_38.png

../_images/notebooks_sampler_docs_20_39.png