Uncertainty Quantification

The following callbacks are designed to help enable uncertainty quantification of the differential equation solutions.

DiffEqCallbacks.ProbIntsUncertaintyFunction
ProbIntsUncertainty(σ, order, save = true)

The ProbInts method for uncertainty quantification involves the transformation of an ODE into an associated SDE where the noise is related to the timesteps and the order of the algorithm.

Arguments

  • σ is the noise scaling factor. It is recommended that σ is representative of the size of the errors in a single step of the equation. If such a value is unknown, it can be estimated automatically in adaptive time-stepping algorithms via AdaptiveProbIntsUncertainty
  • order is the order of the ODE solver algorithm.
  • save is for choosing whether this callback should control the saving behavior. Generally this is true unless one is stacking callbacks in a CallbackSet.

References

Conrad P., Girolami M., Särkkä S., Stuart A., Zygalakis. K, Probability Measures for Numerical Solutions of Differential Equations, arXiv:1506.04592

DiffEqCallbacks.AdaptiveProbIntsUncertaintyFunction
AdaptiveProbIntsUncertainty(order, save = true)

The ProbInts method for uncertainty quantification involves the transformation of an ODE into an associated SDE where the noise is related to the timesteps and the order of the algorithm.

AdaptiveProbIntsUncertainty is a more automated form of ProbIntsUncertainty which uses the error estimate from within adaptive time stepping methods to estimate σ at every step.

Arguments

  • order is the order of the ODE solver algorithm.
  • save is for choosing whether this callback should control the saving behavior. Generally this is true unless one is stacking callbacks in a CallbackSet.

References

Conrad P., Girolami M., Särkkä S., Stuart A., Zygalakis. K, Probability Measures for Numerical Solutions of Differential Equations, arXiv:1506.04592

Example 1: FitzHugh-Nagumo

In this example we will determine our uncertainty when solving the FitzHugh-Nagumo model with the Euler() method. We define the FitzHugh-Nagumo model:

using DiffEqCallbacks, OrdinaryDiffEq, Plots

function fitz(du,u,p,t)
  V,R = u
  a,b,c = p
  du[1] = c*(V - V^3/3 + R)
  du[2] = -(1/c)*(V -  a - b*R)
end
u0 = [-1.0;1.0]
tspan = (0.0,20.0)
p = (0.2,0.2,3.0)
prob = ODEProblem(fitz,u0,tspan,p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 20.0)
u0: 2-element Vector{Float64}:
 -1.0
  1.0

Now we define the ProbInts callback. In this case, our method is the Euler method and thus it is order 1. For the noise scaling, we will try a few different values and see how it changes. For σ=0.2, we define the callback as:

cb = ProbIntsUncertainty(0.2,1)
DiscreteCallback{DiffEqCallbacks.var"#66#67", DiffEqCallbacks.ProbIntsCache{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}(DiffEqCallbacks.var"#66#67"(), DiffEqCallbacks.ProbIntsCache{Float64}(0.2, 1), SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 0])

This is akin to having an error of approximately 0.2 at each step. We now build and solve a EnsembleProblem for 100 trajectories:

ensemble_prob = EnsembleProblem(prob)
sim = solve(ensemble_prob,Euler(),trajectories=100,callback=cb,dt=1/10)
EnsembleSolution Solution of length 100 with uType:
ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Float64, Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.fitz), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Euler, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.fitz), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.EulerCache{Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats}

Now we can plot the resulting Monte Carlo solution:

plot(sim,idxs=(0,1),linealpha=0.4)