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))
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 | 1.001 -0.020 +0.024 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.005 -0.008 +0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -7.533279 |
| total | -7.533279 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 19.772440 |
| BIC | 21.058022 |
| DIC | 19.175325 |
| PDIC | 2.050970 |
[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()
*****************************************************
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.993377
Replacements: 450
Total Samples: 453
Nested Sampling ln(Z): -3173.159415
Acceptance Rate: 0.976562
Replacements: 500
Total Samples: 512
Nested Sampling ln(Z): -2474.194423
Acceptance Rate: 0.966608
Replacements: 550
Total Samples: 569
Nested Sampling ln(Z): -2087.863689
Acceptance Rate: 0.938967
Replacements: 600
Total Samples: 639
Nested Sampling ln(Z): -1900.308464
Acceptance Rate: 0.897790
Replacements: 650
Total Samples: 724
Nested Sampling ln(Z): -1712.254160
Acceptance Rate: 0.866337
Replacements: 700
Total Samples: 808
Nested Sampling ln(Z): -1586.939880
Acceptance Rate: 0.869061
Replacements: 750
Total Samples: 863
Nested Sampling ln(Z): -1482.549408
Acceptance Rate: 0.867380
Replacements: 811
Total Samples: 935
Nested Sampling ln(Z): -1385.647246
Acceptance Rate: 0.850607
Replacements: 911
Total Samples: 1071
Nested Sampling ln(Z): -1182.132259
Acceptance Rate: 0.830731
Replacements: 1011
Total Samples: 1217
Nested Sampling ln(Z): -1021.618714
Acceptance Rate: 0.825409
Replacements: 1111
Total Samples: 1346
Nested Sampling ln(Z): -949.511201
Acceptance Rate: 0.813298
Replacements: 1211
Total Samples: 1489
Nested Sampling ln(Z): -915.326326
Acceptance Rate: 0.799390
Replacements: 1311
Total Samples: 1640
Nested Sampling ln(Z): -765.367513
Acceptance Rate: 0.786072
Replacements: 1411
Total Samples: 1795
Nested Sampling ln(Z): -595.921556
Acceptance Rate: 0.779268
Replacements: 1511
Total Samples: 1939
Nested Sampling ln(Z): -473.177247
Acceptance Rate: 0.773404
Replacements: 1611
Total Samples: 2083
Nested Sampling ln(Z): -368.941302
Acceptance Rate: 0.768989
Replacements: 1711
Total Samples: 2225
Nested Sampling ln(Z): -301.676716
Acceptance Rate: 0.763491
Replacements: 1811
Total Samples: 2372
Nested Sampling ln(Z): -261.134304
Acceptance Rate: 0.761962
Replacements: 1911
Total Samples: 2508
Nested Sampling ln(Z): -222.848932
Acceptance Rate: 0.753748
Replacements: 2011
Total Samples: 2668
Nested Sampling ln(Z): -201.510874
Acceptance Rate: 0.751781
Replacements: 2111
Total Samples: 2808
Nested Sampling ln(Z): -170.630839
Acceptance Rate: 0.750254
Replacements: 2211
Total Samples: 2947
Nested Sampling ln(Z): -144.817610
Acceptance Rate: 0.747171
Replacements: 2311
Total Samples: 3093
Nested Sampling ln(Z): -125.871912
Acceptance Rate: 0.743447
Replacements: 2411
Total Samples: 3243
Nested Sampling ln(Z): -107.593299
Acceptance Rate: 0.742527
Replacements: 2509
Total Samples: 3379
Nested Sampling ln(Z): -92.849763
Acceptance Rate: 0.740880
Replacements: 2559
Total Samples: 3454
Nested Sampling ln(Z): -80.591164
Acceptance Rate: 0.740142
Replacements: 2609
Total Samples: 3525
Nested Sampling ln(Z): -69.847291
Acceptance Rate: 0.739638
Replacements: 2659
Total Samples: 3595
Nested Sampling ln(Z): -63.293008
Acceptance Rate: 0.738751
Replacements: 2709
Total Samples: 3667
Nested Sampling ln(Z): -57.591972
Acceptance Rate: 0.740671
Replacements: 2759
Total Samples: 3725
Nested Sampling ln(Z): -50.554643
Acceptance Rate: 0.741161
Replacements: 2809
Total Samples: 3790
Nested Sampling ln(Z): -45.453364
Acceptance Rate: 0.742212
Replacements: 2859
Total Samples: 3852
Nested Sampling ln(Z): -41.544603
Acceptance Rate: 0.742281
Replacements: 2909
Total Samples: 3919
Nested Sampling ln(Z): -37.397570
Acceptance Rate: 0.743654
Replacements: 2959
Total Samples: 3979
Nested Sampling ln(Z): -34.421386
Acceptance Rate: 0.744802
Replacements: 3009
Total Samples: 4040
Nested Sampling ln(Z): -31.313266
Acceptance Rate: 0.744463
Replacements: 3059
Total Samples: 4109
Nested Sampling ln(Z): -28.812860
Acceptance Rate: 0.743780
Replacements: 3109
Total Samples: 4180
Nested Sampling ln(Z): -26.710340
Acceptance Rate: 0.741549
Replacements: 3159
Total Samples: 4260
Nested Sampling ln(Z): -25.327760
Acceptance Rate: 0.742137
Replacements: 3209
Total Samples: 4324
Nested Sampling ln(Z): -23.986949
Acceptance Rate: 0.743216
Replacements: 3259
Total Samples: 4385
Nested Sampling ln(Z): -22.889399
Acceptance Rate: 0.745102
Replacements: 3309
Total Samples: 4441
Nested Sampling ln(Z): -21.976097
Acceptance Rate: 0.745285
Replacements: 3359
Total Samples: 4507
Nested Sampling ln(Z): -21.055398
Acceptance Rate: 0.746442
Replacements: 3409
Total Samples: 4567
Nested Sampling ln(Z): -20.275976
Acceptance Rate: 0.746923
Replacements: 3459
Total Samples: 4631
Nested Sampling ln(Z): -19.620989
Acceptance Rate: 0.747709
Replacements: 3509
Total Samples: 4693
Nested Sampling ln(Z): -19.089467
Acceptance Rate: 0.748475
Replacements: 3559
Total Samples: 4755
Nested Sampling ln(Z): -18.659327
Acceptance Rate: 0.749844
Replacements: 3609
Total Samples: 4813
Nested Sampling ln(Z): -18.302044
Acceptance Rate: 0.751489
Replacements: 3659
Total Samples: 4869
Nested Sampling ln(Z): -17.985284
Acceptance Rate: 0.752333
Replacements: 3709
Total Samples: 4930
Nested Sampling ln(Z): -17.712361
Acceptance Rate: 0.753005
Replacements: 3759
Total Samples: 4992
Nested Sampling ln(Z): -17.481631
Acceptance Rate: 0.754257
Replacements: 3809
Total Samples: 5050
Nested Sampling ln(Z): -17.287473
Acceptance Rate: 0.755482
Replacements: 3859
Total Samples: 5108
Nested Sampling ln(Z): -17.119134
Acceptance Rate: 0.756678
Replacements: 3909
Total Samples: 5166
Nested Sampling ln(Z): -16.976598
Acceptance Rate: 0.756690
Replacements: 3959
Total Samples: 5232
Nested Sampling ln(Z): -16.855030
Acceptance Rate: 0.756701
Replacements: 4009
Total Samples: 5298
Nested Sampling ln(Z): -16.751132
Acceptance Rate: 0.757842
Replacements: 4059
Total Samples: 5356
Nested Sampling ln(Z): -16.663992
Acceptance Rate: 0.758538
Replacements: 4109
Total Samples: 5417
Nested Sampling ln(Z): -16.590895
Acceptance Rate: 0.759158
Replacements: 4145
Total Samples: 5460
Nested Sampling ln(Z): -16.545333
analysing data from chains/fit-.txt
ln(ev)= -16.190250204928983 +/- 0.13828564699353829
Total Likelihood Evaluations: 5460
Sampling finished. Exiting MultiNest
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.000 +/- 0.020 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.005 +/- 0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -7.533796 |
| total | -7.533796 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 19.773474 |
| BIC | 21.059056 |
| DIC | 18.889493 |
| PDIC | 1.911227 |
| log(Z) | -7.031336 |
[6]:
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%|█████████▉| 3878/3879 [00:01<00:00, 1951.62it/s, +400 | bound: 10 | nc: 1 | ncall: 18256 | eff(%): 23.958 | loglstar: -7.532 | logz: -16.181 +/- 0.138 | dlogz: 0.001 > 0.409]
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.002 -0.025 +0.021 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.005 -0.007 +0.008) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -7.534007 |
| total | -7.534007 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 19.773897 |
| BIC | 21.059479 |
| DIC | 19.097674 |
| PDIC | 2.013462 |
| log(Z) | -7.027256 |
[7]:
[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%|█████████▉| 15032/15035 [00:06<00:00, 2299.84it/s, batch: 8 | bound: 6 | nc: 1 | ncall: 35228 | eff(%): 42.580 | loglstar: -12.266 < -7.533 < -7.704 | logz: -16.172 +/- 0.076 | stop: 0.987]
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.001 +/- 0.022 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.005 +/- 0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -7.533281 |
| total | -7.533281 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 19.772444 |
| BIC | 21.058026 |
| DIC | 19.030898 |
| PDIC | 1.981985 |
| log(Z) | -7.017040 |
[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()
Initialising ensemble of 20 walkers...
Sampling progress : 100%|██████████| 625/625 [00:01<00:00, 363.24it/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: 38
Scale Factor: 1.579642
Mean Integrated Autocorrelation Time: 2.63
Effective Sample Size: 4753.75
Number of Log Probability Evaluations: 64323
Effective Samples per Log Probability Evaluation: 0.073904
None
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.001 -0.023 +0.022 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.005 +/- 0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -7.533209 |
| total | -7.533209 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 19.772300 |
| BIC | 21.057882 |
| DIC | 19.203812 |
| PDIC | 2.068460 |
[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()
sampler set to [blue]ultranest[/blue]
[ultranest] Sampling 400 live points from prior ...
[ultranest] Explored until L=-8
[ultranest] Likelihood function evaluations: 6941
[ultranest] logZ = -16.31 +- 0.1097
[ultranest] Effective samples strategy satisfied (ESS = 985.1, 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.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 | 1.001 -0.022 +0.021 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.005 -0.006 +0.008) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -7.534722 |
| total | -7.534722 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 19.775327 |
| BIC | 21.060909 |
| DIC | 19.111084 |
| PDIC | 2.020415 |
| log(Z) | -7.101731 |
[10]:
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 | 19 | 1 | 4 | 59500 | N/A | 16391 | -16.17
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 | 1.001 -0.023 +0.022 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.005 +/- 0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -7.533288 |
| total | -7.533288 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 19.772458 |
| BIC | 21.058040 |
| DIC | 19.078276 |
| PDIC | 2.005681 |
| log(Z) | -16.174412 |
[11]:
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, 629.56it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 637.70it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 621.34it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 430.48it/s]
checking convergence (iteration 1) ...
acceptance rates: [48 58 60 60 60 60 61 61]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [51 51 54 57 59 59 60 60]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [55 56 58 60 60 60 60 61]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [53 53 54 55 55 56 58 58]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
[autoemcee] rhat chain diagnostic: [1.01672188 1.01625579] (<1.010 is good)
rhat chain diagnostic: [1.01672188 1.01625579] (<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(92)}, L=-7.6
Starting points chosen: {np.int64(92)}, L=-7.6
[autoemcee] Starting at [0.63665937 0.19977256] +- [1.02536660e-04 7.08329971e-06]
Starting at [0.63665937 0.19977256] +- [1.02536660e-04 7.08329971e-06]
100%|██████████| 100/100 [00:00<00:00, 422.86it/s]
100%|██████████| 1000/1000 [00:01<00:00, 626.73it/s]
[autoemcee] Starting points chosen: {np.int64(18)}, L=-7.6
Starting points chosen: {np.int64(18)}, L=-7.6
[autoemcee] Starting at [0.64009971 0.2003788 ] +- [8.54920301e-05 6.98238317e-06]
Starting at [0.64009971 0.2003788 ] +- [8.54920301e-05 6.98238317e-06]
100%|██████████| 100/100 [00:00<00:00, 534.32it/s]
100%|██████████| 1000/1000 [00:01<00:00, 627.18it/s]
[autoemcee] Starting points chosen: {np.int64(83)}, L=-7.6
Starting points chosen: {np.int64(83)}, L=-7.6
[autoemcee] Starting at [0.62408015 0.20276923] +- [1.03194625e-04 6.07994787e-06]
Starting at [0.62408015 0.20276923] +- [1.03194625e-04 6.07994787e-06]
100%|██████████| 100/100 [00:00<00:00, 626.76it/s]
100%|██████████| 1000/1000 [00:01<00:00, 628.39it/s]
[autoemcee] Starting points chosen: {np.int64(13)}, L=-7.6
Starting points chosen: {np.int64(13)}, L=-7.6
[autoemcee] Starting at [0.60624429 0.20276385] +- [9.36711193e-05 6.95752786e-06]
Starting at [0.60624429 0.20276385] +- [9.36711193e-05 6.95752786e-06]
100%|██████████| 100/100 [00:00<00:00, 631.06it/s]
100%|██████████| 1000/1000 [00:01<00:00, 621.01it/s]
[autoemcee] Used 440800 calls in last MCMC run
Used 440800 calls in last MCMC run
checking convergence (iteration 2) ...
acceptance rates: [67 68 68 68 68 68 69 69]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [66 67 67 67 68 68 68 68]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [67 67 67 67 68 68 68 68]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [67 67 67 67 68 68 68 68]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
[autoemcee] rhat chain diagnostic: [1.0004135 1.00014376] (<1.010 is good)
rhat chain diagnostic: [1.0004135 1.00014376] (<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 | 1.001 -0.023 +0.022 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.005 +/- 0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -7.533443 |
| total | -7.533443 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 19.772767 |
| BIC | 21.058350 |
| DIC | 19.088666 |
| PDIC | 2.010908 |
[12]: