bayeshscore

Bayesian model comparison using the Hyvärinen score.
[more details in S. Shao, P.E. Jacob, J. Ding, and V. Tarokh (2017), available here (or alternatively there)]

This package provides functions that compute the Hyvärinen score (and the log-evidence as an aside) of a Bayesian model. This is achieved by using either SMC [e.g. Chopin (2002) and Del Moral, Doucet, Jasra (2006)] or SMC² [cf. Chopin, Jacob, Papaspiliopoulos (2013)] depending on whether the likelihood can be evaluated.

Table of contents

1. Computing the Hyvärinen score
2. Defining a model
3. Setting algorithmic parameters
4. Output description
5. Particle filter
6. Examples

Computing the Hyvärinen score

After installing and loading the package, computing the Hyvärinen score (abbrv. H-score) of a model is done by calling

hscore(observations, model, algorithmic_parameters)

where the inputs consist of

• observations: sequence of dimY-dimensional observations Y₁ , ... , Y_T encoded as a dimY by T matrix
• model: Bayesian model encoded as a list (see Defining a model)
• algorithmic_parameters : algorithmic parameters provided as a list (see Setting algorithmic parameters)

The function hscore is a wrapper that either calls smc(observations, model, algorithmic_parameters) or smc2(observations, model, algorithmic_parameters) depending on whether the model specifies the likelihood or not.
The output is detailed below (see Output description).

Defining a model

The model needs to be provided as a list. A description of the required fields is provided in inst/_template_model.R. Missing fields may be automatically filled-in via the function set_default_model defined in R/util_default.R.

Setting algorithmic parameters

The algorithmic_parameters need to be provided as a list. A complete description of the required fields along with their default values can be found in the function set_default_algorithmic_parameters defined in R/util_default.R.

Output description

The output of hscore, smc, or smc2 is a list. Depending on the specified algorithmic_parameters, some of its fields may be set to NULL. In its most exhaustive form, the output will contain the following objects:

• thetas : last set of particles thetas (dimtheta by Ntheta matrix)
• normw : corresponding normalized weights (vector of length Ntheta)
• byproducts or PFs : corresponding list of byproducts (e.g. particle filters for smc2, see Particle filter)
• logtargetdensities : corresponding evaluations of target log-densities (vector of length Ntheta)
• thetas_history : list of successive sets of particles thetas (one matrix per timestep, starting with the prior, so the length of the list is T+1)
• normw_history : list of successive normalized weights (i.e. list of vector)
• logtargetdensities_history : list of successive target log-densities evaluations (i.e. list of vector)
• byproducts_history of PF_history : list of successive byproducts or particle filters (i.e. list of list)
• logevidence : cumulative logevidence (vector of length T)
• hscore : cumulative H-score using Fisher/Louis type identities (vector of length T)
• hscoreDDE : cumulative H-score using kernel density estimation (vector of length T)
• ESS : successive ESS (vector of length T)
• rejuvenation_times : successive times when rejuvenation occured (vector of random length less than T)
• rejuvenation_rate : associated acceptance rates (vector of random length less than T)
• method : method called (string equal to either "SMC" or "SMC2")
• algorithmic_parameters : list of algorithmic parameters used

Particle filter

The particle filter is implemented in R/conditional_particle_filter.R. The function conditional_particle_filter uses the bootstrap particle filter. It takes an optional argument path when conditioning on a particular path is needed. The output is a list containing the number of x-particles (Nx), the last set of x-particles (X), their respective normalized weights (xnormW), an estimator of the log-likelihood (log_p_y_hat), the last incremental log-likelihood (incremental_ll), a tree encoding the paths and ancestral lineages (tree), and a realized path sampled from the final sets of paths (path).

Examples

Some examples are provided in the folder inst. These scripts reproduce the figures presented in S. Shao, P.E. Jacob, J. Ding, and V. Tarokh (2017) and its supplementary material. Open access to the manuscript and its supplement is also available here.