pytransc.analysis

 12def run_laplace_evidence_approximation(
 13    n_states: int,
 14    n_dims: list[int],
 15    log_posterior: MultiStateDensity | None,
 16    map_models,
 17    ensemble_per_state=None,
 18    log_posterior_ens=None,
 19    verbose=False,
 20    optimize=False,
 21):
 22    """
 23    Function to perform Laplace integration for evidence approximation within each state, using either an input log-posterior function, or posterior ensembles.
 24
 25    Parameters
 26    ----------
 27    log_posterior : function    : Supplied function evaluating the log-posterior function
 28                                    up to a multiplicative constant, for each state.
 29                                    (Not used if ensemble_per_state and log_posterior_ens lists are provided)
 30                                    Calling sequence log_posterior(x,state)
 31    map_models - floats         : List of MAP models in each state where Laplace approximation is evaluated.
 32                                    If optimize=True and a log_posterior() function is supplied, then
 33                                    scipy.minimize is used to find MAP models in each state using map_models as starting guesses.
 34    ensemble_per_state - floats : Optional list of posterior samples in each state, format [i][j][k],(i=1,...,n_samples;j=1,..., n_models;k=1,...,ndim[i]).
 35    log_posterior_ens - floats  : Optional list of log-posterior densities of samples in each state, format [i][j],(i=1,...,n_states;j=1,..., n_samples).
 36    optimize, bool              : Logical to decide whether to use optimization for MAP models (Only relevant if log_posterior()) function supplied.)
 37
 38    Returns:
 39    -------
 40    laplace_log_marginal_likelihoods - n_states*float : list of log-evidence/marginal Likelihoods for each state.
 41    laplace_map_models_per_state - n_states*floats : list of updated MAP models for each state.m if optimize=True.
 42    laplace_map_log_posteriors - n_states*float : list of log-posteriors at MAP models for each state.
 43    laplace_hessians - n_states*NxN : list of negative inverse Hessians if posterior function supplied.
 44
 45    Notes:
 46    Calculates Laplace approximations to evidence integrals in each state. This is equivalent to fitting a Gaussian about the MAP in model
 47    space. Here the Hessian in each state is calculated either by numerical differentiation tools (if a log_posterior function is supplied),
 48    or by taking the covariance of the given ensemble. The MAP model is similarly taken as the maximum of the log-posterior
 49    (if an ensemble is given), or a MAP model in each is estimated by optimization.
 50    Alternatively, if the list `map_models` is given then these are used within each state as MAP models.
 51
 52    This implements equation 10.14 of Raftery (1996) evaluating Laplace's approximation to Bayesian evidence.
 53
 54    Raftery, A.E. (1996) Hypothesis Testing and Model Selection. In: Gilks, W., Richardson, S. and Speigelhalter, D.J., Eds.,
 55    Markov Chain Monte Carlo in Practice, Chapman and Hall, 163-187.
 56    """
 57
 58    if (ensemble_per_state is None) and (log_posterior_ens is not None):
 59        raise InputError(
 60            msg=" In function run_laplace_evidence_approximation: Ensemble probabilities provided as argument without ensemble co-ordinates"
 61        )
 62
 63    if (ensemble_per_state is not None) and (log_posterior_ens is None):
 64        raise InputError(
 65            msg=" In function run_laplace_evidence_approximation: Ensemble co-ordinates provided as argument without ensemble probabilities"
 66        )
 67
 68    if isinstance(ensemble_per_state, list) and isinstance(
 69        log_posterior_ens, list
 70    ):  # we use input ensemble
 71        if verbose:
 72            print(
 73                "run_laplace_evidence_approximation: We are using input ensembles rather than input log_posterior function"
 74            )
 75
 76        hessians, map_models, map_log_posteriors, log_marginal_likelihoods = (
 77            _from_ensemble(
 78                n_dims,
 79                ensemble_per_state,
 80                log_posterior_ens,
 81            )
 82        )
 83    elif log_posterior is not None:
 84        # we are using the supplied log_posterior() function so need Hessian and MAp model
 85        if verbose:
 86            if optimize:
 87                print(
 88                    "run_laplace_evidence_approximation: We are using input log_posterior function with optimization for MAP models"
 89                )
 90            else:
 91                print(
 92                    "run_laplace_evidence_approximation: We are using input log_posterior function with provided MAP models"
 93                )
 94
 95        hessians, map_models, map_log_posteriors, log_marginal_likelihoods = (
 96            _from_log_posterior(
 97                n_states,
 98                n_dims,
 99                log_posterior,
100                map_models,
101                optimize=optimize,
102            )
103        )
104    else:
105        raise InputError(
106            msg="In function run_laplace_evidence_approximation: Either ensemble_per_state or log_posterior must be provided."
107        )
108
109    return (
110        hessians,
111        map_models,
112        map_log_posteriors,
113        log_marginal_likelihoods,
114    )