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.015 -0.012 +0.013 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.002 -0.008 +0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -11.257021 |
| total | -11.257021 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 27.219924 |
| BIC | 28.505506 |
| DIC | 26.383511 |
| PDIC | 1.926752 |
[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.995575
Replacements: 450
Total Samples: 452
Nested Sampling ln(Z): -10414.472465
Acceptance Rate: 0.984252
Replacements: 500
Total Samples: 508
Nested Sampling ln(Z): -8056.700337
Acceptance Rate: 0.948276
Replacements: 550
Total Samples: 580
Nested Sampling ln(Z): -6390.878328
Acceptance Rate: 0.907716
Replacements: 600
Total Samples: 661
Nested Sampling ln(Z): -4935.509639
Acceptance Rate: 0.866667
Replacements: 650
Total Samples: 750
Nested Sampling ln(Z): -3999.070520
Acceptance Rate: 0.852619
Replacements: 700
Total Samples: 821
Nested Sampling ln(Z): -3121.574799
Acceptance Rate: 0.834260
Replacements: 750
Total Samples: 899
Nested Sampling ln(Z): -2562.771690
Acceptance Rate: 0.836820
Replacements: 800
Total Samples: 956
Nested Sampling ln(Z): -2170.078102
Acceptance Rate: 0.834971
Replacements: 850
Total Samples: 1018
Nested Sampling ln(Z): -1918.861065
Acceptance Rate: 0.831793
Replacements: 900
Total Samples: 1082
Nested Sampling ln(Z): -1741.364210
Acceptance Rate: 0.814751
Replacements: 950
Total Samples: 1166
Nested Sampling ln(Z): -1626.135296
Acceptance Rate: 0.808407
Replacements: 1000
Total Samples: 1237
Nested Sampling ln(Z): -1510.282092
Acceptance Rate: 0.807071
Replacements: 1050
Total Samples: 1301
Nested Sampling ln(Z): -1442.241843
Acceptance Rate: 0.805271
Replacements: 1100
Total Samples: 1366
Nested Sampling ln(Z): -1368.727162
Acceptance Rate: 0.803634
Replacements: 1150
Total Samples: 1431
Nested Sampling ln(Z): -1295.471504
Acceptance Rate: 0.796284
Replacements: 1200
Total Samples: 1507
Nested Sampling ln(Z): -1242.695425
Acceptance Rate: 0.794155
Replacements: 1250
Total Samples: 1574
Nested Sampling ln(Z): -1181.973733
Acceptance Rate: 0.790390
Replacements: 1316
Total Samples: 1665
Nested Sampling ln(Z): -1071.323797
Acceptance Rate: 0.784488
Replacements: 1416
Total Samples: 1805
Nested Sampling ln(Z): -897.019468
Acceptance Rate: 0.778500
Replacements: 1557
Total Samples: 2000
Nested Sampling ln(Z): -815.803559
Acceptance Rate: 0.768227
Replacements: 1707
Total Samples: 2222
Nested Sampling ln(Z): -687.596553
Acceptance Rate: 0.763883
Replacements: 1857
Total Samples: 2431
Nested Sampling ln(Z): -532.031788
Acceptance Rate: 0.755648
Replacements: 2007
Total Samples: 2656
Nested Sampling ln(Z): -417.584926
Acceptance Rate: 0.751568
Replacements: 2157
Total Samples: 2870
Nested Sampling ln(Z): -335.144200
Acceptance Rate: 0.752233
Replacements: 2274
Total Samples: 3023
Nested Sampling ln(Z): -273.896178
Acceptance Rate: 0.750079
Replacements: 2374
Total Samples: 3165
Nested Sampling ln(Z): -233.106096
Acceptance Rate: 0.746530
Replacements: 2474
Total Samples: 3314
Nested Sampling ln(Z): -195.466029
Acceptance Rate: 0.744576
Replacements: 2574
Total Samples: 3457
Nested Sampling ln(Z): -162.604405
Acceptance Rate: 0.738470
Replacements: 2674
Total Samples: 3621
Nested Sampling ln(Z): -130.963743
Acceptance Rate: 0.734034
Replacements: 2724
Total Samples: 3711
Nested Sampling ln(Z): -110.417567
Acceptance Rate: 0.728466
Replacements: 2774
Total Samples: 3808
Nested Sampling ln(Z): -97.265611
Acceptance Rate: 0.727460
Replacements: 2824
Total Samples: 3882
Nested Sampling ln(Z): -83.564293
Acceptance Rate: 0.727227
Replacements: 2874
Total Samples: 3952
Nested Sampling ln(Z): -72.122321
Acceptance Rate: 0.723762
Replacements: 2924
Total Samples: 4040
Nested Sampling ln(Z): -60.305608
Acceptance Rate: 0.719923
Replacements: 2974
Total Samples: 4131
Nested Sampling ln(Z): -52.651646
Acceptance Rate: 0.716079
Replacements: 3024
Total Samples: 4223
Nested Sampling ln(Z): -46.797414
Acceptance Rate: 0.713888
Replacements: 3074
Total Samples: 4306
Nested Sampling ln(Z): -42.491331
Acceptance Rate: 0.712916
Replacements: 3124
Total Samples: 4382
Nested Sampling ln(Z): -39.434256
Acceptance Rate: 0.706589
Replacements: 3174
Total Samples: 4492
Nested Sampling ln(Z): -36.342821
Acceptance Rate: 0.703930
Replacements: 3224
Total Samples: 4580
Nested Sampling ln(Z): -33.564350
Acceptance Rate: 0.705148
Replacements: 3274
Total Samples: 4643
Nested Sampling ln(Z): -31.016710
Acceptance Rate: 0.706182
Replacements: 3324
Total Samples: 4707
Nested Sampling ln(Z): -29.128115
Acceptance Rate: 0.707189
Replacements: 3374
Total Samples: 4771
Nested Sampling ln(Z): -27.605557
Acceptance Rate: 0.707584
Replacements: 3424
Total Samples: 4839
Nested Sampling ln(Z): -26.317484
Acceptance Rate: 0.708257
Replacements: 3474
Total Samples: 4905
Nested Sampling ln(Z): -25.328915
Acceptance Rate: 0.708484
Replacements: 3524
Total Samples: 4974
Nested Sampling ln(Z): -24.533095
Acceptance Rate: 0.708986
Replacements: 3574
Total Samples: 5041
Nested Sampling ln(Z): -23.934644
Acceptance Rate: 0.709059
Replacements: 3624
Total Samples: 5111
Nested Sampling ln(Z): -23.406288
Acceptance Rate: 0.709266
Replacements: 3674
Total Samples: 5180
Nested Sampling ln(Z): -22.933380
Acceptance Rate: 0.710010
Replacements: 3724
Total Samples: 5245
Nested Sampling ln(Z): -22.508575
Acceptance Rate: 0.711538
Replacements: 3774
Total Samples: 5304
Nested Sampling ln(Z): -22.155118
Acceptance Rate: 0.711839
Replacements: 3824
Total Samples: 5372
Nested Sampling ln(Z): -21.862750
Acceptance Rate: 0.713312
Replacements: 3874
Total Samples: 5431
Nested Sampling ln(Z): -21.623490
Acceptance Rate: 0.714754
Replacements: 3924
Total Samples: 5490
Nested Sampling ln(Z): -21.426570
Acceptance Rate: 0.714620
Replacements: 3974
Total Samples: 5561
Nested Sampling ln(Z): -21.253328
Acceptance Rate: 0.715378
Replacements: 4024
Total Samples: 5625
Nested Sampling ln(Z): -21.107860
Acceptance Rate: 0.714988
Replacements: 4074
Total Samples: 5698
Nested Sampling ln(Z): -20.986795
Acceptance Rate: 0.716097
Replacements: 4124
Total Samples: 5759
Nested Sampling ln(Z): -20.888266
Acceptance Rate: 0.716567
Replacements: 4174
Total Samples: 5825
Nested Sampling ln(Z): -20.806817
Acceptance Rate: 0.717032
Replacements: 4214
Total Samples: 5877
Nested Sampling ln(Z): -20.752418
analysing data from chains/fit-.txt
ln(ev)= -20.401990503785903 +/- 0.14290981923284479
Total Likelihood Evaluations: 5877
Sampling finished. Exiting MultiNest
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.015 -0.012 +0.014 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.003 +/- 0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -11.256572 |
| total | -11.256572 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 27.219027 |
| BIC | 28.504609 |
| DIC | 26.474342 |
| PDIC | 1.976409 |
| log(Z) | -8.860472 |
[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%|█████████▉| 4083/4084 [00:05<00:00, 700.28it/s, +400 | bound: 13 | nc: 1 | ncall: 20106 | eff(%): 22.749 | loglstar: -11.243 | logz: -20.399 +/- 0.143 | dlogz: 0.001 > 0.409]
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.015 -0.012 +0.015 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.002 -0.007 +0.008) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -11.258207 |
| total | -11.258207 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 27.222296 |
| BIC | 28.507878 |
| DIC | 26.459033 |
| PDIC | 1.965552 |
| log(Z) | -8.859311 |
[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%|█████████▉| 16378/16380 [00:20<00:00, 809.97it/s, batch: 8 | bound: 5 | nc: 1 | ncall: 40073 | eff(%): 40.595 | loglstar: -16.112 < -11.241 < -11.665 | logz: -20.474 +/- 0.074 | stop: 0.885]
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.015 -0.012 +0.014 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.003 -0.008 +0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -11.256641 |
| total | -11.256641 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 27.219165 |
| BIC | 28.504747 |
| DIC | 26.525430 |
| PDIC | 1.999048 |
| log(Z) | -8.892031 |
[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:05<00:00, 119.00it/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: 25
Scale Factor: 0.977129
Mean Integrated Autocorrelation Time: 2.77
Effective Sample Size: 4504.99
Number of Log Probability Evaluations: 67505
Effective Samples per Log Probability Evaluation: 0.066736
None
Maximum a posteriori probability (MAP) point:
| result | unit | |
|---|---|---|
| parameter | ||
| demo.spectrum.main.Sin.K | 1.015 -0.011 +0.015 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.003 +/- 0.008) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -11.256613 |
| total | -11.256613 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 27.219109 |
| BIC | 28.504691 |
| DIC | 26.550852 |
| PDIC | 2.010378 |
[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=-1e+01
[ultranest] Likelihood function evaluations: 13710
[ultranest] logZ = -20.56 +- 0.1024
[ultranest] Effective samples strategy satisfied (ESS = 978.5, need >400)
[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-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.14 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 | 1.015 -0.012 +0.015 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.003 -0.008 +0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -11.256689 |
| total | -11.256689 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 27.219261 |
| BIC | 28.504843 |
| DIC | 26.609641 |
| PDIC | 2.037707 |
| log(Z) | -8.928279 |
[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 | 20 | 1 | 4 | 62500 | N/A | 16403 | -20.41
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.015 -0.012 +0.015 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.002 +/- 0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -11.256728 |
| total | -11.256728 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 27.219338 |
| BIC | 28.504920 |
| DIC | 26.436160 |
| PDIC | 1.954825 |
| log(Z) | -20.407000 |
[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, 255.00it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 253.80it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 256.39it/s]
global sampling for starting point ...
100%|██████████| 100/100 [00:00<00:00, 259.37it/s]
checking convergence (iteration 1) ...
acceptance rates: [53 53 56 56 56 57 57 59]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [48 48 48 53 54 55 56 57]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [54 55 56 57 59 59 59 59]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [46 49 51 51 53 54 55 55]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
[autoemcee] rhat chain diagnostic: [1.03408822 1.0338745 ] (<1.010 is good)
rhat chain diagnostic: [1.03408822 1.0338745 ] (<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(82)}, L=-12.2
Starting points chosen: {np.int64(82)}, L=-12.2
[autoemcee] Starting at [0.62795592 0.19932823] +- [6.29672587e-05 6.73309814e-06]
Starting at [0.62795592 0.19932823] +- [6.29672587e-05 6.73309814e-06]
100%|██████████| 100/100 [00:00<00:00, 252.95it/s]
100%|██████████| 1000/1000 [00:03<00:00, 256.85it/s]
[autoemcee] Starting points chosen: {np.int64(0)}, L=-12.2
Starting points chosen: {np.int64(0)}, L=-12.2
[autoemcee] Starting at [0.636291 0.20079998] +- [4.75661705e-05 6.49863506e-06]
Starting at [0.636291 0.20079998] +- [4.75661705e-05 6.49863506e-06]
100%|██████████| 100/100 [00:00<00:00, 259.40it/s]
100%|██████████| 1000/1000 [00:03<00:00, 257.22it/s]
[autoemcee] Starting points chosen: {np.int64(72)}, L=-12.2
Starting points chosen: {np.int64(72)}, L=-12.2
[autoemcee] Starting at [0.63447799 0.19938494] +- [6.74443105e-05 7.23551060e-06]
Starting at [0.63447799 0.19938494] +- [6.74443105e-05 7.23551060e-06]
100%|██████████| 100/100 [00:00<00:00, 256.19it/s]
100%|██████████| 1000/1000 [00:03<00:00, 255.46it/s]
[autoemcee] Starting points chosen: {np.int64(82)}, L=-12.2
Starting points chosen: {np.int64(82)}, L=-12.2
[autoemcee] Starting at [0.65641769 0.20047938] +- [5.24600607e-05 7.61979786e-06]
Starting at [0.65641769 0.20047938] +- [5.24600607e-05 7.61979786e-06]
100%|██████████| 100/100 [00:00<00:00, 257.36it/s]
100%|██████████| 1000/1000 [00:03<00:00, 255.63it/s]
[autoemcee] Used 440800 calls in last MCMC run
Used 440800 calls in last MCMC run
checking convergence (iteration 2) ...
acceptance rates: [66 67 68 68 68 69 69 69]% (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
acceptance rates: [67 67 67 68 68 68 68 68]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
acceptance rates: [67 68 68 68 68 68 68 69]% (worst few)
autocorrelation length: tau=inf -> 0x lengths
[autoemcee] rhat chain diagnostic: [1.00016213 1.00020747] (<1.010 is good)
rhat chain diagnostic: [1.00016213 1.00020747] (<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.015 -0.012 +0.014 | 1 / (keV s cm2) |
| demo.spectrum.main.Sin.f | (1.003 -0.008 +0.007) x 10^-1 | rad / keV |
Values of -log(posterior) at the minimum:
| -log(posterior) | |
|---|---|
| demo | -11.256693 |
| total | -11.256693 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | 27.219268 |
| BIC | 28.504850 |
| DIC | 26.383567 |
| PDIC | 1.928229 |
[12]:
