bristlecone


Tutorial

Quick Start: Predator-Prey Dynamics

Here we use the classic example of snowshoe hare and lynx predator-prey dynamics, to demonstrate the basic functions of Bristlecone. The full example is in /samples/1. lynx-hare.fsx.

In an F# script or interactive session, first load and open the Bristlecone library.

#load "bristlecone.fsx"
open Bristlecone
open Bristlecone.ModelSystem

Defining the model

A model system is defined as a series of differential equations, parameters, and a likelihood function. The model system is here defined as a pair of ordinary differential equations.

/// Number of snowshoe hares
let dhdt' hare lynx alpha beta =
    alpha * hare - beta * hare * lynx

/// Number of lynx
let dldt' lynx hare delta gamma =
    - gamma * lynx + delta * hare * lynx

The parameters required for the model are defined as a CodedMap<Parameter>. For the predator-prey model, they are defined as:

[ // Natural growth rate of hares in absence of predation
  ShortCode.create "alpha",  Parameter.create  Unconstrained 0.10 0.60
  // Death rate per encounter of hares due to predation
  ShortCode.create "beta",   Parameter.create  Unconstrained 0.001 0.0135
  // Efficiency of turning predated hares into lynx
  ShortCode.create "delta",  Parameter.create  Unconstrained 0.001 0.0135 
  // Natural death rate of lynx in the absence of food
  ShortCode.create "gamma",  Parameter.create  Unconstrained 0.10 0.60
  ] |> Map.ofList

The likelihood function

For this example, we are simply using sum of squares as our measure of goodness of fit. Bristlecone includes sum of squares, and some neagtive log likelihood functions within the Bristlecone.Likelihood module.

ModelLibrary.Likelihood.sumOfSquares ["hare"; "lynx"]

Putting it all together

The above models, parameters, and likelihood function must be tied together in a ModelSystem type. A complete model system for this problem is defined here:

let ``predator-prey`` =

    let dhdt p _ x (e:Environment) =
        dhdt' x 
            (e.[ShortCode.create "lynx"]) 
            (p |> Pool.getEstimate "alpha") 
            (p |> Pool.getEstimate "beta")
    
    let dldt p _ y (e:Environment) =
        dldt' y 
            (e.[ShortCode.create "hare"]) 
            (p |> Pool.getEstimate "delta") 
            (p |> Pool.getEstimate "gamma")

    { Equations =  [ code "hare", dhdt
                     code "lynx", dldt ] |> Map.ofList
      Measures  =  [] |> Map.ofList
      Parameters = [ code "alpha",  parameter PositiveOnly 0.10 0.60
                     code "beta",   parameter PositiveOnly 0.001 0.0135
                     code "delta",  parameter PositiveOnly 0.001 0.0135
                     code "gamma",  parameter PositiveOnly 0.10 0.60
                   ] |> Map.ofList
      Likelihood = ModelLibrary.Likelihood.sumOfSquares ["hare"; "lynx"] }

The intermediate functions dhdt and dldt connect the raw ODEs to the functional signature required by Bristlecone. This feeds the current parameter estimates and variable values into the model.

The modelling protocol is defined by an EstimationEngine, which can be reused for any number of analyses.

let engine = 
    Bristlecone.mkContinuous
    |> Bristlecone.withConditioning RepeatFirstDataPoint

Does my model and method work?

It is recommended that any new model and method combination is tested, to verify that - given a known set of parameters - the same parameters are estimated. Bristlecone can draw many sets of random parameters and test that these can be correctly estimated:

let startValues = 
    [ code "lynx", 30.09
      code "hare", 19.58 ] |> Map.ofList

//``predator-prey`` |> Bristlecone.testModel engine Options.testSeriesLength startValues Options.iterations Options.burn

Fitting to observation data

In the example script, we load data using FSharp.Data, from a CSV file. You can choose to load data using any method. The raw data is parsed into a Map of TimeSeries using the TimeSeries.createVarying function.

type PopulationData = FSharp.Data.CsvProvider<"../samples/data/lynx-hare.csv">

let data = 
    let csv = PopulationData.Load "../samples/data/lynx-hare.csv"
    [ code "hare", 
        TimeSeries.fromObservations (csv.Rows |> Seq.map(fun r -> r.Year, float r.Hare))
      code "lynx", 
        TimeSeries.fromObservations (csv.Rows |> Seq.map(fun r -> r.Year, float r.Lynx)) 
    ] |> Map.ofList

To get model fitting results:

let result = 
    ``predator-prey``
    //|> Bristlecone.fit engine Options.iterations Options.burn data

Time Modes

Bristlecone can run models in either discrete time, or continuous time. When running models in continuous time, an integration function is required:

type TimeMode<'data, 'time> =
| Discrete
| Continuous of Integrate<'data, 'time>

and Integrate<'data,'time> = 'time -> 'time -> 'time -> CodedMap<'data> -> CodedMap<ODE> -> CodedMap<'data[]>

Currently only fixed timesteps are supported, but variable timestep support (e.g. for sediment core data) is planned.

Two integration functions are included:

Solver

Function

Description

Runge-Kutta 4 (MathNet Numerics)

Integration.MathNet.integrate

A fourth-order Runge Kutta method to provide approximate solutions to ODE systems.

Runge-Kutta 547M (Open Solving Library for ODEs - Microsoft Research)

Integration.MsftOslo.integrate

A method based on classic Runge-Kutta, but with automatic error and step size control. See the documentation.

Optimisation

Bristlecone supports optimisation functions that have the following type signature:

type Optimise<'data> = int -> int -> Domain -> ('data[] -> 'data) -> ('data * 'data []) list

There are two optimsation techniques currently built-in:

Method

Function

Description

Amoeba (Nelder-Mead)

Optimisation.Amoeba.solve

A gradient-descent method.

MCMC Random Walk

Optimisation.MonteCarlo.randomWalk

A method based on classic Runge-Kutta, but with automatic error and step size control. See the documentation.

Estimation Engines

To use Bristlecone functions requires a configured EstimationEngine. The easiest way is with the helper functions within the Bristlecone module:

Function

Type

Description

mkContinuous

EstimationEngine

Default continuous engine

mkDiscrete

EstimationEngine

Default discrete model engine

withContinuousTime

t : Integrate<'a,'b> -> engine: EstimationEngine<'a,'b> -> EstimationEngine<'a,'b>

Transforms engine to continuous time mode, using the given integrator function.

withConditioning

c: Conditioning -> engine: EstimationEngine<'a,'b> -> EstimationEngine<'a,'b>

Choose how the start point is chosen when solving the model system.

*

namespace Bristlecone
Multiple items
module Bristlecone from Bristlecone

--------------------
namespace Bristlecone
module ModelSystem from Bristlecone
<summary> Represents an ordinary differential equation model system and its likelihood as as objective function that may be optimised. </summary>
val dhdt' : hare:float -> lynx:float -> alpha:float -> beta:float -> float
 Number of snowshoe hares
val hare : float
val lynx : float
val alpha : float
val beta : float
val dldt' : lynx:float -> hare:float -> delta:float -> gamma:float -> float
 Number of lynx
val delta : float
val gamma : float
module ShortCode from Bristlecone
val create : str:string -> ShortCode.ShortCode option
module Parameter from Bristlecone
val create : con:Parameter.Constraint -> bound1:float -> bound2:float -> Parameter.Parameter option
<summary> Create an estimatable `Parameter` that may be constrained or unconstrained to a certain value range. Bristlecone draws an initial value from the given bounds. To retrieve the estimated parameter value, use `Parameter.finalise`. </summary>
Multiple items
module Map from Bristlecone

--------------------
module Map from Microsoft.FSharp.Collections
<summary>Contains operations for working with values of type <see cref="T:Microsoft.FSharp.Collections.Map`2" />.</summary>

--------------------
type Map<'Key,'Value (requires comparison)> = interface IReadOnlyDictionary<'Key,'Value> interface IReadOnlyCollection<KeyValuePair<'Key,'Value>> interface IEnumerable interface IComparable interface IEnumerable<KeyValuePair<'Key,'Value>> interface ICollection<KeyValuePair<'Key,'Value>> interface IDictionary<'Key,'Value> new : elements:seq<'Key * 'Value> -> Map<'Key,'Value> member Add : key:'Key * value:'Value -> Map<'Key,'Value> member Change : key:'Key * f:('Value option -> 'Value option) -> Map<'Key,'Value> ...
<summary>Immutable maps based on binary trees, where keys are ordered by F# generic comparison. By default comparison is the F# structural comparison function or uses implementations of the IComparable interface on key values.</summary>
<remarks>See the <see cref="T:Microsoft.FSharp.Collections.MapModule" /> module for further operations on maps. All members of this class are thread-safe and may be used concurrently from multiple threads.</remarks>


--------------------
new : elements:seq<'Key * 'Value> -> Map<'Key,'Value>
val ofList : elements:('Key * 'T) list -> Map<'Key,'T> (requires comparison)
<summary>Returns a new map made from the given bindings.</summary>
<param name="elements">The input list of key/value pairs.</param>
<returns>The resulting map.</returns>
module ModelLibrary
<summary> Pre-built model parts for use in Bristlecone. </summary>
module Likelihood from ModelLibrary
<summary> Likelihood functions to represent a variety of distributions and data types. </summary>
val sumOfSquares : keys:string list -> 'a -> data:CodedMap<PredictedSeries> -> float
val ( predator-prey ) : ModelSystem
val dhdt : ('a -> 'b -> float -> Environment -> float)
val p : 'a
val x : float
val e : Environment
type Environment = CodedMap<float>
val dldt : ('a -> 'b -> float -> Environment -> float)
val y : float
type Likelihood = Parameter.Pool -> CodedMap<PredictedSeries> -> float
<summary> A function that computes the likelihood of a set of parameters. </summary>
val engine : EstimationEngine.EstimationEngine<float,float>
val mkContinuous : EstimationEngine.EstimationEngine<float,float>
<summary> A standard `EstimationEngine` for ordinary differential equation models. </summary>
val withConditioning : c:Conditioning.Conditioning<'a> -> engine:EstimationEngine.EstimationEngine<'a,'b> -> EstimationEngine.EstimationEngine<'a,'b>
<summary> Choose how the start point is chosen when solving the model system </summary>
val startValues : Map<System.IComparable,float>
type PopulationData = obj
namespace Microsoft.FSharp
namespace Microsoft.FSharp.Data
val data : Map<System.IComparable,obj>
val csv : obj
Multiple items
module Seq from Bristlecone

--------------------
module Seq from Microsoft.FSharp.Collections
<summary>Contains operations for working with values of type <see cref="T:Microsoft.FSharp.Collections.seq`1" />.</summary>
val map : mapping:('T -> 'U) -> source:seq<'T> -> seq<'U>
<summary>Builds a new collection whose elements are the results of applying the given function to each of the elements of the collection. The given function will be applied as elements are demanded using the <c>MoveNext</c> method on enumerators retrieved from the object.</summary>
<remarks>The returned sequence may be passed between threads safely. However, individual IEnumerator values generated from the returned sequence should not be accessed concurrently.</remarks>
<param name="mapping">A function to transform items from the input sequence.</param>
<param name="source">The input sequence.</param>
<returns>The result sequence.</returns>
<exception cref="T:System.ArgumentNullException">Thrown when the input sequence is null.</exception>
Multiple items
val float : value:'T -> float (requires member op_Explicit)
<summary>Converts the argument to 64-bit float. This is a direct conversion for all primitive numeric types. For strings, the input is converted using <c>Double.Parse()</c> with InvariantCulture settings. Otherwise the operation requires an appropriate static conversion method on the input type.</summary>
<param name="value">The input value.</param>
<returns>The converted float</returns>


--------------------
[<Struct>] type float = System.Double
<summary>An abbreviation for the CLI type <see cref="T:System.Double" />.</summary>
<category>Basic Types</category>


--------------------
type float<'Measure> = float
<summary>The type of double-precision floating point numbers, annotated with a unit of measure. The unit of measure is erased in compiled code and when values of this type are analyzed using reflection. The type is representationally equivalent to <see cref="T:System.Double" />.</summary>
<category index="6">Basic Types with Units of Measure</category>
val result : ModelSystem
type TimeMode<'data,'time> = | Discrete | Continuous of Integrate<'data,'time>
union case TimeMode.Discrete: TimeMode<'data,'time>
union case TimeMode.Continuous: Integrate<'data,'time> -> TimeMode<'data,'time>
type Integrate<'data,'time> = 'time -> 'time -> 'time -> CodedMap<'data> -> obj -> CodedMap<'data []>
type CodedMap<'T> = Map<ShortCode.ShortCode,'T>
type Optimise<'data> = int -> int -> obj -> ('data [] -> 'data) -> ('data * 'data []) list
Multiple items
val int : value:'T -> int (requires member op_Explicit)
<summary>Converts the argument to signed 32-bit integer. This is a direct conversion for all primitive numeric types. For strings, the input is converted using <c>Int32.Parse()</c> with InvariantCulture settings. Otherwise the operation requires an appropriate static conversion method on the input type.</summary>
<param name="value">The input value.</param>
<returns>The converted int</returns>


--------------------
[<Struct>] type int = int32
<summary>An abbreviation for the CLI type <see cref="T:System.Int32" />.</summary>
<category>Basic Types</category>


--------------------
type int<'Measure> = int
<summary>The type of 32-bit signed integer numbers, annotated with a unit of measure. The unit of measure is erased in compiled code and when values of this type are analyzed using reflection. The type is representationally equivalent to <see cref="T:System.Int32" />.</summary>
<category>Basic Types with Units of Measure</category>
type 'T list = List<'T>
<summary>The type of immutable singly-linked lists. </summary>
<remarks>See the <see cref="T:Microsoft.FSharp.Collections.ListModule" /> module for further operations related to lists. Use the constructors <c>[]</c> and <c>::</c> (infix) to create values of this type, or the notation <c>[1; 2; 3]</c>. Use the values in the <c>List</c> module to manipulate values of this type, or pattern match against the values directly. See also <a href="https://docs.microsoft.com/dotnet/fsharp/language-reference/lists">F# Language Guide - Lists</a>. </remarks>