Fitting and comparison

Fitting, model comparison and model-recovery studies, common to every family.

• fit() fits one model class by multi-start optimisation (maximum likelihood or least squares).

• compare() fits several classes and ranks them by BIC or AIC.

• recovery() checks that the comparison selects the model that generated simulated data.

• fit_individuals() and compare_individuals() fit every person separately, which tests whether one set of parameters describes everyone.

Examples

>>> from quantum_mind import fit, compare
>>> from quantum_mind.families.order_effects import QuantumOrderModel4D, BayesOrderModel
>>> data = QuantumOrderModel4D().sample(None, 500)
>>> ranking = compare([QuantumOrderModel4D, BayesOrderModel], data)
>>> ranking[0].summary()['model']
'QuantumOrderModel4D'

Classes

FitResult	Result of fit().
IndividualFits	Per-person fits of one model class, returned by fit_individuals().

Functions

fit	Fit a model class to data by multi-start optimisation.
compare	Fit several model classes to the same data and rank them.
recovery	Model-recovery study.
kl	Kullback-Leibler divergence.
tvd	Total variation distance.
fit_individuals	Fit a model class separately to each person's data.
compare_individuals	Fit every candidate to every person and compare them.

class FitResult(model, loss, loglik, k, n, options=<factory>)

Bases: object

Result of fit().

Variables:

• model (Model) – Fitted model instance.

• loss (float) – Negative log-likelihood (multinomial models) or sum of squared errors ('sse' models).

• loglik (float) – Log-likelihood. For 'sse' models this is the Gaussian log-likelihood with the error variance profiled out, so that AIC and BIC are comparable across models.

• k (int) – Number of fitted parameters.

• n (int) – Number of observations.

• options (dict) – Structural options of the selected fit.

Parameters:

• model (object)

• loss (float)

• loglik (float)

• k (int)

• n (int)

• options (dict)

property aic

Akaike information criterion \(2k - 2\log L\).

Type:

float

property bic

Bayesian information criterion \(k\log n - 2\log L\).

Type:

float

summary()

Summary of the fit.

Returns:

dict – Keys model, params, options, loglik, k, n, aic and bic.

class IndividualFits(model_name, results)

Bases: object

Per-person fits of one model class, returned by fit_individuals().

Every person has their own parameters, so group criteria are sums over people.

Variables:

• model_name (str) – Name of the fitted class.

• results (dict) – {person: FitResult}.

Parameters:

• model_name (str)

• results (dict)

property loglik

sum of the per-person log-likelihoods.

Type:

float

property k

total number of fitted parameters (per person times number of people).

Type:

int

property n

total number of observations.

Type:

int

property bic

group BIC, the sum of the per-person BICs.

Type:

float

property aic

group AIC, the sum of the per-person AICs.

Type:

float

parameters()

Fitted parameter values across people.

Returns:

dict – {parameter name: numpy.ndarray over people}.

fit(model_cls, data, design=None, restarts=8, rng=None, structures=None, method='Nelder-Mead', x0=None, **options)

Fit a model class to data by multi-start optimisation.

Parameters:

• model_cls (type) – Subclass of Model.

• data (dict) – {condition: counts} (multinomial models) or {condition: values} ('sse' models).

• design (iterable, optional) – Conditions passed to predict; most families accept None for their default design.

• restarts (int, optional) – Number of random starting points per structure.

• rng (numpy.random.Generator, optional) – Random number generator for the starting points (default: seed 0, so fits are reproducible).

• structures (list of dict, optional) – Discrete structural choices (for example projector ranks). Every structure is fitted and the best one kept.

• method (str, optional) – scipy.optimize.minimize() method.

• x0 (array_like, optional) – Unconstrained starting point for the first restart of every structure (the other restarts start at random). bootstrap() uses it to start each refit at the original estimate.

• **options – Fixed options passed to the model constructor.

Returns:

FitResult – Best fit over all structures and restarts.

compare(candidates, data, design=None, restarts=8, rng=None, criterion='bic')

Fit several model classes to the same data and rank them.

Parameters:

• candidates (list) – Model classes, or (class, fit_kwargs) tuples.

• data (dict) – Data passed to fit().

• design (iterable, optional) – Conditions passed to predict.

• restarts (int, optional) – Restarts per fit.

• rng (numpy.random.Generator, optional) – Random number generator.

• criterion ({'bic', 'aic'}, optional) – Ranking criterion (lower is better).

Returns:

list of FitResult – Sorted from best to worst.

recovery(generators, candidates, design, n, reps=20, rng=None, criterion='bic', restarts=8)

Model-recovery study.

Data are simulated from each generating model, every candidate is fitted, and the criterion’s choice is counted. A good design selects the generating model most of the time.

Parameters:

• generators (dict) – {name: model instance} used to simulate data.

• candidates (list) – Model classes (or (class, fit_kwargs)) to compare.

• design (iterable) – Conditions to simulate.

• n (int) – Observations per condition.

• reps (int, optional) – Simulated data sets per generator.

• rng (numpy.random.Generator, optional) – Random number generator.

• criterion ({'bic', 'aic'}, optional) – Selection criterion.

• restarts (int, optional) – Restarts per fit.

Returns:

dict – {generator: {candidate class name: times selected}}.

kl(p, q)

Kullback-Leibler divergence.

Parameters:

• p (array_like) – Probability vectors (clipped at 1e-12).

• q (array_like) – Probability vectors (clipped at 1e-12).

Returns:

float – \(\mathrm{KL}(p\,\Vert\,q) = \sum_i p_i \log(p_i / q_i)\) in nats.

tvd(p, q)

Total variation distance.

Parameters:

• p (array_like) – Probability vectors.

• q (array_like) – Probability vectors.

Returns:

float – \(\tfrac12 \sum_i |p_i - q_i|\), between 0 and 1.

fit_individuals(model_cls, data_by_person, design=None, restarts=4, rng=None, **fit_kwargs)

Fit a model class separately to each person’s data.

Aggregate fitting assumes that everyone shares one set of parameters; per-person fits test that assumption.

Parameters:

• model_cls (type) – Model class.

• data_by_person (dict) – {person: data}.

• design (iterable, optional) – Conditions passed to predict.

• restarts (int, optional) – Restarts per fit.

• rng (numpy.random.Generator, optional) – Random number generator.

• **fit_kwargs – Passed to fit().

Returns:

IndividualFits

Raises:

ValueError – If data_by_person is empty.

compare_individuals(candidates, data_by_person, design=None, restarts=4, rng=None, criterion='bic')

Fit every candidate to every person and compare them.

Parameters:

• candidates (list) – Model classes, or (class, fit_kwargs) tuples.

• data_by_person (dict) – {person: data}.

• design (iterable, optional) – Conditions passed to predict.

• restarts (int, optional) – Restarts per fit.

• rng (numpy.random.Generator, optional) – Random number generator.

• criterion ({'bic', 'aic'}, optional) – Comparison criterion (lower is better).

Returns:

dict – 'group': [(model name, group criterion)] sorted from best; 'best_counts': {model name: number of people it fits best}; 'fits': {model name: IndividualFits}.