PyIBS: Frequently Asked Questions

Contents

PyIBS: Frequently Asked Questions#

This FAQ covers PyIBS 1.5 and is curated by Luigi Acerbi. It adapts the MATLAB IBS FAQ and adds questions about the Python package.

For a tutorial with detailed examples, see the Jupyter notebook examples.

For questions not covered here, ask in the lab’s Discussions forum.

PyIBS estimates log-likelihoods from a model’s simulator. PyBADS can use these estimates to fit the parameters, and PyVBMC can use them to infer the posterior and model evidence. PyIBS, PyBADS and PyVBMC are among the lab’s tools for fitting models to data.

Acknowledgments: Many of the questions answered here originated in a live Q&A session with the Ma lab. We thank Hsin-Hung Li for taking notes.

In the answers, [1] is the IBS paper: B. van Opheusden, L. Acerbi and W. J. Ma (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.

The snippets below use np for NumPy and IBS for the estimator class:

import numpy as np
from pyibs import IBS

The snippets use sample_from_model for your simulator, R for the observed responses, S for the design, and theta for a parameter vector. Both R and S have one row per trial; they are the response_matrix and design_matrix arguments of IBS. Define these inputs for your model before running a snippet.

The fitting snippets also use a starting point x0, hard bounds lb and ub, and plausible bounds plb and pub. Each is a one-dimensional NumPy array with one element per parameter. The PyVBMC snippet additionally requires a prior. The example notebooks provide complete worked setups.

Table of contents#

General#

Which kind of problems is PyIBS suited for?#

Use PyIBS when you can simulate your model but cannot evaluate its likelihood. The problem must meet these requirements ([1], Section 2.2):

  • the responses are discrete (choices, categories, counts), or you define a matching event by binning or an approximate matching rule (see below);

  • the simulator can generate a response for any trial given its context, such as the stimulus and any relevant earlier observed stimuli and responses (conditional simulation);

  • every observed response, or the matching event you define, has a positive probability under the model. For IBS to be practical, these probabilities must also be large enough at the parameters of interest: a response of probability p takes 1 / p samples on average ([1], Section 4.2).

For each trial, IBS simulates responses until one matches the observation. For discrete responses, the number of draws gives an unbiased estimate of the trial’s log-likelihood. Summing these estimates gives an unbiased estimate of the data’s log-likelihood. IBS also estimates the variance of its estimate ([1], Section 4.3), and that variance estimate is unbiased too: its expectation, Li₂(1 − p) for a response of probability p, is the variance itself. With continuous responses, the same guarantee applies to the log-probability of the matching event, with approximation error from discretization. The answer on checking estimates explains how to assess the uncertainty and its limitations.

These estimates let you use likelihood-based methods for fitting and model comparison. PyBADS supports maximum-likelihood and maximum-a-posteriori estimation; PyVBMC approximates the posterior and model evidence. Both accept noisy targets.

Use a closed-form likelihood or an accurate analytical or numerical approximation when one is tractable. IBS can then help check its implementation and accuracy ([1], Section 6.4).

When should I use IBS rather than amortized simulation-based inference?#

Both approaches can fit models with an intractable likelihood. The better choice depends on the simulator, the trial contexts, how many datasets you will fit, and the accuracy you need.

What IBS gives. IBS estimates the log-likelihood of one dataset at one parameter vector by simulating each trial in its own context. For discrete responses, the estimator is unbiased for each dataset and parameter vector, provided the observed responses have positive probabilities, sampling runs to completion, and no call is discarded or kept for its outcome. It also estimates the variance ([1], Section 4.3), and that estimate is unbiased too. IBS requires no training or summary statistics. Its estimates can serve as a target for methods that support noisy log-likelihoods, including PyBADS for optimization and PyVBMC for posterior and evidence inference.

When amortized inference is the better choice. Amortized methods train a neural network on simulations, then reuse it across datasets. In neural posterior estimation, a trained network can produce a posterior approximation for a new dataset almost immediately (Li et al., 2026, Section 1). The training cost is shared across all the datasets it serves. This is often the better choice when you fit one model to many datasets and its simulations are cheap.

IBS, by contrast, needs new simulations for every dataset and every parameter vector evaluated during a fit, and a fit typically needs hundreds of evaluations. Libraries such as sbi and BayesFlow implement amortized inference (Li et al., 2026, Section 2.3).

What favors IBS.

  • The trials’ contexts are many and richly structured. An amortized estimator must learn the model’s behavior across the contexts it may encounter. IBS simulates only the contexts in the observed data. For example, a model of human board-game play chooses a move given the current board position, which may occur only once in the dataset ([1], Section 5.4; B. van Opheusden et al., 2023, “Expertise increases planning depth in human gameplay”, Nature 618: 1000–1005, https://doi.org/10.1038/s41586-023-06124-2). In the four-in-a-row game studied in [1], Section 5.4, the datasets used positions drawn from 5,482 positions recorded in human play. The model’s distribution over moves could not be computed even numerically, so IBS estimated its log-likelihood from simulations.

  • Guarantees on each dataset matter. “A given pre-trained amortized neural estimator may be perfectly suitable for some real datasets while it is utterly untrustworthy for others” (C. Li, A. Vehtari, P.-C. Bürkner, S. T. Radev, L. Acerbi and M. Schmitt, 2026, “Amortized Bayesian Workflow”, Transactions on Machine Learning Research, https://openreview.net/forum?id=osV7adJlKD, Section 2.2). Their workflow checks each dataset’s results with diagnostics, then uses importance sampling and MCMC when the checks fail. IBS provides unbiased log-likelihood estimates for each dataset without training, together with variance estimates. These guarantees concern the likelihood estimator. The optimizer or inference method can still make errors, so assess the fit itself, for example by comparing independent runs.

Its costs and requirements. A response of probability p takes 1 / p samples on average ([1], Section 4.2), so an improbable response can be very expensive to match. A lapse rate can put a lower bound on response probabilities, and a likelihood threshold can stop sampling at poor parameters, at the cost of bias ([1], Section 6.4 and Appendix C.1).

Responses must be discrete or have a matching event defined by binning or an approximate matching rule ([1], Section 6.3). The simulator must also generate an individual trial’s response conditioned on its observed context. A model with latent dynamics that can generate whole sequences, but cannot simulate one response conditional on an observed history, does not meet this requirement ([1], Section 2.2).

How does IBS compare with approximate Bayesian computation (ABC) or synthetic likelihood?#

ABC and synthetic likelihood commonly compare summary statistics of observed and simulated data ([1], Section 1, citing Beaumont et al., 2002, for ABC and Wood, 2010, for synthetic likelihood). Such statistics may discard information relevant to inference. IBS instead estimates the log-likelihood of the full data, trial by trial ([1], Sections 1 and 6.3). Its output can be used by likelihood-based methods that support noisy targets; ABC uses dedicated inference algorithms ([1], Section 6.3).

A synthetic likelihood gives a noisy estimate of the log-likelihood of the summary statistics. PyVBMC’s FAQ notes that a synthetic likelihood could also serve as PyVBMC’s target (Can I use any technique to estimate a noisy log-likelihood?). For continuous responses, [1] proposes approximate IBS, which borrows ABC’s tolerance but compares the full responses (see below).

What is ibs_basic for?#

ibs_basic is a minimal implementation for teaching and reading the algorithm. It follows ibs_basic.m of MATLAB IBS and the loop in [1], Section 4.1: sample one trial at a time, one response per simulator call, until it matches. The function returns one log-likelihood estimate, not its negative as IBS does by default. It has no averaging over repeats, variance estimate, sampling cap or likelihood threshold. Use IBS for analyses.

Installing PyIBS#

Where can I download PyIBS?#

Install or upgrade PyIBS with pip:

python -m pip install --upgrade pyibs

Or install with Conda:

conda install --channel=conda-forge pyibs

PyIBS 1.5 requires NumPy 2.0 or newer, and its conda-forge package requires Python 3.11 or newer. In an environment with NumPy 1.x, a package that requires it, or Python 3.10, conda may select PyIBS 0.1.0 without warning. Require version 1.5 or newer explicitly:

conda install --channel=conda-forge "pyibs>=1.5"

Conda then installs PyIBS 1.5 or reports the conflict.

See the installation instructions for more details, and the GitHub repository for the source code.

Which external packages does PyIBS require?#

PyIBS requires NumPy 2.0 or newer and SciPy 1.13 or newer. These are its only runtime dependencies, and pip or conda installs them with PyIBS. MATLAB is not required.

For fitting, install PyBADS 1.5.1 or newer, or PyVBMC 1.5 or newer. Both accept PyIBS’s estimate and estimated SD as a noisy target. Running the example notebooks also requires Jupyter; see the installation instructions.

Which version of Python do I need?#

PyIBS requires Python 3.10 or newer, and its conda-forge package Python 3.11 or newer.

How do I know whether a newer version of PyIBS exists?#

Run pyibs.check_for_updates(). It queries PyPI and prints a message. If a newer release is available, the message includes an update command: python -m pip install --upgrade pyibs for pip, or conda update --channel=conda-forge pyibs for conda. The conda-forge package can appear a few days after the PyPI release. The function also returns the installed version, the latest release and whether an update is available.

This function is PyIBS’s only network access, and it runs only when you call it. PyIBS prints no automatic messages. Use pyibs.__version__ to read the installed version without a network request. See the check_for_updates API for details.

How do I check that PyIBS works on my computer?#

The test suite ships with PyIBS. Install the test extra to add pytest, then run the tests:

python -m pip install --upgrade "pyibs[test]"
python -m pytest --pyargs pyibs

If PyBADS 1.5.1 or newer or PyVBMC 1.5 or newer is installed, the suite also fits the example model with the available package. These integration tests take a few minutes. To omit them, run python -m pytest --pyargs pyibs -m "not integration".

I am having trouble installing PyIBS. Can you help?#

Ask in the Discussions forum. Include your operating system, Python version, installation command and the full error message so that we can reproduce the problem.

The simulator and the data#

How do I write the simulator?#

Write a function sample_from_model(theta, design_rows). If it also has a parameter named rng that accepts a keyword, IBS calls it as sample_from_model(theta, design_rows, rng=rng) and supplies the object’s random number generator. Use that generator for reproducible runs.

theta is the parameter vector passed to ibs. design_rows contains the design rows for the trials requested in this simulator call. If design_matrix=None, it contains their 0-based trial indices instead. Return one simulated response per requested row:

  • For r requested rows, return an array of shape (r,) or (r, 1) for single-column responses, or (r, C) for C response columns. A response matches only when all its columns match the observation.

  • Generate each row independently. A trial may be requested more than once, so each occurrence needs a fresh draw. Sharing random draws across rows violates the sampling assumptions.

  • Use compatible response types. Numbers and booleans can match each other, but NumPy does not equate them with text or bytes.

For example, this observer makes a noisy measurement of stimulus S, reports right (1) when the measurement is at least a bias and left (-1) otherwise, and responds at random on a fraction lapse of trials:

def sample_from_model(theta, S, rng):
    log_sigma, bias, lapse = theta
    x = S + np.exp(log_sigma) * rng.standard_normal(S.shape)
    r = np.where(x >= bias, 1, -1)
    guess = rng.random(S.shape) < lapse
    r[guess] = rng.choice([-1, 1], size=np.count_nonzero(guess))
    return r


ibs = IBS(sample_from_model, R, S, random_seed=1)
neg_logl = ibs(theta)

This is the model used in the notebooks. PyIBS installs its simulator as pyibs.examples.psycho_model.psycho_generator and its closed-form negative log-likelihood as psycho_neg_logl. The IBS API describes the simulator interface and all settings.

My simulator needs additional data or inputs. How do I pass them to IBS?#

IBS passes the parameter vector, the requested design rows and, if the simulator accepts it, the random generator. Bind additional inputs to a wrapper function or use functools.partial. Unlike MATLAB ibslike, IBS has no varargin argument.

For a simulator with signature simulator(theta, design_rows, rng, data), you can write:

def sample_from_model(theta, design_rows, rng):
    return simulator(theta, design_rows, rng, data)

Define data before calling this wrapper. Alternatively, bind it with functools.partial:

from functools import partial

ibs = IBS(partial(simulator, data=data), R, S)

Both approaches retain a parameter named rng, so IBS still passes its generator to the simulator.

The responses of my model depend both on the current trial and on responses or stimuli from previous trials. How can I tell this to IBS?#

IBS estimates each response’s probability conditional on the observed history ([1], Section 2.2). Your simulator must therefore generate the requested trial’s response given the earlier stimuli and observed responses. It must not replace that history with newly simulated trials.

To give the simulator access to the full history:

  • Leave design_matrix=None: IBS then calls the simulator with the 0-based indices of the requested trials;

  • Bind the full stimulus and response arrays, and any other data, to the simulator (see the previous question);

  • For each requested index, generate a response using that trial’s stimulus and the relevant observed history.

For example, this model shifts the decision towards the previous observed response. Here S and R are one-dimensional arrays:

def make_simulator(S, R):
    def sample_from_model(theta, idx, rng):
        log_sigma, bias, pull = theta
        r_prev = np.where(idx > 0, R[idx - 1], 0)  # 0 before the first trial
        x = S[idx] + pull * r_prev
        x = x + np.exp(log_sigma) * rng.standard_normal(idx.size)
        return np.where(x >= bias, 1, -1)

    return sample_from_model


ibs = IBS(make_simulator(S, R), R, random_seed=1)  # no design: trial indices

You can also include the history in the design when each trial’s inputs fit in one row. If r_prev contains the previous observed response for each trial, use IBS(sample_from_model, R, np.column_stack([S, r_prev])) and write the simulator to read the two columns.

My model uses data structures which are not easily converted to numerical arrays. Can I still use IBS?#

Yes. The responses must support comparison, but the simulator’s other inputs can be arbitrary Python objects.

  • Encode discrete responses as numbers, booleans, text (such as "left" and "right") or bytes; IBS compares them with ==. For example, encode a board-game move as the index of its square. If a response combines numbers and text, give both the observed and simulated responses as object arrays (dtype=object).

  • Other inputs need not be numerical arrays. Bind them to the simulator and request trials by index, as explained in the question above. Alternatively, pass them as design_matrix, which accepts an array with one row per trial, including an object array of board positions or other Python objects. The simulator receives the requested entries. A list of dictionaries becomes an object array directly. For objects that NumPy may interpret as sequences, such as lists of different lengths, create an array with np.empty(N, dtype=object) and fill it element by element.

Can I use IBS with continuous responses?#

A simulated continuous response has probability zero of matching the observed value exactly, so ordinary IBS cannot estimate its density ([1], Section 6.3). You can instead define discrete matching events in either of two ways.

The simplest approach is binning. Map both observed and simulated responses to bins with the same edges:

edges = np.arange(-10, 10.5, 0.5)  # bins 0.5 wide
R = np.digitize(y, edges)  # y, the observed continuous responses


def sample_from_model(theta, S, rng):
    return np.digitize(simulate(theta, S, rng), edges)

where simulate draws continuous responses of the model.

The second approach is approximate IBS ([1], Section 6.3), inspired by approximate Bayesian computation. Define a match as a simulated response within eps of the observation. For a one-dimensional response, the approximate density is p_match / (2 * eps), where p_match is the probability of this event. IBS estimates log(p_match); subtract log(2 * eps) to estimate the log density. In more dimensions, replace 2 * eps by the volume of the matching region.

Implement this by returning a boolean for each requested trial, with every observed “response” set to True:

eps = 0.25


def within_eps(theta, idx, rng):
    return np.abs(simulate(theta, S[idx], rng) - y[idx]) <= eps


ibs = IBS(within_eps, np.ones(len(y), dtype=bool))
neg_logl, sd = ibs(theta, num_reps=10, additional_output="std")
neg_logl += len(y) * np.log(2 * eps)  # the volume term

For a fixed tolerance, the volume correction is constant across parameter vectors. It therefore changes neither the optimum nor the posterior, but it matters for model evidence or a comparison with another calculation of the log-likelihood.

Approximate IBS centers the matching region on each observation and avoids the choice of bin boundaries. Binning is often adequate if the bins are narrow compared with the noise in the data. Both approaches introduce a discretization approximation; roughly, they smooth the model’s responses on the scale of eps or half the bin width. Make this scale much smaller than the noise you wish to model. Narrower regions take more samples to match, and the expected sample count diverges as eps tends to zero ([1], Section 6.3).

What does a call of IBS return?#

ibs(theta) returns a Python float estimating the negative log-likelihood at theta. It averages num_reps independent IBS repeats, 10 by default. Set additional_output to request more information:

  • "var": the tuple (neg_logl, neg_logl_var), containing the estimate and its estimated variance;

  • "std": the tuple (neg_logl, neg_logl_std), containing the estimate and the square root of the variance estimate;

  • "full": an EstimateResult, a dictionary whose entries can also be read as attributes: neg_logl, neg_logl_var, neg_logl_std, exit_flag and its message, elapsed_time, num_samples_per_trial (the average sample count per trial), fun_count (the number of simulator calls), and the per-trial estimates neg_logl_trials and neg_logl_var_trials.

PyBADS and PyVBMC take the "std" tuple as a noisy target. return_positive=True changes the returned value to the log-likelihood. The variance estimate and the per-trial arrays keep their values and signs. trial_weights weights the trials’ log-likelihoods. A zero-weight trial is still sampled; remove a trial from the data if you want to exclude it from sampling.

With "full", exit_flag reports how sampling ended:

  • 0: every repeat completed without reaching the likelihood threshold;

  • 1: the likelihood threshold ended at least one repeat. Each ended repeat counts exactly neg_logl_threshold in the negative log-likelihood, which biases the estimate, and all per-trial estimates are NaN;

  • 2: max_time stopped sampling. The call returns a potentially biased estimate and issues a warning.

Without a likelihood threshold and with the default max_time=np.inf, every call that returns has flag 0 and an unbiased estimate. A finite max_time can bias even the calls that end with flag 0, since calls that need fewer samples, which give higher log-likelihoods, are likelier to finish in time (max_time).

PyBADS and PyVBMC#

How do I use PyIBS with PyBADS?#

PyBADS minimizes the target function, so return the negative log-likelihood estimate and its estimated SD as a tuple. Set options={"specify_target_noise": True} to tell PyBADS that the target supplies its noise level.

from pybads import BADS

ibs = IBS(sample_from_model, R, S, random_seed=1)


def target(theta):
    return ibs(theta, num_reps=100, additional_output="std")


bads = BADS(target, x0, lb, ub, plb, pub, options={"specify_target_noise": True})
optimize_result = bads.optimize()

  • Return the tuple directly from IBS. PyBADS requires a Python tuple of length two; returning a list or array raises ValueError.

  • For maximum-a-posteriori estimation, return neg_logl - log_prior(theta) with the same SD, since the log prior adds no noise.

  • PyBADS calls the target with one one-dimensional parameter vector in the original parameter space, using the same coordinates as your bounds.

  • Choose num_reps for an SD of about 1 near the optimum (see How do I choose num_reps?), and re-evaluate the solution with more repeats (below).

  • Run PyBADS from several starting points, as the PyBADS FAQ explains.

How do I use PyIBS with PyVBMC?#

PyVBMC approximates the parameter posterior and estimates the model evidence. When you give it a prior, its target should return the log-likelihood estimate (return_positive=True) and its estimated SD. PyVBMC adds the log prior itself. Set options={"specify_target_noise": True} so it uses the supplied noise level.

from pyvbmc import VBMC

ibs = IBS(sample_from_model, R, S, random_seed=1)


def log_likelihood(theta):
    return ibs(theta, num_reps=100, additional_output="std", return_positive=True)


vbmc = VBMC(log_likelihood, x0, lb, ub, plb, pub, prior=prior, options={"specify_target_noise": True})
vp, results = vbmc.optimize()

Here prior is a PyVBMC prior object, such as pyvbmc.priors.Trapezoidal(lb, plb, pub, ub). You can instead pass a function with log_prior=. Without either argument, return the log joint yourself: add log_prior(theta) to the log-likelihood estimate and keep the same SD. PyVBMC calls the target with one one-dimensional parameter vector in the original parameter space, using your bounds’ coordinates.

results["elbo"] estimates the evidence lower bound (ELBO), which PyVBMC uses as an approximation to the log model evidence. results["elbo_sd"] reports uncertainty in that estimate. Aim for target noise with an SD of about 1, “and probably not larger than ~3” in the region containing most posterior mass (the PyVBMC FAQ, How large can the noise be?).

Why does the target return the SD of the estimate, and not its variance?#

PyBADS and PyVBMC interpret the second element as a standard deviation (SD). Use additional_output="std". With "var", they would silently interpret the variance as an SD: a variance of 4 means an SD of 2, but they would read it as an SD of 4. MATLAB ibslike returns the variance unless ReturnStd is set, so check this when porting code.

Why is the SD of the estimate zero, and why do PyBADS and PyVBMC refuse it?#

Why it is zero. For a trial that first matches at sample K, IBS estimates the variance as ψ₁(1) − ψ₁(K), where ψ₁ is the trigamma function ([1], Sections 2.4 and 4.3). At K = 1, both the log-likelihood estimate and the variance estimate are zero. If every positive-weight trial matches immediately in every repeat, the whole call returns zero for both quantities.

This is a valid outcome of the estimator. When you request a variance, SD or full result, IBS returns it and issues a UserWarning linking this answer:

The IBS variance estimate is 0, as it is when every trial of positive weight matched its response at its first sample. PyBADS and PyVBMC refuse an SD of 0 for a noisy target: see https://acerbilab.org/pyibs/faq.html#faq-why-is-the-sd-of-the-estimate-zero-and-why-do-pybads-and-pyvbmc-refuse-it

It is common when the model predicts every observed response with probability close to one. For 100 trials with matching probability 0.999, the probability of a zero variance estimate is 0.37 at num_reps=10 and 0.000045 at num_reps=100.

Why PyBADS and PyVBMC refuse it. Their Gaussian-process models use the supplied SD as the noise level of the evaluation. Zero would declare the value exact, although immediate matches do not establish that it is. With positive trial weights, IBS has zero true variance only when every observed response has probability one. Otherwise a zero estimated variance is a possible sampling outcome.

Both packages require a finite, strictly positive SD for a noisy target. PyBADS raises ValueError with a message beginning The target function returned the noise SD 0.0. PyVBMC’s ValueError says the estimated SD must be a finite, positive real-valued scalar.

What to do.

  • More repeats make zero variance estimates exponentially less likely as num_reps increases. They do not rule them out: over a fit’s hundreds of evaluations, a small chance per call adds up.

  • A lapse rate in the model puts an upper bound below one on matching probabilities. If the lapse rate is at least lapse and a lapse chooses uniformly among k responses, then p ≤ 1 − lapse * (k − 1) / k. For N positive-weight trials, the probability of zero estimated variance is at most (1 - lapse * (k - 1) / k) ** (N * num_reps). This is below 1e-13 for 600 binary trials, a lapse rate of 0.01 and ten repeats. The IBS paper also recommends a positive lapse rate to control sampling cost ([1], Section 6.4).

Do I need to give PyBADS or PyVBMC the SD of every estimate?#

Supplying the SD is recommended for both packages, but neither strictly requires it:

  • For optimization with PyBADS, you can return the value alone and set options={"uncertainty_handling": True}. PyBADS then estimates the noise. This can work, since IBS’s variance stays bounded across the parameter space (see How do I choose num_reps?). Still, PyIBS computes its SD without additional simulations, and supplying it gives PyBADS more information. The PyBADS FAQ recommends this too (Should I provide an estimate of the noise associated with each evaluation?).

  • For Bayesian inference with PyVBMC, supply it. Bayesian inference is very sensitive to noisy estimates of the log-likelihood, so the inference algorithm needs all available information about the magnitude of the noise. PyVBMC can infer a noise level from a scalar-valued target with options={"uncertainty_handling": True}, but supplying the SD of each evaluation is preferable. Its FAQ recommends this (Does VBMC automatically detect that the target function is noisy?), and PyIBS provides it without additional simulations.

I have several questions about using PyBADS to optimize the log-likelihood. Can you help?#

PyBADS is a robust optimizer that works well with stochastic targets, and in particular with the noisy estimates of IBS. The PyBADS FAQ covers optimization settings, bounds, starting points and interpretation of results. Start with its section on noisy objective functions, but read the others too: all of them are relevant.

What if I want to use IBS to perform Bayesian posterior or model inference?#

We recommend Variational Bayesian Monte Carlo (PyVBMC) for approximating the posterior over model parameters and estimating the marginal likelihood, or model evidence. It supports noisy log-likelihood estimates such as IBS; see its FAQ on noisy target functions. VBMC performed well with IBS in the empirical benchmark of Acerbi, 2020.

How do I use PyIBS with PyVBMC? gives the target and settings.

Precision and IBS repeats#

A repeat is an independent run of IBS ([1], Section 4.4). A PyIBS call averages num_reps repeats, ten by default. Increasing this number makes the estimate more precise at the cost of more simulations.

How do I choose num_reps?#

For fitting, aim for an estimated SD of about 1 in the region that matters: near the optimum for PyBADS, or where most posterior mass lies for PyVBMC. The PyBADS FAQ says “a standard deviation of order 1 or less should work” (Can PyBADS handle any arbitrary amount of noise in the objective?). The PyVBMC FAQ says “ideally you want the SD of the noise in the log-likelihood to be around 1, and probably not larger than ~3” (How large can the noise be?).

For ordinary IBS, variance falls as 1 / num_reps and SD as 1 / sqrt(num_reps). The expected sample count grows in proportion to num_reps. A pilot evaluation near the expected solution, perhaps from a preliminary fit, gives a starting choice:

_, sd = ibs(theta, num_reps=10, additional_output="std")
num_reps = max(1, int(np.ceil(10 * sd**2)))  # for an SD of about 1

This uses the pilot’s variance estimate to target an SD of about 1. That estimate is noisy: use more pilot repeats, or average variance estimates from several independent calls, for a steadier guide. A pilot SD of zero does not establish that the target is deterministic; see the zero-SD answer.

One fixed num_reps usually serves a whole fit, since IBS’s variance is bounded: one trial contributes at most π²/6 ≈ 1.645 to the variance of one repeat ([1], Section 4.3), so with unit weights the SD is at most sqrt(π²/6 * N / num_reps) for N trials, wherever the fit goes. Larger noise away from the solution is usually acceptable during optimization.

If you cannot bring the SD down to about 1 within your computational budget, be as precise as you can afford.

Once the optimizer has found the best parameters, should I evaluate the log-likelihood there with more repeats?#

Yes. This is standard practice whenever a noisy target is optimized, with IBS or without: re-evaluate the candidate solution with an independent, more precise estimate. Selecting a point because its noisy value was low can make the reported minimum too optimistic ([1], Section 4.1). With PyBADS:

x_best = optimize_result["x"]
neg_logl, sd = ibs(x_best, num_reps=1000, additional_output="std")

Use an IBS object without a likelihood threshold for this evaluation. Compare models using these independent estimates and their SDs.

In an ideal world, would you let the number of repeats depend on how close the optimization algorithm thinks it is to the maximum?#

Adaptive precision can be useful and is a topic of research. PyIBS accepts num_reps separately for each call, so your target can choose a precision and return the corresponding estimated SD. PyBADS and PyVBMC call the target with the parameters alone; they do not choose its repeat count. Choose the count before drawing that evaluation’s samples. Changing the count in response to those same samples can introduce selection bias.

As I increase num_reps, the SD of the estimate goes down slowly, but the computational time increases linearly. Is this normal?#

Yes. The expected sample count grows in proportion to num_reps, while the SD falls as 1 / sqrt(num_reps). Halving the SD therefore takes four times as many repeats.

Elapsed time also depends on batching and simulator overhead. With vectorized=True, IBS requests larger batches as sampling proceeds, so the number of simulator calls can grow much more slowly than the sample count. If most time is spent on a fixed cost per call, additional repeats may therefore cost little extra time (see Should I set vectorized?).

What are trial-dependent repeats, and can I use them with PyIBS?#

Trial-dependent repeats assign a separate repeat count to each trial to minimize variance for a fixed expected sample budget ([1], Appendix C.2). A fixed allocation preserves unbiasedness as long as every trial gets at least one repeat. The optimal allocation favors probabilities near 1/2. It allocates fewer repeats to responses near probability zero, which are expensive to match, and near probability one, whose estimates are already precise.

The allocation depends on unknown trial probabilities. The paper recommends computing it from a pilot estimate with many repeats at a representative parameter vector, and retaining it for the fit.

For this to work well, the trial likelihoods should be correlated across reasonable regions of parameter space, so that an allocation computed at one parameter vector remains beneficial at the others that the fit evaluates. This seems to hold often in practice, since the improbable trials are often those where something unexpected occurred, such as a lapse; more empirical studies are needed. In the paper’s example of 500 trials with probabilities drawn uniformly between zero and one, the optimal allocation gave a median precision gain of about 1.6 over equal repeats at the same budget ([1], Appendix C.2).

PyIBS takes one num_reps for all trials in a call. To use different counts, group trials with the same count and create an IBS object for each group. Sum the independent groups’ estimates and variance estimates:

# reps: the number of repeats of each trial, from the allocation of [1]
groups = [(np.flatnonzero(reps == k), k) for k in np.unique(reps)]
ibs_groups = [(IBS(sample_from_model, R[g], S[g]), k) for g, k in groups]


def target(theta):
    neg_logl, var = 0.0, 0.0
    for ibs_g, k in ibs_groups:
        value_g, var_g = ibs_g(theta, num_reps=k, additional_output="var")
        neg_logl += value_g
        var += var_g
    return neg_logl, float(np.sqrt(var))

Without a design, each group’s simulator receives indices local to that group, starting at zero. To preserve the original trial indices, pass them as the design: IBS(sample_from_model, R[g], g).

A group whose responses nearly always match can return zero estimated variance, with the zero-SD warning. An individual group’s zero is harmless for the fitting interface when the sum of the groups’ variances is positive.

Cost, limits and the likelihood threshold#

How many samples does an estimate take?#

If trial i has matching probability p_i, it takes 1 / p_i samples per repeat on average ([1], Section 4.2). With no early stopping, a call therefore needs an expected num_reps * sum(1 / p_i) samples. Batching with vectorized=True also generates surplus samples after a trial’s last needed match, which are discarded.

Use additional_output="full" to inspect the cost: num_samples_per_trial is the total number of simulated responses divided by the number of trials, including surplus samples; fun_count counts simulator calls; and elapsed_time gives the call’s duration.

Since 1 / p = exp(-log p), the cost of matching an individual response can grow rapidly as its log-probability falls ([1], Appendix C.1). A positive lapse rate can bound this expected cost, and a likelihood threshold can stop sampling at poor parameters ([1], Section 6.4). The sampling cap max_iter raises an error when a trial takes too many samples, including when its response is impossible (see A call raises IBSSamplingError).

Should I set vectorized?#

An explicit choice can improve speed, and removes a decision based on timing from reproducibility. vectorized controls how IBS batches simulator requests:

  • True requests several samples of each trial that still needs matches. The batch size grows by the acceleration factor, 1.5 by default, subject to the batch limits. This reduces simulator calls but can generate samples beyond the last match needed.

  • False requests one sample per unfinished trial per call. It generates no surplus, but can require many calls while waiting for rare responses.

  • None, the default, times one simulation of all trials at the object’s first call with num_reps > 1. It selects False if that simulation takes at least vectorized_threshold seconds (0.1 by default), and True otherwise. The object retains the choice, readable as ibs.vectorized.

A call with num_reps=1 always requests one sample per trial per simulator call. It warns if you explicitly set vectorized=True.

To choose explicitly, consider where your simulator spends its time:

  • Expensive individual responses: False avoids surplus samples and is often preferable when runtime is proportional to sample count.

  • Expensive calls: choose True when most cost is incurred once per call, for example when starting a process, loading a model or moving data to a GPU. A slow initial simulation would make None choose False, although fewer, larger calls would be faster.

The difference can be large. For the example model with 100 to 10,000 trials and a fixed cost of 0.2 s per call, the validation of PyIBS found that 100 repeats took 15.5 to 27 minutes with False. Requesting batches with True required only 9 to 12 calls, corresponding to about 2 s of simulator overhead.

For a cheap simulator, batching can also save Python and sampler overhead. If both call overhead and per-response cost matter, time both schedules at representative parameters.

A simulator that is slow only on its first call, such as one compiled just in time, can distort the automatic decision. Warm it up before the first IBS call, or set vectorized explicitly.

What does the likelihood threshold neg_logl_threshold do, and how do I choose it?#

The threshold saves simulations at poor parameter vectors. During each repeat, IBS maintains a running lower bound on the repeat’s final negative log-likelihood estimate. Once the bound exceeds T, the repeat stops and contributes exactly T to the negative log-likelihood estimate ([1], Appendix C.1). If any repeat stops this way, the call reports exit flag 1, unless a time limit subsequently stops the call with exit flag 2.

T applies to the weighted negative log-likelihood. With return_positive=True, the same rule applies, but the stopped repeat contributes -T to the returned log-likelihood estimate.

A usual choice during optimization is the chance level: the negative log-likelihood of choosing uniformly among each trial’s possible responses. It is sum_i w_i log(k_i) for k_i responses and weight w_i on trial i, or N log(2) for N binary trials with unit weights:

ibs_fit = IBS(sample_from_model, R, S, neg_logl_threshold=len(R) * np.log(2))

Clipping changes the estimator:

  • It biases log-likelihood estimates upwards and negative log-likelihood estimates downwards. The bias becomes negligible when the true log-likelihood is sufficiently above -T ([1], Appendix C.1). The threshold is intended to identify poor parameters during optimization.

  • The variance estimate of a stopped repeat uses its counts at stopping, not the variance of the clipped value, so it is not calibrated when the threshold acts.

  • All per-trial arrays returned by additional_output="full" are NaN if any repeat was thresholded. The clipped contribution belongs to the repeat as a whole and has no per-trial decomposition.

You can use the threshold during optimization, but evaluate the final solution with an object that has no threshold: neg_logl_threshold=np.inf, the default (see above).

For Bayesian posterior or evidence inference, leave the threshold off, unless runs with and without it show that clipping leaves the posterior and evidence unchanged.

If each trial’s number of possible responses is hard to count, an average can provide a practical optimization threshold. For the four-in-a-row game, where the number of legal moves depends on the board, [1] used N log(20) (Appendix C.1).

Is it okay to stop the IBS algorithm for one trial after a fixed number of samples (e.g., 20)?#

No. Returning a value from a trial that has not matched truncates its sampling distribution and introduces bias, losing the benefit of IBS over fixed sampling ([1], Section 3). For optimization, the paper instead proposes an early-stopping threshold on the whole dataset’s estimate ([1], Appendix C.1), implemented by the likelihood threshold. The sampling cap max_iter therefore raises an error instead of returning an estimate from incomplete counts (see A call raises IBSSamplingError).

What does max_time do?#

max_time limits each call’s duration in seconds. The limit is checked after every simulator call, so a simulator call can run past it. If the limit is reached while sampling is unfinished, each trial’s estimate averages the repeats it completed. The result has exit flag 2 and a warning. A trial without a completed repeat causes IBSSamplingError. If some repeats were thresholded, they retain their clipped contributions.

A finite time limit can bias estimates because draws needing fewer samples are more likely to finish. Keep max_time=np.inf, its default, for ordinary estimation and fitting. Use the sampling cap to raise an error on excessive cost, or consider the likelihood threshold for optimization; see their respective answers for the statistical consequences. A finite time limit also makes the result depend on timing, so a seed alone no longer reproduces the run.

Troubleshooting#

A call raises IBSSamplingError. What do I do?#

IBSSamplingError usually means a trial exceeded the sampling cap. A call allows max_iter * num_reps samples for each trial, with max_iter=10**5 by default. It returns no estimate from an incomplete draw. The message gives the affected trials’ 0-based indices and sample counts:

In a draw of 10 repeats, trial 1 drew 1135 samples, more than max_iter * num_reps = 1000. 1 of 3 trials still need matches, and IBS returns no estimate for an incomplete repeat. Check that the simulator can produce every observed response at this parameter vector, or raise max_iter.

The simulator may be unable to produce an observed response, or that response may be extremely improbable at the current parameters:

  • Check the response coding and order. For example, a simulator returning 0 and 1 cannot match data encoded as -1 and 1. Return responses in the order of the requested design rows.

  • A lapse component with a small positive lower bound, such as 0.005, can make every response possible and bound the expected cost ([1], Section 6.4). It should be part of the model you intend to fit.

  • For optimization, consider a likelihood threshold to stop sampling poor parameter vectors before they reach the cap.

Raise max_iter when the response probabilities justify the sample cost and that cost is acceptable. Do not retry a failed call until one succeeds: the calls that finish under a cap are those with shorter draws, which give higher log-likelihood estimates. Investigate failures and revise the simulator or sampling settings instead.

IBSSamplingError is also raised if max_time stops a call before a trial has completed a repeat.

IBS calls my simulator with different subsets of trials. What does this mean for my simulator?#

Each requested row must have the same response distribution regardless of which other rows are requested or their order. Requests may contain a subset of trials and may repeat a trial several times.

For example, if S[:, 0] holds stimulus contrast, then np.min(design_rows[:, 0]) is the minimum over the current request. It is not the dataset’s minimum and may change on every call. If the model needs the dataset’s minimum or other global information, compute it from the full data and bind it to the simulator.

Generate each requested response independently. Drawing one random number and sharing it across all rows makes samples dependent and can invalidate the likelihood estimate or its estimated variance.

A call raises a ValueError or a TypeError about the simulator’s output. What is wrong?#

  • ValueError: The simulator was asked for ... and returned an array of shape ...: return one response per requested row. For r rows, the array must have shape (r,) or (r, 1) for one response column, or (r, C) for C columns. Check that you use the requested rows, rather than returning a response for every trial in the full dataset.

  • TypeError: The simulator returned responses of dtype ..., which NumPy never finds equal to the ... among the observed responses: the simulator returned a type that cannot match the observed responses, such as text for numerical observations. Use compatible types. For responses mixing numbers and text, use object arrays (dtype=object) for both the data and simulator output.

ValueError: response_matrix holds a NaN occurs during construction and names the affected trials. NaN is unequal to itself, so it cannot match. If it denotes a response category, encode that category with a value the simulator can return. If it denotes missing data, handle the missingness in your analysis or remove the trial.

How do I make a run reproducible?#

Seed the IBS object and draw all simulator randomness from its supplied generator:

def sample_from_model(theta, S, rng):
    ...  # every random draw from rng


ibs = IBS(sample_from_model, R, S, vectorized=True, random_seed=1)

random_seed follows PyBADS’s convention:

  • An integer or numpy.random.SeedSequence seeds a new generator.

  • A numpy.random.Generator is used directly.

  • None, the default, derives a new generator from NumPy’s global random state. Calling np.random.seed before constructing the object therefore also fixes its seed.

The object stores its generator as ibs.rng and passes it to simulators with a keyword-compatible parameter named rng.

Two objects with the same seed reproduce the same sequence of estimates when given the same calls and simulator, provided timing does not change the sampling:

  • Set vectorized=True or False. Automatic None also reproduces draws if it makes the same decision in both runs, but that decision depends on timing (see Should I set vectorized?).

  • Keep acceleration_threshold=None and max_time=np.inf, their defaults. Neither then uses simulator runtime to alter sampling.

A simulator using global functions such as np.random.normal is not controlled by the seed passed to IBS. Reproducibility across computers or different Python, NumPy and SciPy versions is also not guaranteed.

For fitting, give separate seeds to PyIBS and the optimizer or inference method:

ibs = IBS(sample_from_model, R, S, vectorized=True, random_seed=1)
bads = BADS(target, x0, lb, ub, plb, pub, options={"specify_target_noise": True, "random_seed": 2})
vbmc = VBMC(log_likelihood, x0, lb, ub, plb, pub, prior=prior, options={"specify_target_noise": True}, seed=2)

Create and seed the IBS object once for the whole fit. Recreating it with the same seed at every evaluation reuses the random stream and freezes the noise, which can bias fitting (the PyBADS FAQ, Can I make a noisy objective function deterministic by fixing the noise process?).

How do I check that the estimates are right for my model?#

Check the simulator against a likelihood computed independently wherever possible ([1], Section 6.4). This may be a closed-form or accurate numerical likelihood for the full model, a simplified model, or a few carefully chosen trials and parameter vectors. PyIBS itself has been validated across models and settings, but that does not check your simulator.

Draw many independent IBS estimates at a fixed parameter vector, without a likelihood threshold or time limit. Compare their mean with the reference value, allowing for the mean’s sampling error. Also compare their empirical variance with the mean of the reported variance estimates. Ordinary IBS gives unbiased estimates of both log-likelihood and variance when its draws are completed and not selected.

If the estimates are approximately normal and the reported SDs are stable and positive, inspect the z-scores (neg_logl - exact) / sd. They should be approximately centered on zero with SD near one. This rests on a normal approximation, which unbiasedness does not guarantee. If nearly all trials match on their first samples, estimated variances can be zero and z-scores can behave poorly even with many repeats; see the zero-SD answer.

The installed example model, pyibs.examples.psycho_model, provides a simulator and the closed-form negative log-likelihood psycho_neg_logl. The first example notebook shows this comparison on 600 trials, where the standardized errors are close to normal.

Check the fit separately, for example by comparing independent runs of PyBADS or PyVBMC from different starting points.

Miscellanea#

I used ibslike in MATLAB. What is different in PyIBS?#

PyIBS implements the same IBS estimator as ibslike.m 0.96 of MATLAB IBS, with a Python interface. In PyIBS:

  • Create an IBS object with the simulator, responses and design: ibs = IBS(sample_from_model, response_matrix, design_matrix, ...). Evaluate it with ibs(params, num_reps, trial_weights, additional_output, return_positive). The corresponding MATLAB call is ibslike(fun, params, respMat, designMat, options, varargin). Extra simulator arguments are bound to it rather than passed as varargin.

  • MATLAB’s Nreps, TrialWeights and ReturnPositive become the call arguments num_reps, trial_weights and return_positive. ReturnStd becomes additional_output="std". additional_output="full" combines the exit flag and output structure in one result. Its fields use Python names: funcCount becomes fun_count, NsamplesPerTrial becomes num_samples_per_trial, and nlogL_trials and nlogLvar_trials become neg_logl_trials and neg_logl_var_trials.

  • Other options become object settings: Vectorized becomes vectorized ('auto' becomes None), Acceleration becomes acceleration, NsamplesPerCall becomes num_samples_per_call, MaxIter becomes max_iter, MaxTime becomes max_time, and NegLogLikeThreshold becomes neg_logl_threshold (np.inf turns it off). MATLAB’s hard-coded MaxSamples, AccelerationThreshold, VectorizedThreshold and MaxMem are also settings: max_samples, acceleration_threshold, vectorized_threshold and max_mem.

  • With no design, the simulator receives 0-based trial indices.

  • The simulator can take the object’s random generator, rng, for reproducible runs.

The sampling and stopping rules also have deliberate differences:

  • vectorized=None makes its decision once per object.

  • The number of samples per call grows after every call by default. Set acceleration_threshold=0.1 to use ibslike’s timing rule.

  • The likelihood threshold checks each repeat separately and assigns a stopped repeat exactly the threshold value ([1], Appendix C.1).

  • The sample cap counts samples separately for each trial and raises IBSSamplingError when exceeded.

  • A time-limited estimate averages each trial’s completed repeats and issues a warning.

For the tests corresponding to ibslike('test'), run pytest --pyargs pyibs. The catalogue of deliberate differences lists all differences from ibslike.m and the didactic ibs_basic.m, with their reasons. The two implementations draw random numbers differently, so runs do not match draw for draw even with the same simulator.

I used PyIBS 0.1.0. What do I need to change?#

PyIBS 1.5 keeps the main calling convention: construct IBS(...) and evaluate it with ibs(params, num_reps, ...). Check an existing script for the following changes:

  • Results differ because of the new sampler and fixes to 0.1.0, including its default max_iter=15, which biased estimates for improbable responses. Calling np.random.seed does not reproduce 0.1.0’s results.

  • Settings, data and call arguments are checked more strictly. Values accepted by 0.1.0 may raise an exception.

  • Exceeding the sample cap raises IBSSamplingError instead of returning exit flag 3.

  • The example model is imported from pyibs.examples.psycho_model.

The “Upgrading from 0.1.0” list in the changelog gives the full set of checks, including changes to returned values and object settings.

Contents