Fitting models to data#

The optimisers in cpm.optimisation estimate the free parameters of a model for every participant separately. They take a Wrapper, the data grouped by participant, and an objective function, and handle the rest: running the model, comparing its predictions with the data, and collecting the results.

We use the model and data of Models, for 10 participants:

from cpm.optimisation import FminBound, minimise

fit = FminBound(
    model=wrapper,
    data=data.groupby("ppt"),
    minimisation=minimise.LogLikelihood.bernoulli,
    prior=True,
    number_of_starts=2,
    ppt_identifier="ppt",
    parallel=False,
    display=False,
    approx_grad=True,  # passed on to scipy.optimize.fmin_l_bfgs_b
)
fit.optimise()

The arguments that all optimisers share:

model: the Wrapper. The bounds of the free parameters come from their Value definitions, so you do not pass them again.
data: the data, grouped by participant.
minimisation: the objective function, see below.
prior: with True, the optimiser minimises the negative log posterior, the negative log-likelihood plus the negative log prior, so the estimates are maximum a posteriori. With False, they are maximum-likelihood estimates.
number_of_starts: how many random starting points to try per participant. The best fit is kept, which protects against local minima.
ppt_identifier: the column that identifies participants, used to label the results.
parallel: fit participants in parallel, see the how-to guide on Parallelise model fitting.
metrics: extra measures of fit to compute at the estimates, such as PenalisedLikelihoods.BIC.

All other keyword arguments are passed on to the underlying scipy.optimize function.

The results#

export() returns one row per participant:

fit.export().head()

	x_0	x_1	grad_0	grad_1	task	funcalls	nit	hessian_0	hessian_1	hessian_2	hessian_3	log_likelihood	log_prior	ppt	fun
0	0.352648	9.230891	0.000000e+00	0.000000e+00	CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL	36	9	185.612488	0.614965	0.614965	0.321324	-16.243776	-2.880474	1	19.124250
1	0.046537	4.897143	9.947598e-06	7.105427e-07	CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL	51	12	1470.910516	20.002881	20.002881	0.411809	-46.677172	-2.920614	2	49.597787
2	0.453396	2.236194	-7.105427e-07	0.000000e+00	CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL	63	15	64.013221	5.909899	5.909899	2.208668	-44.505584	-1.903203	3	46.408787
3	0.689269	6.363270	7.105427e-07	-3.552714e-07	CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL	36	10	112.864865	-1.114420	-1.114420	0.436221	-27.330719	-1.709999	4	29.040718
4	0.388322	5.253323	-5.684342e-06	-7.105427e-07	CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL	33	8	126.933727	2.787430	2.787430	0.807702	-30.633954	-1.379647	5	32.013601

x_0, x_1, … are the estimates of the free parameters, in the order of free(),
fun is the value of the objective function at the estimates,
log_likelihood and log_prior split it into its two parts when prior=True,
the remaining columns come from the optimiser, such as the gradient and the Hessian.

fit.parameters holds the same estimates as a list with one dictionary per participant, which you can pass straight to a Simulator:

fit.parameters[:3]

[{'alpha': np.float64(0.35264764658046016),
  'temperature': np.float64(9.230890556986454)},
 {'alpha': np.float64(0.046536868534096804),
  'temperature': np.float64(4.897142755490239)},
 {'alpha': np.float64(0.4533962586738872),
  'temperature': np.float64(2.2361935716634216)}]

Objective functions#

cpm.optimisation.minimise has the objective functions. They compare the dependent output of the model with the observed column of the data:

Objective	Use it when the model predicts
`LogLikelihood.bernoulli`	the probability of a binary choice
`LogLikelihood.categorical`	the probabilities of more than two options
`LogLikelihood.continuous`	a continuous response
`Distance`	a value to match, with squared-error measures
`Discrete`	frequencies of discrete outcomes

Choosing an optimiser#

Optimiser	Method	Bounds	Hessian (for `cpm.hierarchical`)
`FminBound`	L-BFGS-B, gradient-based	yes	yes
`Fmin`	Nelder–Mead simplex	no	yes
`Minimize`	any method of `scipy.optimize.minimize()`	passed on to the method	yes
`DifferentialEvolution`	global, population-based	yes	no
`Bads`	Bayesian adaptive direct search	yes	yes

FminBound is a good default: it is fast, respects the bounds, and works with the hierarchical methods. If the likelihood surface has many local minima, try more starting points, or a global optimiser such as DifferentialEvolution. For noisy or expensive objective functions, Bads is often more reliable.

Next steps#

Work through the Tutorials, starting with 1. Build and fit your first model: the two-armed bandit.