VariationalBayes#

class VariationalBayes(optimiser=None, objective='minimise', iteration=50, tolerance_lme=0.001, tolerance_param=0.001, chain=4, hyperpriors=None, convergence='parameters', quiet=False, start='random', **kwargs)[source]#

Bases: object

Performs hierarchical Bayesian estimation of a given model using variational (approximate) inference methods, a reduced version of the Hierarchical Bayesian Inference (HBI) algorithm proposed by Piray et al. (2019), to exclude model comparison and selection.

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 50.
tolerance_lmefloat, optional: The tolerance for convergence with respect to the log model evidence. Default is 1e-3.
tolerance_paramfloat, optional: The tolerance for convergence with respect to the “normalized” means of parameters. Default is 1e-3.
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.
hyperpriors: dict, optional: A dictionary of given parameter values of the prior distributions on the population-level parameters (means mu and precisions tau). See Notes for details. Default is None.
convergencestr, optional: The convergence criterion. Default is ‘parameters’. Options are ‘lme’ and ‘parameters’.
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 in the HBI implementation of Piray et al. (2019), and the others are drawn at random. Default is “random”.

Notes

The hyperprios are as follows:

a0array-like
Vector of means of the normal prior on the population-level means, mu.
bfloat
Scalar value that is multiplied with population-level precisions, tau, to determine the standard deviations of the normal prior on the population-level means, mu.
vfloat
Scalar value that is used to determine the shape parameter (nu) of the gamma prior on population-level precisions, tau.
sarray-like
Vector of values that serve as lower bounds on the scale parameters (sigma) of the gamma prior on population-level precisions, tau.

With the number of parameters as N, the default values are as follows:

a0 : np.zeros(N)
b : 1
v : 0.5
s : np.repeat(0.01, N)

The convergence criterion can be set to ‘lme’ or ‘parameters’. If set to ‘lme’, the algorithm will stop when the log model evidence converges. If set to ‘parameters’, the algorithm will stop when the “normalized” means of the population-level parameters converge.

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

Piray, P., Dezfouli, A., Heskes, T., Frank, M. J., & Daw, N. D. (2019). Hierarchical Bayesian inference for concurrent model fitting and comparison for group studies. PLoS computational biology, 15(6), e1007043.

Examples

>>> from cpm.applications.reinforcement_learning import RLRW
>>> from cpm.datasets import load_bandit_data
>>> from cpm.hierarchical import VariationalBayes
>>> 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"),
...     minimisation=minimise.LogLikelihood.bernoulli,
...     prior=True,
...     number_of_starts=2,
...     ppt_identifier="ppt",
...     display=False,
... )
>>> vb = VariationalBayes(optimiser=optimiser, iteration=50, chain=4)
>>> vb.optimise()

See the hierarchical estimation tutorials in the documentation for a full walkthrough.

get_lme(log_post, hessian)[source]#

Function to approximate the participant-wise log model evidence using Laplace’s approximation.

Parameters:

log_postarray_like: Participant-wise value of log posterior density function at the mode (i.e., MAP parameter estimates).
hessianarray_like: Participant-wise Hessian matrix of log posterior density function evaluated at the mode (i.e., MAP parameter estimates).

Returns:

lmearray: Participant-wise log model evidence.
lme_sumfloat: Summed log model evidence.

update_population(param, hessian, lme, iter_idx=0, chain_idx=0)[source]#

Function to update the population-level parameters based on the results of participant-wise optimisation.

Parameters:

paramarray_like: Participant-wise parameter estimates.
hessianarray_like: Participant-wise Hessian matrices of the log posterior density function evaluated at the mode (i.e., MAP parameter estimates).
lmefloat: Summed log model evidence.
iter_idxint, optional: The iteration index. Default is 0.
chain_idxint, optional: The chain index. Default is 0.

Returns:

population_updatesdict: Dictionary of updated population-level parameters.
param_snrarray_like: Standardised estimates of population-level means.

check_convergence(lme_new, lme_old, param_snr_new, param_snr_old, iter_idx=0, use_lme=True, use_param=True)[source]#

Function to check if the algorithm has converged.

Parameters:

lme_newfloat: The new log model evidence.
lme_oldfloat: The old log model evidence.
param_snr_newarray_like: The new standardised estimates of population-level means.
param_snr_oldarray_like: The old standardised estimates of population-level means.
iter_idxint, optional: The iteration index. Default is 0.
use_lmebool, optional: Whether to use the log model evidence for checking convergence. Default is True.
use_parambool, optional: Whether to use the standardised estimates of population-level means for checking convergence. Default is True.

Returns:

convergencebool: Whether the algorithm has converged.

run_vb(chain_index=0)[source]#

Run the hierarchical Bayesian inference algorithm.

Parameters:

chain_indexint, optional: The chain index. Default is 0.

Returns:

outputdict: Dictionary of results.

optimise()[source]#: Run the Variational Bayes algorithm for multiple chains.

ttest(null=None)[source]#

Perform a one-sample Student’s t-test on the estimated values of population-level means with respect to given null hypothesis values.

Parameters:

nulldict or pd.DataFrame: The null hypothesis values for the population-level means for each parameters.

Returns:

t_dfpd.DataFrame: The results of the t-test.

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.