Diagnostics

Diagnostic statistics are essential to ensure that the model-fitting and model-selection results are valid and complete.

Convergence Statistics

For Monte-Carlo based optimisation techniques (such as the random walk MCMC or the Filzbach-style MCMC) convergence statistics are useful to assess whether multiple independent chains have identified the same minimum -log likelihood in parameter space, which may be the global minimum.

Bristlecone includes the Gelman-Rubin convergence statistic (Rhat), which can be calculated on a per-parameter basis.

open Bristlecone
open Bristlecone.Diagnostics
open Bristlecone.Data

let saveDirectory = "/some/save/dir"

fun listOfResults ->

    // Calculate Gelman-Rubin statistic for all parameters
    // across multiple model runs
    let stat =
        Convergence.gelmanRubinAll 1000 (fun _ -> "some subject") (fun _ -> "some hypothesis") listOfResults

    // Save results to a single file
    Convergence.save saveDirectory stat

Logging individual model components

Sometimes it is helpful for diagnostic or visualisation purposes to see the values of individual parts of a mathematical model during computation, or at the minimum -log likelihood.

Bristlecone includes the ability to log out the values of individual model components through time. An worked example is given below of how to obtain computed values across the time-series for components within a model system. First, you must set up a model that has a parameter of the IComponentLogger<float>) type:

[ NB TODO: This function must be refactored to work with the new v2 Bristlecone Language, and as such is currently disabled. ]

// open Bristlecone
// open Bristlecone.Language
// open Bristlecone.Diagnostics.ModelComponents

// let hypothesis (cLog:IComponentLogger<ModelExpression>) =

//     let ``dN/dt`` =
//         Parameter "r" * This * cLog.StoreValue "log this thing" (Constant 1. - (This / Parameter "K"))

//     Model.empty
//     |> Model.addEquation       "N"      ``dN/dt``
//     |> Model.estimateParameter "r"      noConstraints 0.10 5.00
//     |> Model.estimateParameter "K"      noConstraints 50.0 150.
//     |> Model.useLikelihoodFunction      (ModelLibrary.Likelihood.sumOfSquares [ "N" ])
//     |> Model.compile


// fun engine singleResult ->
//     // A wrapper around the Bristlecone.fit function
//     let fitFn ts m e = Bristlecone.fit e (Optimisation.EndConditions.afterIteration 1) ts m
//     let logs = ModelComponents.calculateComponents fitFn engine singleResult
//     ()

Multiple items
module Bristlecone from Bristlecone
<namespacedoc><summary>The core library of Bristlecone, containing model-fitting functions.</summary></namespacedoc>
Main functionality of Bristlecone, including functions to scaffold `ModelSystem`s and for model-fitting (tests and real fits).

--------------------
namespace Bristlecone

module Diagnostics from Bristlecone
<summary>Diagnostic techniques for determining the suitability of results obtained with Bristlecone.</summary>

namespace Bristlecone.Data

val saveDirectory: string

val listOfResults: ModelSelection.ResultSet.ResultSet<'a,'b> seq

val stat: Convergence.ConvergenceStatistic seq

Multiple items
module Convergence from Bristlecone.Data

--------------------
module Convergence from Bristlecone.Diagnostics
<summary>Convergence diagnostics for monte-carlo markov chain (MCMC) analyses.</summary>

val gelmanRubinAll: nMostRecent: int -> subject: ('subject -> string) -> hypothesis: ('hypothesis -> string) -> results: ModelSelection.ResultSet.ResultSet<'subject,'hypothesis> seq -> Convergence.ConvergenceStatistic seq
<summary>Calculate the Gelman-Rubin statistic for each parameter in all of the given `ResultSet`. The statistic tends downwards to one, with one indicating perfect convergence between all chains.</summary>
<param name="nMostRecent">How many recent iterations to use from the trace.</param>
<param name="subject">A function to retrieve a subject ID from a subject</param>
<param name="hypothesis">A function to retrieve a hypothesis ID from a hypothesis</param>
<param name="result">A result set (of 1 .. many results) for a particular subject and hypothesis</param>
<returns>If more than one replicate, the R-hat convergence statistic across replicates</returns>

val save: directory: string -> result: Convergence.ConvergenceStatistic seq -> unit