PyIBS Example 2: Maximum-likelihood estimation with PyBADS#
Fit a model using only its simulator: this notebook passes IBS estimates of a negative log-likelihood to PyBADS, then checks the fitted parameters using a closed-form likelihood.
We use the orientation discrimination model and data set introduced in the first notebook. This is the second of three examples; the third performs Bayesian inference with PyVBMC. The code is also available as a Python script.
See the installation guide for PyIBS and the packages needed to run this notebook.
import sys
if "google.colab" in sys.modules: # Colab lacks PyIBS and PyBADS: install them
%pip install "pyibs>=1.5" "pybads>=1.5.1"
import math
import numpy as np
from pybads import BADS
from scipy.optimize import minimize
from pyibs import IBS
from pyibs.examples.psycho_model import psycho_generator, psycho_neg_logl
0. Choosing an optimizer for a noisy likelihood#
Repeated IBS evaluations at the same parameters give different values. An optimizer needs to distinguish that sampling noise from a useful improvement. Bayesian Adaptive Direct Search (BADS) supports noisy objectives; its Python implementation, PyBADS, can use the estimated SD of each evaluation when the target supplies it.
1. The model and its search region#
The parameter vector is \(\theta=(\eta,\mu,\gamma)\): log sensory-noise SD, response bias and lapse rate. We generate the same 600 trials as in the first notebook: normally distributed orientations of SD 3°, and responses at \(\theta_\text{true}=(\log 1,0.2,0.03)\).
Seeds spawned from a common root give separate streams for data generation, IBS, the starting point and PyBADS. The data stream is shared across this series; the other streams keep sampling and optimization reproducible.
# Independent seeds for the data, the IBS estimates, the start and PyBADS
seed_data, seed_ibs, seed_start, seed_bads = np.random.SeedSequence(
2026
).spawn(4)
theta_true = np.array([math.log(1.0), 0.2, 0.03]) # log(sigma), bias, lapse
n_trials = 600
rng = np.random.default_rng(seed_data)
S = 3 * rng.standard_normal((n_trials, 1)) # Orientation of each trial
R = psycho_generator(theta_true, S, rng) # Response of each trial
The bounds have two roles:
Hard bounds
LBandUBdefine where PyBADS may evaluate the model. The lapse rate is at least 0.01, so either response has probability at least 0.005 on every trial. This limits the cost of IBS: a near-impossible response would otherwise require many samples before a match.Plausible bounds
PLBandPUBdescribe a likely region for the solution. PyBADS uses that box to place initial points and scale the search.
These bounds and the data-generation recipe follow MATLAB IBS’s ibs_example.m.
LB = np.array([math.log(0.1), -2.0, 0.01]) # Lower bounds
UB = np.array([math.log(10.0), 2.0, 1.0]) # Upper bounds
PLB = np.array([math.log(0.2), -1.0, 0.02]) # Plausible lower bounds
PUB = np.array([math.log(5.0), 1.0, 0.2]) # Plausible upper bounds
For a reference fit, we minimize the closed-form negative log-likelihood with L-BFGS-B. This target is deterministic and smooth. Its numerical minimum gives us a benchmark for the later fit with IBS.
exact_fit = minimize(
psycho_neg_logl,
(PLB + PUB) / 2,
args=(S, R),
method="L-BFGS-B",
bounds=list(zip(LB, UB)),
)
theta_ml, neg_logl_min = exact_fit.x, exact_fit.fun
print(f"Reference maximum-likelihood point: {theta_ml.round(4)}")
print(f"Reference minimum of the negative log-likelihood: {neg_logl_min:.2f}")
Reference maximum-likelihood point: [0.0584 0.2828 0.0184]
Reference minimum of the negative log-likelihood: 151.78
2. Defining the IBS target#
PyBADS minimizes a target value. We return the negative log-likelihood estimate and its estimated SD as the tuple (value, sd), using additional_output="std" and options={"specify_target_noise": True}.
The second component is measured in the same units as the target value. For example, variance 9 corresponds to SD 3; passing the variance would give PyBADS the wrong noise scale.
We use 100 repeats, which give an SD near 1 for this data set. PyBADS works best with noise of about 1 or less near the solution, though larger noise farther away is usually acceptable (see the PyBADS FAQ).
PyBADS requires a finite, positive SD. With the unit weights used here, IBS reports zero when every trial matches on its first sample in every repeat. That is extremely unlikely for this data set, but IBS warns if it happens. The PyIBS FAQ explains how to handle it.
ibs = IBS(psycho_generator, R, S, vectorized=True, random_seed=seed_ibs)
def target(theta):
"""IBS estimate of the negative log-likelihood, and its SD."""
return ibs(theta, num_reps=100, additional_output="std")
We draw a starting point uniformly within the plausible box, as ibs_example.m does, and compare one IBS evaluation there with the closed-form value. This gives an initial check that the target has the intended sign and noise scale.
x0 = PLB + np.random.default_rng(seed_start).random(3) * (PUB - PLB)
value, sd = target(x0)
print(f"Starting point: {x0.round(4)}")
print(f"Target at the start: {value:.2f} +/- {sd:.2f}")
print(f"Exact value: {psycho_neg_logl(x0, S, R):.2f}")
Starting point: [ 0.1637 -0.5321 0.1351]
Target at the start: 192.54 +/- 1.24
Exact value: 194.46
3. Running the optimization#
The BADS object receives the target, starting point and bounds. specify_target_noise tells it to read the target’s (value, sd) output; random_seed fixes PyBADS’s own random stream. The separate seed of the IBS object controls the simulator draws.
optimize() prints progress as it searches. It may also print a documentation tip; setting options["show_tips"] = False disables tips.
bads = BADS(
target,
x0,
LB,
UB,
PLB,
PUB,
options={"specify_target_noise": True, "random_seed": seed_bads},
)
optimize_result = bads.optimize()
Variables (index) internally transformed to log coordinates: [2].
Beginning optimization of a STOCHASTIC objective function (specified noise)
Tip: For angles and other periodic parameters, list their indices (from 0) in options['periodic_vars'], with the period as their hard bounds (e.g. -np.pi and np.pi for an angle in radians), usually as their plausible bounds too. PyBADS then treats the two bounds as the same point and can move across them.
https://acerbilab.github.io/pybads/faq.html#faq-does-pybads-support-periodic-variables-such-as-angles
Iteration f-count E[f(x)] SD[f(x)] MeshScale Method Actions
0 1 194.397 1.24537 1
0 33 156.692 1.24537 1 Initial mesh Initial points
1 38 156.692 1.10304 0.5 Refine grid Train
2 39 155.125 1.05164 0.5 Successful search (ES-wcm)
2 40 150.852 0.708976 0.5 Successful search (ES-wcm)
2 52 150.852 0.708976 0.25 Refine grid
3 53 151.446 0.324128 0.25 Successful search (ES-ell)
3 66 151.446 0.324128 0.125 Refine grid Train
4 67 151.805 0.233659 0.125 Incremental search (ES-wcm)
4 76 151.805 0.233659 0.0625 Refine grid Train
5 77 151.829 0.188398 0.0625 Successful search (ES-wcm)
5 78 151.813 0.186596 0.0625 Successful search (ES-wcm)
5 90 151.813 0.186596 0.03125 Refine grid Train
6 91 151.949 0.156924 0.03125 Successful search (ES-ell)
6 104 151.949 0.156924 0.015625 Refine grid Train
7 114 151.973 0.1346 0.03125 Successful poll Train
8 118 152.028 0.12801 0.015625 Refine grid
9 119 152.049 0.124006 0.015625 Successful search (ES-wcm)
9 122 152.022 0.125891 0.015625 Successful search (ES-wcm)
9 132 152.022 0.125891 0.0078125 Refine grid
10 134 152.021 0.116963 0.0078125 Successful search (ES-wcm)
10 135 152.016 0.123814 0.0078125 Successful search (ES-wcm)
10 139 152.007 0.125444 0.0078125 Successful search (ES-wcm)
10 141 152.007 0.130882 0.0078125 Incremental search (ES-ell)
10 144 152 0.13896 0.0078125 Successful search (ES-wcm)
10 146 151.992 0.145943 0.0078125 Successful search (ES-wcm)
10 158 151.992 0.145943 0.00390625 Refine grid Train
11 161 151.96 0.0990209 0.00390625 Successful search (ES-wcm)
11 172 151.956 0.0949202 0.0078125 Successful poll
12 173 151.95 0.095887 0.0078125 Successful search (ES-wcm)
12 177 151.941 0.0986764 0.0078125 Successful search (ES-ell)
12 190 151.938 0.093241 0.015625 Successful poll Train
13 193 151.941 0.0888689 0.015625 Successful search (ES-ell)
13 204 151.941 0.0888689 0.0078125 Refine grid
14 205 151.95 0.0858407 0.0078125 Successful search (ES-wcm)
14 218 151.95 0.0858407 0.00390625 Refine grid
15 219 151.946 0.085222 0.00390625 Successful search (ES-ell)
15 220 151.918 0.0911543 0.00390625 Successful search (ES-ell)
15 222 151.913 0.0963444 0.00390625 Successful search (ES-ell)
15 224 151.904 0.10288 0.00390625 Successful search (ES-ell)
15 225 151.897 0.110846 0.00390625 Successful search (ES-ell)
15 229 151.892 0.111758 0.00390625 Successful search (ES-ell)
15 230 151.888 0.116576 0.00390625 Successful search (ES-ell)
15 240 151.888 0.116576 0.00195312 Refine grid Train
16 243 151.849 0.124858 0.00195312 Successful search (ES-wcm)
16 244 151.835 0.127956 0.00195312 Successful search (ES-ell)
16 246 151.827 0.133615 0.00195312 Successful search (ES-ell)
16 247 151.803 0.138115 0.00195312 Successful search (ES-ell)
16 250 151.764 0.15164 0.00195312 Successful search (ES-wcm)
16 262 151.764 0.15164 0.000976562 Refine grid
17 266 151.787 0.107299 0.000976562 Successful search (ES-wcm)
17 267 151.771 0.108201 0.000976562 Successful search (ES-ell)
17 268 151.757 0.107718 0.000976562 Successful search (ES-wcm)
17 269 151.664 0.148904 0.000976562 Successful search (ES-wcm)
17 271 151.648 0.150077 0.000976562 Successful search (ES-ell)
17 284 151.648 0.150077 0.000488281 Refine grid Train
18 294 151.737 0.110366 0.00012207 Refine grid
19 299 151.703 0.100431 0.000244141 Successful poll Train
20 300 151.707 0.100021 0.000244141 Successful search (ES-wcm)
20 301 151.696 0.100531 0.000244141 Successful search (ES-ell)
20 302 151.694 0.101579 0.000244141 Successful search (ES-ell)
20 303 151.683 0.101221 0.000244141 Successful search (ES-wcm)
20 304 151.676 0.101861 0.000244141 Successful search (ES-ell)
20 306 151.675 0.101999 0.000244141 Incremental search (ES-ell)
20 317 151.675 0.101999 0.00012207 Refine grid Train
21 327 151.69 0.143423 6.10352e-05 Refine grid
22 337 151.732 0.140097 3.05176e-05 Refine grid
23 347 151.72 0.133224 1.52588e-05 Refine grid Train
24 357 151.732 0.095507 7.62939e-06 Refine grid
Optimization terminated: change in the function value less than options['tol_fun'].
Estimated function value at minimum: 151.47354237613908 ± 0.3335538477172075 (precision-weighted mean ± SEM from 10 samples)
4. Assessing the solution#
optimize_result["x"] contains the fitted parameters. fval is the precision-weighted mean of PyBADS’s final target evaluations there, and fsd is the standard error of that mean. We compare the parameters with the generating vector and the reference fit.
theta_bads = optimize_result["x"]
fval, fsd = optimize_result["fval"], optimize_result["fsd"]
print(f"PyBADS estimate at its solution: {fval:.2f} +/- {fsd:.2f}")
print(f"Target evaluations: {optimize_result['func_count']}")
print()
print(f"{'':17}{'eta':>8}{'bias':>8}{'lapse':>8}")
for name, theta in [
("Generating", theta_true),
("Reference ML", theta_ml),
("PyBADS solution", theta_bads),
]:
print(f"{name:17}" + "".join(f"{x:8.4f}" for x in theta))
PyBADS estimate at its solution: 151.47 +/- 0.33
Target evaluations: 368
eta bias lapse
Generating 0.0000 0.2000 0.0300
Reference ML 0.0584 0.2828 0.0184
PyBADS solution 0.0823 0.2913 0.0209
A fresh, more precise evaluation helps assess a solution chosen using noisy values. We use 1,000 IBS repeats, ten times the number used during the search, as ibs_example.m does, and compare that estimate with the closed-form likelihood at the fitted parameters. We also measure how far the fitted value lies above the reference minimum.
neg_logl, neg_logl_sd = ibs(theta_bads, num_reps=1000, additional_output="std")
neg_logl_exact = psycho_neg_logl(theta_bads, S, R)
print(
"IBS estimate at the solution (1,000 repeats): "
f"{neg_logl:.2f} +/- {neg_logl_sd:.2f}"
)
print(f"Exact value at the solution: {neg_logl_exact:.2f}")
print(f"Reference minimum: {neg_logl_min:.2f}")
print(f"Difference: {neg_logl_exact - neg_logl_min:.2f}")
IBS estimate at the solution (1,000 repeats): 152.05 +/- 0.33
Exact value at the solution: 151.85
Reference minimum: 151.78
Difference: 0.08
The 1,000-repeat estimate lies within two of its SDs of the exact value at the solution, and the fitted negative log-likelihood is close to the reference minimum relative to the roughly unit noise of each search evaluation. This run illustrates that a noisy target can still locate a useful solution; the closed form lets us check the fit independently.
The reference maximum-likelihood point differs from the generating parameters because a finite data set has sampling variation, which shrinks with more trials. The PyBADS solution also differs from the reference point, since PyBADS searched a noisy target.
For a model whose likelihood cannot be checked directly, run PyBADS from several starting points and compare the solutions with fresh, more precise evaluations. The PyBADS FAQ describes this workflow.
5. Moving to Bayesian inference#
The IBS target supplies both a negative log-likelihood estimate and its uncertainty. PyBADS uses those evaluations to search for maximum-likelihood parameters. The next notebook uses the same simulator to infer a distribution over parameters and estimate the model evidence with PyVBMC.
References#
van Opheusden, B., Acerbi, L. & Ma, W. J. (2020). Unbiased and efficient log-likelihood estimation with inverse binomial sampling. PLOS Computational Biology 16(12): e1008483. https://doi.org/10.1371/journal.pcbi.1008483
Singh, G. S. & Acerbi, L. (2024). PyBADS: Fast and robust black-box optimization in Python. Journal of Open Source Software 9(94): 5694. https://doi.org/10.21105/joss.05694
Acerbi, L. & Ma, W. J. (2017). Practical Bayesian optimization for model fitting with Bayesian adaptive direct search. In Advances in Neural Information Processing Systems 30: 1834-1844. https://arxiv.org/abs/1705.04405