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.85 -0.16 +0.17) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.969 -0.028 +0.026) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.217231
total	-10.217231

Values of statistical measures:

	statistical measures
AIC	25.140344
BIC	26.425926
DIC	24.487552
PDIC	2.025729

[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.984683
Replacements:                                450
Total Samples:                               457
Nested Sampling ln(Z):             -11606.820042
Acceptance Rate:                        0.972763
Replacements:                                500
Total Samples:                               514
Nested Sampling ln(Z):              -8226.023072
Acceptance Rate:                        0.945017
Replacements:                                550
Total Samples:                               582
Nested Sampling ln(Z):              -6498.624097
Acceptance Rate:                        0.918836
Replacements:                                600
Total Samples:                               653
Nested Sampling ln(Z):              -5224.030897
Acceptance Rate:                        0.885559
Replacements:                                650
Total Samples:                               734
Nested Sampling ln(Z):              -4358.821334
Acceptance Rate:                        0.846433
Replacements:                                700
Total Samples:                               827
Nested Sampling ln(Z):              -3727.351460
Acceptance Rate:                        0.786988
Replacements:                                750
Total Samples:                               953
Nested Sampling ln(Z):              -3206.863717
Acceptance Rate:                        0.756144
Replacements:                                800
Total Samples:                              1058
Nested Sampling ln(Z):              -2751.995920
Acceptance Rate:                        0.750221
Replacements:                                850
Total Samples:                              1133
Nested Sampling ln(Z):              -2449.562432
Acceptance Rate:                        0.745033
Replacements:                                900
Total Samples:                              1208
Nested Sampling ln(Z):              -2222.159729
Acceptance Rate:                        0.747443
Replacements:                                950
Total Samples:                              1271
Nested Sampling ln(Z):              -2064.434959

 MultiNest Warning!
 Parameter            2  of mode            3  is converging towards the edge of the prior.
Acceptance Rate:                        0.747573
Replacements:                               1001
Total Samples:                              1339
Nested Sampling ln(Z):              -1898.195562

 MultiNest Warning!
 Parameter            2  of mode            3  is converging towards the edge of the prior.
Acceptance Rate:                        0.746441
Replacements:                               1101
Total Samples:                              1475
Nested Sampling ln(Z):              -1704.584570
Acceptance Rate:                        0.749198
Replacements:                               1168
Total Samples:                              1559
Nested Sampling ln(Z):              -1696.323670
Acceptance Rate:                        0.748771
Replacements:                               1219
Total Samples:                              1628
Nested Sampling ln(Z):              -1584.185906
Acceptance Rate:                        0.753232
Replacements:                               1340
Total Samples:                              1779
Nested Sampling ln(Z):               -951.149846
Acceptance Rate:                        0.761759
Replacements:                               1490
Total Samples:                              1956
Nested Sampling ln(Z):               -672.440522
Acceptance Rate:                        0.764212
Replacements:                               1640
Total Samples:                              2146
Nested Sampling ln(Z):               -483.759708
Acceptance Rate:                        0.764657
Replacements:                               1826
Total Samples:                              2388
Nested Sampling ln(Z):               -393.860237
Acceptance Rate:                        0.759494
Replacements:                               1980
Total Samples:                              2607
Nested Sampling ln(Z):               -291.371032
Acceptance Rate:                        0.756124
Replacements:                               2130
Total Samples:                              2817
Nested Sampling ln(Z):               -226.734439
Acceptance Rate:                        0.752972
Replacements:                               2280
Total Samples:                              3028
Nested Sampling ln(Z):               -146.734040
Acceptance Rate:                        0.746544
Replacements:                               2430
Total Samples:                              3255
Nested Sampling ln(Z):               -113.620172
Acceptance Rate:                        0.737800
Replacements:                               2555
Total Samples:                              3463
Nested Sampling ln(Z):                -90.912073
Acceptance Rate:                        0.733628
Replacements:                               2655
Total Samples:                              3619
Nested Sampling ln(Z):                -72.157376
Acceptance Rate:                        0.730769
Replacements:                               2755
Total Samples:                              3770
Nested Sampling ln(Z):                -58.252620
Acceptance Rate:                        0.724987
Replacements:                               2855
Total Samples:                              3938
Nested Sampling ln(Z):                -48.007014
Acceptance Rate:                        0.720556
Replacements:                               2955
Total Samples:                              4101
Nested Sampling ln(Z):                -40.919668
Acceptance Rate:                        0.716463
Replacements:                               3055
Total Samples:                              4264
Nested Sampling ln(Z):                -34.968582
Acceptance Rate:                        0.713962
Replacements:                               3155
Total Samples:                              4419
Nested Sampling ln(Z):                -31.300802
Acceptance Rate:                        0.712878
Replacements:                               3255
Total Samples:                              4566
Nested Sampling ln(Z):                -28.434121
Acceptance Rate:                        0.710956
Replacements:                               3355
Total Samples:                              4719
Nested Sampling ln(Z):                -26.224727
Acceptance Rate:                        0.705966
Replacements:                               3455
Total Samples:                              4894
Nested Sampling ln(Z):                -24.503092
Acceptance Rate:                        0.704938
Replacements:                               3555
Total Samples:                              5043
Nested Sampling ln(Z):                -23.375323
Acceptance Rate:                        0.702479
Replacements:                               3655
Total Samples:                              5203
Nested Sampling ln(Z):                -22.530961
Acceptance Rate:                        0.700821
Replacements:                               3755
Total Samples:                              5358
Nested Sampling ln(Z):                -21.866954
Acceptance Rate:                        0.696477
Replacements:                               3855
Total Samples:                              5535
Nested Sampling ln(Z):                -21.355088
Acceptance Rate:                        0.693981
Replacements:                               3955
Total Samples:                              5699
Nested Sampling ln(Z):                -20.966384
Acceptance Rate:                        0.690612
Replacements:                               4009
Total Samples:                              5805
Nested Sampling ln(Z):                -20.679999
Acceptance Rate:                        0.685874
Replacements:                               4059
Total Samples:                              5918
Nested Sampling ln(Z):                -20.464404
Acceptance Rate:                        0.683580
Replacements:                               4109
Total Samples:                              6011
Nested Sampling ln(Z):                -20.307994
Acceptance Rate:                        0.681221
Replacements:                               4150
Total Samples:                              6092
Nested Sampling ln(Z):                -20.211928
  analysing data from chains/fit-.txt
 ln(ev)=  -19.858843975459777      +/-  0.14700605119366805
 Total Likelihood Evaluations:         6092
 Sampling finished. Exiting MultiNest

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.85 -0.15 +0.17) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.971 -0.028 +0.025) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.219245
total	-10.219245

Values of statistical measures:

	statistical measures
AIC	25.144373
BIC	26.429955
DIC	24.305536
PDIC	1.934353
log(Z)	-8.624586

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%|█████████▉| 4489/4490 [00:04<00:00, 1074.11it/s, +400 | bound: 14 | nc: 1 | ncall: 19592 | eff(%): 25.474 | loglstar: -10.232 | logz: -20.404 +/-  0.152 | dlogz:  0.001 >  0.409]

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.85 +/- 0.16) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.969 -0.028 +0.029) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.217342
total	-10.217342

Values of statistical measures:

	statistical measures
AIC	25.140566
BIC	26.426148
DIC	24.597402
PDIC	2.081255
log(Z)	-8.861378

[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%|█████████▉| 17218/17221 [00:13<00:00, 1245.46it/s, batch: 8 | bound: 5 | nc: 1 | ncall: 40415 | eff(%): 42.509 | loglstar: -15.196 < -10.233 < -10.448 | logz: -20.597 +/-  0.078 | stop:  0.879]

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.85 +/- 0.16) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.969 -0.027 +0.028) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.217262
total	-10.217262

Values of statistical measures:

	statistical measures
AIC	25.140406
BIC	26.425988
DIC	24.466930
PDIC	2.015762
log(Z)	-8.943592

[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:02<00:00, 210.63it/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: 24
Scale Factor: 1.177629
Mean Integrated Autocorrelation Time: 2.73
Effective Sample Size: 4578.43
Number of Log Probability Evaluations: 65104
Effective Samples per Log Probability Evaluation: 0.070325

None

Maximum a posteriori probability (MAP) point:

	result	unit
parameter
demo.spectrum.main.Sin.K	(9.85 -0.16 +0.15) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.969 -0.027 +0.028) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.217238
total	-10.217238

Values of statistical measures:

	statistical measures
AIC	25.140358
BIC	26.425940
DIC	24.448633
PDIC	2.006792

[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: 14680
[ultranest]   logZ = -20.1 +- 0.1036
[ultranest] Effective samples strategy satisfied (ESS = 966.4, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.45+-0.05 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.

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.84 -0.14 +0.16) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.970 -0.026 +0.027) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.220217
total	-10.220217

Values of statistical measures:

	statistical measures
AIC	25.146317
BIC	26.431899
DIC	24.236155
PDIC	1.900334
log(Z)	-8.727218

[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        | 66100    | N/A    | 16641 | -20.34

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

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.217244
total	-10.217244

Values of statistical measures:

	statistical measures
AIC	25.140371
BIC	26.425953
DIC	24.426432
PDIC	1.995801
log(Z)	-20.337102

[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, 396.26it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 383.84it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 397.76it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 391.38it/s]
checking convergence (iteration 1) ...
acceptance rates: [51 54 56 57 57 57 58 59]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [46 50 51 51 52 53 53 54]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [51 52 53 54 54 55 56 56]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [49 50 53 54 54 54 55 55]% (worst few)
autocorrelation length: tau=inf -> 0x lengths

[autoemcee] rhat chain diagnostic: [1.04493995 1.05166901] (<1.010 is good)

rhat chain diagnostic: [1.04493995 1.05166901] (<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(34)}, L=-11.6

Starting points chosen: {np.int64(34)}, L=-11.6

[autoemcee] Starting at [0.60222867 0.19949939] +- [6.91184026e-05 2.33190549e-06]

Starting at [0.60222867 0.19949939] +- [6.91184026e-05 2.33190549e-06]
100%|██████████| 100/100 [00:00<00:00, 400.44it/s]
100%|██████████| 1000/1000 [00:02<00:00, 390.91it/s]

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


Starting points chosen: {np.int64(20)}, L=-11.6

[autoemcee] Starting at [0.61469362 0.19943942] +- [7.30915968e-05 2.45700901e-06]

Starting at [0.61469362 0.19943942] +- [7.30915968e-05 2.45700901e-06]
100%|██████████| 100/100 [00:00<00:00, 389.90it/s]
100%|██████████| 1000/1000 [00:02<00:00, 394.67it/s]

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


Starting points chosen: {np.int64(55)}, L=-11.6

[autoemcee] Starting at [0.63116213 0.19888821] +- [7.5048799e-05 2.6561787e-06]

Starting at [0.63116213 0.19888821] +- [7.5048799e-05 2.6561787e-06]
100%|██████████| 100/100 [00:00<00:00, 391.13it/s]
100%|██████████| 1000/1000 [00:02<00:00, 388.47it/s]

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


Starting points chosen: {np.int64(62)}, L=-11.6

[autoemcee] Starting at [0.62182179 0.19920611] +- [7.61168769e-05 2.94366722e-06]

Starting at [0.62182179 0.19920611] +- [7.61168769e-05 2.94366722e-06]
100%|██████████| 100/100 [00:00<00:00, 390.60it/s]
100%|██████████| 1000/1000 [00:02<00:00, 397.33it/s]

[autoemcee] Used 440800 calls in last MCMC run


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

[autoemcee] rhat chain diagnostic: [1.00013631 1.00026119] (<1.010 is good)

rhat chain diagnostic: [1.00013631 1.00026119] (<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.85 +/- 0.16) x 10^-1	1 / (keV s cm2)
demo.spectrum.main.Sin.f	(9.969 -0.026 +0.028) x 10^-2	rad / keV

Values of -log(posterior) at the minimum:

	-log(posterior)
demo	-10.217306
total	-10.217306

Values of statistical measures:

	statistical measures
AIC	25.140495
BIC	26.426077
DIC	24.540023
PDIC	2.052329

[12]:

../_images/notebooks_sampler_docs_20_37.png

../_images/notebooks_sampler_docs_20_38.png

../_images/notebooks_sampler_docs_20_39.png