IBS

Note

The IBS class uses inverse binomial sampling to estimate a model’s log-likelihood from simulated discrete responses.

Create an IBS object with the model’s simulator, the observed responses and, optionally, the design of each trial; then call it with a parameter vector. Each call draws a new estimate of the negative log-likelihood (the log-likelihood with return_positive=True). Request its estimated variance with additional_output="var", its estimated standard deviation with "std", or an EstimateResult with "full". Exceeding the sample cap raises IBSSamplingError. A time limit can return a partial estimate with a warning; it raises the error if any trial has no completed repeat.

With additional_output="std", a call returns the pair that PyBADS and PyVBMC take from a noisy target: the quick start shows both. The catalogue of deliberate differences lists where IBS differs from ibslike.m of MATLAB IBS, the reference implementation, and why.

class pyibs.IBS(sample_from_model, response_matrix, design_matrix=None, vectorized=None, acceleration=1.5, num_samples_per_call=0, max_iter=100000, max_time=inf, max_samples=10000, acceleration_threshold=None, vectorized_threshold=0.1, max_mem=None, neg_logl_threshold=inf, *, random_seed=None)

Estimate a model’s log-likelihood by inverse binomial sampling.

An IBS object holds a simulator, the observed responses and the sampling settings. Call it with a parameter vector to estimate the negative log-likelihood of the data and, optionally, its variance. IBS simulates each trial until a response matches the observation ([1]). The estimator is unbiased when sampling is allowed to complete; a likelihood threshold or time limit can introduce bias.

Parameters:

sample_from_modelcallable

Called as sample_from_model(params, design_rows). If it has a parameter named rng that accepts a keyword, IBS also passes its random generator as rng=rng.

design_rows contains the requested rows of design_matrix, or their 0-based trial indices when no design is supplied. A trial can appear more than once in a call. Return an independent draw for every requested row: an array of shape (r,) or (r, 1) for single-column responses, or (r, C) for C > 1 response columns, where r is the number of requested rows.

response_matrixarray_like of shape (N,) or (N, C)

Observed responses, one row per trial; a scalar represents one trial. A simulated response matches only if every column agrees. Responses must be discrete and comparable with ==: numbers, booleans, text, bytes, or objects. The simulator must return the same kind of response: text and bytes do not match numbers. For responses that mix numbers and text, give both the observed and the simulated responses as object arrays (dtype=object): otherwise NumPy converts the numbers to text, and responses can fail to match.

design_matrixarray_like of shape (N, …), optional

Experimental conditions or other simulator inputs, one row per trial. With None (the default), the simulator receives trial indices instead.

vectorizedbool or None, optional

The sampling schedule. True requests several samples per open trial in each simulator call; acceleration controls how this number grows. False requests one sample per open trial per call. None (the default) chooses a schedule on the first call with num_reps > 1 by timing one simulation of all trials. It chooses False if that simulation takes at least vectorized_threshold seconds, and True otherwise. The choice is retained for later calls and exposed through the vectorized attribute.

This timing cannot distinguish a fixed cost per simulator call from a cost per response. It can choose False for a simulator with a large fixed cost per call, which True spreads over more samples; set True for such a simulator. If its first call compiles code or performs other setup, warm it up before IBS times it. With num_reps=1, IBS uses the False schedule and warns if True was explicitly requested.

accelerationfloat, optional

Factor by which the requested samples per trial grow between simulator calls. Must be finite and >= 1. Default 1.5.

num_samples_per_callint, optional

Initial number of samples per trial per simulator call, subject to max_samples and max_mem. The default, 0, starts at num_reps samples. Unused with vectorized=False.

max_iterint, optional

Sample cap per trial, expressed per repeat. A call requesting num_reps repeats raises IBSSamplingError if any trial draws more than max_iter * num_reps samples. Default 10**5. A very large value, such as 10**18, effectively removes the cap; sampling can then continue indefinitely if an observed response has zero probability under the model.

max_timefloat, optional

Time limit of each call in seconds, > 0. Checked after each simulator call. The default, np.inf, imposes no limit. If the limit stops sampling, IBS averages each trial’s completed repeats, returns exit flag 2 and warns. A trial with no completed repeat raises IBSSamplingError. Repeats ended by the likelihood threshold contribute neg_logl_threshold to the negative log-likelihood before averaging.

A finite limit can bias the returned estimates, since repeats that need fewer samples, which give higher log-likelihoods, are likelier to finish in time.

max_samplesint, optional

Maximum samples per trial in one simulator call. Default 10**4.

acceleration_thresholdfloat or None, optional

None (the default) grows the requested samples after every simulator call. A positive time in seconds grows them only after calls faster than that time. With a finite threshold, wall-clock timing affects the draws, so a seed alone cannot reproduce a run.

vectorized_thresholdfloat, optional

Time in seconds at or above which vectorized=None selects False. Must be > 0. Default 0.1.

max_memint or None, optional

Approximate limit on samples in one simulator call. Each open trial receives at most ceil(max_mem / n_open) samples, where n_open is the number of trials still being sampled. Rounding can exceed the limit by fewer than n_open samples. None (the default) sets max(min(N, 10**4), 10) * 100.

neg_logl_thresholdfloat, optional

Threshold T > 0 for the weighted negative log-likelihood. A repeat stops once its accumulating estimate exceeds T and contributes T to the negative log-likelihood (or -T to the log-likelihood). This saves simulations at poor parameter vectors but biases the log-likelihood estimate upwards ([1], Appendix C.1).

For optimization, a common choice is the chance-level negative log-likelihood: sum_i w_i log(k_i) for k_i equally likely responses and trial weight w_i, or N log 2 for N binary trials with unit weights. For Bayesian inference, leave the threshold disabled unless its effect on the posterior and model evidence has been assessed. The default, np.inf, disables it.

random_seedNone, int, numpy.random.SeedSequence or numpy.random.Generator, optional

Seed for rng; keyword only. None (the default) derives a new generator from NumPy’s global random state, so calling np.random.seed before construction fixes its initial state. An integer or SeedSequence seeds a new generator. An existing Generator is used directly.

Attributes:

rngnumpy.random.Generator

Random generator passed to sample_from_model when it has a parameter named rng that accepts a keyword.

vectorizedbool or None

Selected sampling schedule.

Raises:

TypeError

If sample_from_model is not callable or a setting has the wrong type, such as a boolean for a count or a time.

ValueError

If response_matrix is empty, has more than two dimensions, or contains a NaN, an element not equal to itself, which no simulated response matches; if design_matrix does not have N rows; or if a setting is out of range.

Notes

Settings. Constructor parameters other than random_seed are available as read-only attributes. The response and design arrays are read-only copies. max_mem contains the resolved sample limit; vectorized contains the selected schedule. Counts take integers or whole-number floats, such as 1e5, and are stored as integers; other numeric settings are stored as floats.

Reproducibility. The simulator must draw from the rng it receives. Two objects created with the same seed then reproduce the same sequence of estimates when timing does not affect the sampling: set vectorized explicitly to True or False, and keep acceleration_threshold and max_time at their defaults. Automatic schedule selection also reproduces the draws if it makes the same choice in both runs.

Cost. If the observed response has probability p, its trial takes 1 / p samples per repeat on average, before any stopping rule. The sample cap stops runs with exceptionally rare or impossible responses. A likelihood threshold can reduce the cost at poor parameter vectors.

MATLAB reference. pyibs/README.md catalogues deliberate differences from MATLAB ibslike.m and explains their reasons.

References

[1] (1,2)

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

__call__(params, num_reps=10, trial_weights=None, additional_output=None, return_positive=False)

Estimate the negative log-likelihood of params.

Parameters:

paramsarray_like

The parameter vector, passed to the simulator as given.

num_repsint, optional

Number of independent IBS repeats to average, at least 1. Whole-number floats are accepted as integers. Default 10.

trial_weightsNone, float or array_like of shape (N,), optional

Non-negative, finite real weights of the trials’ log-likelihoods. Booleans and strings are rejected. A scalar applies the same weight to every trial; None (the default) uses unit weights. Zero-weight trials are still sampled and can reach the sample cap or time limit. Remove trials from the data to exclude them entirely.

additional_outputNone or str, optional

None or "none" returns only the estimate. "var" adds its estimated variance; "std" adds the square root of that variance. "full" returns an EstimateResult with diagnostics and per-trial estimates.

return_positivebool, optional

Whether to return the log-likelihood rather than the negative log-likelihood; the variance estimate and the per-trial arrays of "full" are unchanged. Default False.

Returns:

neg_loglfloat

The negative log-likelihood estimate (the log-likelihood with return_positive=True), when additional_output is None.

(neg_logl, neg_logl_var)tuple of two floats

Estimate and estimated variance, with "var".

(neg_logl, neg_logl_std)tuple of two floats

Estimate and estimated standard deviation, with "std". PyBADS and PyVBMC accept this pair from a noisy target.

resultEstimateResult

Estimate and diagnostics, with "full".

Raises:

IBSSamplingError

If a trial draws more than max_iter * num_reps samples, or max_time stops the sampling before a trial has completed a repeat.

ValueError

If num_reps, trial_weights or additional_output is invalid, or the simulator returns an array with the wrong shape.

TypeError

If num_reps, trial_weights or return_positive has the wrong type, or the simulator returns responses of a kind that cannot match the observations, such as text for numbers.

Warns:

UserWarning

When max_time stops the sampling (exit flag 2); when a call returns a variance estimate of 0, which gives an SD that PyBADS and PyVBMC refuse; or when vectorized=True is used with num_reps=1.

property acceleration

Growth factor for samples requested per trial per simulator call.

property acceleration_threshold

Simulator-call time in seconds below which sample requests grow.

None grows the requests after every call.

property design_matrix

Simulator inputs per trial as a read-only array, or None.

property max_iter

The cap on the samples of one trial per repeat.

property max_mem

Resolved sample limit per simulator call, subject to rounding.

property max_samples

Maximum samples per trial in one simulator call.

property max_time

The time limit of a call, in seconds.

property neg_logl_threshold

Weighted negative log-likelihood threshold; inf disables it.

property num_samples_per_call

Initial samples per trial per simulator call; 0 uses num_reps.

property response_matrix

The observed responses, a read-only array.

property sample_from_model

The simulator.

property vectorized

Selected sampling schedule. With automatic selection, remains None until the first call with num_reps > 1.

property vectorized_threshold

Time in seconds at or above which automatic scheduling selects vectorized=False.