Bristlecone has been designed to enable swift model composition, whereby
smaller models may be composed into larger systems of models.
To get started, open the Bristlecone libraries.
open Bristlecone // Opens Bristlecone core library and estimation engine
open Bristlecone.Language // Open the language for writing Bristlecone models
open FSharp.Data.UnitSystems.SI.UnitNames
One key use of compositional modelling is to understand which (ecological) theories
best represent one or many particular time-series data. A nested model is one where
one or more processes may be substituted for different theoretical forms and competed
to determine which combination of model components best represents the given data.
Using Bristlecone, we can simply define a nested model as a Bristlecone model but where
the definition takes some additional model components as function parameters.
We can define a nestable model, in which a number of components
may be swapped for alternative mathematical forms based on (ecological) theory. Here, we will
apply a simple model of non-linear plant growth. To do this, we
first define a base model as a ModelSystem that takes a number of interchangable components:
let r = parameter "r" noConstraints 0.01< / Time.month> 1.00< / Time.month>
let m = state<kilogram> "m"
let baseModel growthLimit lossRate =
Model.empty<Time.month>
|> Model.addRateEquation m (P r * growthLimit This<kilogram> - lossRate This<kilogram>)
|> Model.estimateParameter r
|> Model.useLikelihoodFunction (ModelLibrary.Likelihood.sumOfSquares [ Require.state m ])
You can see from the type signatures that growthLimit returns a kilogram value, whereas
the loss rate returns a rate (kilogram per month).
Importantly, unlike when normally building models, we do not call Model.compile at the end.
This step will be done later.
We can then design the growth limit and loss rate components that will fit into this model.
These are made using the modelComponent and subComponent functions as follows:
open Bristlecone.Language.Components
let K = parameter "K" notNegative 10.0<kilogram> 50.0<kilogram>
let growthLimit =
modelComponent
"Limitation to growth"
[
subComponent "Linear" (fun _ -> Constant 1.<kilogram>)
subComponent "Monomolecular" (fun (m: ModelExpression<kilogram>) -> P K - m)
|> estimateParameter K
subComponent "Gompertz" (fun (m: ModelExpression<kilogram>) -> m * Logarithm(P K / m))
|> estimateParameter K
subComponent "Logistic (three parameter)" (fun (m: ModelExpression<kilogram>) -> m * (Constant 1. - m / P K))
|> estimateParameter K ]
If a component requires additional parameters over the base model, these may be added by piping into
the estimateParameter function as above.
We can similarly define the other required component for the base model - the loss rate:
let alpha = parameter "alpha" notNegative 0.001< / Time.month> 0.10< / Time.month>
let lossRate =
modelComponent
"Biomass loss rate"
[
subComponent "None" (fun _ -> Constant 0.<kilogram / Time.month>)
subComponent "Density-dependent" (fun m -> m * P alpha)
|> estimateParameter alpha ]
Once the components are set up, we can compile the nested model by making all possible
combinations of the model components. We do this using the Hypotheses module, which will
compile a list of model hypotheses for us to test.
let hypotheses =
baseModel
|> Hypotheses.createFromModel
|> Hypotheses.apply growthLimit
|> Hypotheses.apply lossRate
|> Hypotheses.compile
The resultant constructed hypotheses are shown in the following printout:
[Hypothesis
({ Parameters =
Pool
(map
[(ShortCode "r",
{ Name = "r"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> })])
EnvironmentKeys = []
Equations =
DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
Measures = map []
Initialisers = map []
NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
Evaluate = <fun:Pipe #1 input at line 27@27> } },
[{ Component = "Limitation to growth"
Implementation = "Linear" }; { Component = "Biomass loss rate"
Implementation = "None" }]);
Hypothesis
({ Parameters =
Pool
(map
[(ShortCode "K",
{ Name = "K"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> });
(ShortCode "r",
{ Name = "r"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> })])
EnvironmentKeys = []
Equations =
DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
Measures = map []
Initialisers = map []
NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
Evaluate = <fun:Pipe #1 input at line 27@27> } },
[{ Component = "Limitation to growth"
Implementation = "Monomolecular" }; { Component = "Biomass loss rate"
Implementation = "None" }]);
Hypothesis
({ Parameters =
Pool
(map
[(ShortCode "K",
{ Name = "K"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> });
(ShortCode "r",
{ Name = "r"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> })])
EnvironmentKeys = []
Equations =
DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
Measures = map []
Initialisers = map []
NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
Evaluate = <fun:Pipe #1 input at line 27@27> } },
[{ Component = "Limitation to growth"
Implementation = "Gompertz" }; { Component = "Biomass loss rate"
Implementation = "None" }]);
Hypothesis
({ Parameters =
Pool
(map
[(ShortCode "K",
{ Name = "K"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> });
(ShortCode "r",
{ Name = "r"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> })])
EnvironmentKeys = []
Equations =
DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
Measures = map []
Initialisers = map []
NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
Evaluate = <fun:Pipe #1 input at line 27@27> } },
[{ Component = "Limitation to growth"
Implementation = "Logistic (three parameter)" };
{ Component = "Biomass loss rate"
Implementation = "None" }]);
Hypothesis
({ Parameters =
Pool
(map
[(ShortCode "alpha",
{ Name = "alpha"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> });
(ShortCode "r",
{ Name = "r"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> })])
EnvironmentKeys = []
Equations =
DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
Measures = map []
Initialisers = map []
NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
Evaluate = <fun:Pipe #1 input at line 27@27> } },
[{ Component = "Limitation to growth"
Implementation = "Linear" }; { Component = "Biomass loss rate"
Implementation = "Density-dependent" }]);
Hypothesis
({ Parameters =
Pool
(map
[(ShortCode "K",
{ Name = "K"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> });
(ShortCode "alpha",
{ Name = "alpha"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> });
(ShortCode "r",
{ Name = "r"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> })])
EnvironmentKeys = []
Equations =
DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
Measures = map []
Initialisers = map []
NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
Evaluate = <fun:Pipe #1 input at line 27@27> } },
[{ Component = "Limitation to growth"
Implementation = "Monomolecular" };
{ Component = "Biomass loss rate"
Implementation = "Density-dependent" }]);
Hypothesis
({ Parameters =
Pool
(map
[(ShortCode "K",
{ Name = "K"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> });
(ShortCode "alpha",
{ Name = "alpha"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> });
(ShortCode "r",
{ Name = "r"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> })])
EnvironmentKeys = []
Equations =
DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
Measures = map []
Initialisers = map []
NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
Evaluate = <fun:Pipe #1 input at line 27@27> } },
[{ Component = "Limitation to growth"
Implementation = "Gompertz" }; { Component = "Biomass loss rate"
Implementation = "Density-dependent" }]);
Hypothesis
({ Parameters =
Pool
(map
[(ShortCode "K",
{ Name = "K"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> });
(ShortCode "alpha",
{ Name = "alpha"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> });
(ShortCode "r",
{ Name = "r"
ToTensorRealIO = <fun:make@119-1>
FromTensorRealIO = <fun:make@124-2>
GetConstraint = <fun:make@130-3>
TryGetReal = <fun:make@131-4>
TryGetBounds = <fun:make@133-6> })])
EnvironmentKeys = []
Equations =
DifferentialEqs (map [(ShortCode "m", <fun:ToFSharpFunc@3936-2>)])
Measures = map []
Initialisers = map []
NegLogLikelihood = { RequiredCodes = [State ShortCode "m"]
Evaluate = <fun:Pipe #1 input at line 27@27> } },
[{ Component = "Limitation to growth"
Implementation = "Logistic (three parameter)" };
{ Component = "Biomass loss rate"
Implementation = "Density-dependent" }])]
|
For further analysis, you may choose to use the orchestration functionality to
run many hypotheses and replicates (multiple fits per hypothesis) in parallel.
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 Language
from Bristlecone
<summary>
An F# Domain Specific Language (DSL) for scripting with
Bristlecone.
</summary>
namespace Microsoft.FSharp
namespace Microsoft.FSharp.Data
namespace Microsoft.FSharp.Data.UnitSystems
namespace Microsoft.FSharp.Data.UnitSystems.SI
namespace Microsoft.FSharp.Data.UnitSystems.SI.UnitNames
val r: IncludedParameter</Time.month>
Multiple items
val parameter: code: string -> con: Parameter.Constraint -> lower: float<'u> -> upper: float<'u> -> IncludedParameter<'u>
<summary>
Define an estimatable parameter for a Bristlecone model.
</summary>
--------------------
[<Measure>]
type parameter
val noConstraints: Parameter.Constraint
Multiple items
val Time<'t> : ModelExpression<'t>
--------------------
module Time
from Bristlecone
[<Measure>]
type month
val m: StateId<kilogram>
val state: name: string -> StateId<'u>
[<Measure>]
type kilogram
val baseModel: growthLimit: (ModelExpression<kilogram> -> ModelExpression<kilogram>) -> lossRate: (ModelExpression<kilogram> -> ModelExpression<kilogram/Time.month>) -> ModelBuilder.ModelBuilder<Time.month>
val growthLimit: (ModelExpression<kilogram> -> ModelExpression<kilogram>)
val lossRate: (ModelExpression<kilogram> -> ModelExpression<kilogram/Time.month>)
module Model
from Bristlecone.Language
<summary>
Terms for scaffolding a model system for use with Bristlecone.
</summary>
val empty<'time> : ModelBuilder.ModelBuilder<'time>
val addRateEquation: name: StateId<'state> -> expr: ModelExpression<'state/'time> -> mb: ModelBuilder.ModelBuilder<'time> -> ModelBuilder.ModelBuilder<'time>
val P: name: IncludedParameter<'u> -> ModelExpression<'u>
val This<'u> : ModelExpression<'u>
val estimateParameter: p: IncludedParameter<'u> -> builder: ModelBuilder.ModelBuilder<'time> -> ModelBuilder.ModelBuilder<'time>
val useLikelihoodFunction: likelihoodFn: ModelSystem.Likelihood<ModelSystem.state> -> builder: ModelBuilder.ModelBuilder<'u> -> ModelBuilder.ModelBuilder<'u>
namespace Bristlecone.ModelLibrary
module Likelihood
from Bristlecone.ModelLibrary
<summary>Likelihood functions to represent a variety of distributions and data types.</summary>
<namespacedoc><summary>Pre-built model parts for use in Bristlecone</summary></namespacedoc>
val sumOfSquares: keys: ModelSystem.LikelihoodRequirement list -> ModelSystem.Likelihood<'u>
<summary>
Residual sum of squares. Provides a simple metric of distance between
observed data and model predictions.
</summary>
module Require
from Bristlecone.Language
val state: s: StateId<'u> -> ModelSystem.LikelihoodRequirement
module Components
from Bristlecone.Language
val K: IncludedParameter<kilogram>
val notNegative: Parameter.Constraint
val growthLimit: ModelComponent<(ModelExpression<kilogram> -> ModelExpression<kilogram>)>
val modelComponent: name: string -> list: SubComponent<'a> list -> ModelComponent<'a>
val subComponent: name: string -> transform: 'a -> SubComponent<'a>
<summary>
Creates a nested component that can be inserted into a base model.
</summary>
val Constant: x: float<'u> -> ModelExpression<'u>
val m: ModelExpression<kilogram>
type ModelExpression<'u> =
private | ME of ModelExpressionUntyped
static member (%) : ModelExpression<'u> * ModelExpression<'u> -> ModelExpression<'u>
static member ( * ) : ModelExpression<'u1> * ModelExpression<'u2> -> ModelExpression<'u1 'u2>
static member (+) : ModelExpression<'u> * ModelExpression<'u> -> ModelExpression<'u>
static member (-) : ModelExpression<'u> * ModelExpression<'u> -> ModelExpression<'u>
static member (.<) : ModelExpression<'u> * ModelExpression<'u> -> BoolExpression
static member (.>) : ModelExpression<'u> * ModelExpression<'u> -> BoolExpression
static member (/) : ModelExpression<'u1> * ModelExpression<'u2> -> ModelExpression<'u1/'u2>
static member (=.) : ModelExpression<'u> * ModelExpression<'u> -> BoolExpression
static member Environment: sid: StateId<'u> -> ModelExpression<'u>
static member Exp: ModelExpression<1> -> ModelExpression<1>
...
val estimateParameter: p: IncludedParameter<'u> -> comp: SubComponent<'a> -> SubComponent<'a>
val Logarithm: expr: ModelExpression<1> -> ModelExpression<1>
val alpha: IncludedParameter</Time.month>
val lossRate: ModelComponent<(ModelExpression<kilogram> -> ModelExpression<kilogram/Time.month>)>
val hypotheses: Hypotheses.Hypothesis<Time.month> list
module Hypotheses
from Bristlecone.Language
<summary>Types to represent a hypothesis, given that a hypothesis
is a model system that contains some alternate formulations of certain
components.</summary>
val createFromModel: baseModel: ('a -> 'rest) -> Writer.Writer<('a -> 'rest),(Hypotheses.ComponentName * CodedMap<Parameter.Pool.AnyParameter>)> list
val apply: comp: ModelComponent<'a> -> builders: Writer.Writer<('a -> 'rest),(Hypotheses.ComponentName * CodedMap<Parameter.Pool.AnyParameter>)> list -> Writer.Writer<'rest,(Hypotheses.ComponentName * CodedMap<Parameter.Pool.AnyParameter>)> list
val compile: builders: Writer.Writer<ModelBuilder.ModelBuilder<'u>,(Hypotheses.ComponentName * CodedMap<Parameter.Pool.AnyParameter>)> list -> Hypotheses.Hypothesis<'u> list
<summary>
Compile: run the writer(s), add parameters into the model builder, and wrap in Hypothesis
</summary>