EmpiricalBayes#

class EmpiricalBayes(optimiser=None, objective='minimise', iteration=1000, tolerance=1e-06, chain=4, quiet=False, start='random', **kwargs)[source]#

Bases: object

Implements an Expectation-Maximisation algorithm for the optimisation of the group-level distributions of the parameters of a model from subject-level parameter estimations.

Parameters:

optimiserobject: The initialized Optimiser object. It must use an optimisation algorithm that also returns the Hessian matrix.
objectivestr: The objective of the optimisation, either ‘maximise’ or ‘minimise’. Default is ‘minimise’. Only affects how we arrive at the participant-level a posteriori parameter estimates.
iterationint, optional: The maximum number of iterations. Default is 1000.
tolerancefloat, optional: The tolerance for convergence. Default is 1e-6.
chainint, optional: The number of random parameter initialisations. Default is 4. The first chain starts from the priors of the model, and every later chain from priors drawn at random within the bounds of the parameters; numpy.random.seed makes them reproducible.
quietbool, optional: Whether to suppress the output. Default is False.
start{“random”, “prior_mean”}, optional: Where the participant-wise fits start on every iteration. With “random”, every start is drawn at random within the bounds of the parameters. With “prior_mean”, the first start is the mean (location) of the current prior of each parameter, clipped into its bounds, as the cbm toolbox does, and the others are drawn at random. Default is “random”.

Notes

The EmpiricalBayes class implements an Expectation-Maximisation algorithm for the optimisation of the group-level distributions of the parameters of a model from subject-level parameter estimations. For the complete description of the method, please see Gershman (2016).

The fitting function must return the Hessian matrix of the optimisation. The Hessian matrix is then used in establishing the within-subject variance of the parameters. It is also important to note that we will require the Hessian matrix of second derivatives of the negative log posterior (Gershman, 2016, p. 3). This requires us to minimise or maximise the log posterior density as opposed to a simple log likelihood, when estimating participant-level parameters.

In the current implementation, we try to calculate the second derivative of the negative log posterior density function according to the following algorithm:

Attempt to use Cholesky decomposition.
If fails, attempt to use LU decomposition.
If fails, attempt to use QR decomposition.
If the result is a complex number with zero imaginary part, keep the real part.

In addition, because the the Hessian matrix should correspond to the precision matrix, hence its inverse is the variance-covariance matrix, we will use its inverse to calculate the within-subject variance of the parameters. If the algorithm fails to calculate the inverse of the Hessian matrix, it will use the Moore-Penrose pseudoinverse instead.

The current implementation also controls for some edge-cases that are not covered by the algorithm above:

When calculating the within-subject variance via the Hessian matrix, the algorithm clips the variance to a minimum value of 1e-6 to avoid numerical instability.
When calculating the within-subject variance via the Hessian matrix, the algorithm sets any non-finite or non-positive values to NaN.
If the second derivative of the negative log posterior density function is not finite, we set the log determinant to -1e6.

An estimate on a bound of a parameter can have a Hessian matrix that is not positive definite. The estimate is kept, as the bounds require. The positive_definite column of fit records whether each participant’s Hessian matrix was positive definite on each iteration, and the not_positive_definite column of hyperparameters how many were not.

References

Gershman, S. J. (2016). Empirical priors for reinforcement learning models. Journal of Mathematical Psychology, 71, 1-6.

Examples

>>> from cpm.applications.reinforcement_learning import RLRW
>>> from cpm.datasets import load_bandit_data
>>> from cpm.hierarchical import EmpiricalBayes
>>> from cpm.optimisation import FminBound, minimise
>>> data = load_bandit_data()
>>> data["observed"] = data["response"]
>>> model = RLRW(data=data[data.ppt == 1], dimensions=4)
>>> optimiser = FminBound(
...     model=model,
...     data=data.groupby("ppt"),
...     initial_guess=None,
...     number_of_starts=2,
...     minimisation=minimise.LogLikelihood.bernoulli,
...     parallel=False,
...     prior=True,
...     ppt_identifier="ppt",
...     display=False,
...     maxiter=200,
...     approx_grad=True,
... )
>>> eb = EmpiricalBayes(optimiser=optimiser, iteration=1000, tolerance=1e-6, chain=4)
>>> eb.optimise()

stair(chain_index=0)[source]#

The main function that runs the Expectation-Maximisation algorithm for the optimisation of the group-level distributions of the parameters of a model from subject-level parameter estimations. This is essentially a single chain.

Returns:

dict: A dictionary containing the log model evidence, the hyperparameters of the group-level distributions, and the parameters of the model.

optimise()[source]#: This method runs the Expectation-Maximisation algorithm for the optimisation of the group-level distributions of the parameters of a model from subject-level parameter estimations. This is essentially the main function that runs the algorithm for multiple chains with random starting points for the priors.

parameters()[source]#

Returns the estimated individual-level parameters for each iteration and chain.

Returns:

pandas.DataFrame: The estimated individual-level parameters for each iteration and chain.

diagnostics(show=True, save=False, path=None)[source]#

Returns the convergence diagnostics plots for the group-level hyperparameters.

Parameters:

showbool, optional: Whether to show the plots. Default is True.
savebool, optional: Whether to save the plots. Default is False.
pathstr, optional: The path to save the plots. Default is None.

Notes

The convergence diagnostics plots show the convergence of the log model evidence, the means, and the standard deviations of the group-level hyperparameters. It also shows the distribution of the means and the standard deviations of the group-level hyperparameters sampled for each chain.