I’d like to do a Bayesian inversion of parameters to an ODE model. I would like to add some uncertainty on the input time series. Is there any way to do this? I’ve tried to add a parameter tau ~ Uniform(-60, 60)
and used it in the saveat = times .+ tau
part of the ODE specification, but this gets me in trouble.
Here’s a slightly big MWE based on the Lotka-Volterra example code. First I rewrote the example to use the new function(x) | (y = data)
syntax to see if I understood that correctly. (I had to add a if !ismissing(x)
block around stuff that relied on x for the calculation to allow me to sample from the prior).
My first error was the forwarddiff dual cache size issue I’ve run into a few times now, while the second was a dimension mismatch (might be the same?)
I had first posted this on the Turing slack but someone referred me to the DifferentialEquations side of things.
using Turing
using DifferentialEquations
using FillArrays
# Load StatsPlots for visualizations and diagnostics.
using StatsPlots
using LinearAlgebra
# Set a seed for reproducibility.
using Random
Random.seed!(14);
# Define Lotka-Volterra model.
function lotka_volterra(du, u, p, t)
# Model parameters.
α, β, γ, δ = p
# Current state.
x, y = u
# Evaluate differential equations.
du[1] = (α - β * y) * x # prey
du[2] = (δ * x - γ) * y # predator
return nothing
end
# Define initial-value problem.
u0 = [1.0, 1.0]
p = [1.5, 1.0, 3.0, 1.0]
tspan = (0.0, -20.0) # slightly larger than extrema(times) because of tau
prob = ODEProblem(lotka_volterra, u0, tspan, p)
# # Plot simulation.
plot(solve(prob, Tsit5()))
# in my use-case time goes backwards, wanted to ensure this isn't
# the cause of the problem
times = 0:(-0.1):(-10)
sol = solve(prob, Tsit5(); saveat=times)
odedata = Array(sol) + 0.8 * randn(size(Array(sol)))
# # Plot simulation and noisy observations.
# plot(sol; alpha=0.3)
# scatter!(sol.t, odedata'; color=[1 2], label="")
# used the second version, which only has predator data
@model function fitlv(data::AbstractVector, times, prob)
# Prior distributions.
σ ~ InverseGamma(2, 3)
α ~ truncated(Normal(1.5, 0.5); lower=0.5, upper=2.5)
β ~ truncated(Normal(1.2, 0.5); lower=0, upper=2)
γ ~ truncated(Normal(3.0, 0.5); lower=1, upper=4)
δ ~ truncated(Normal(1.0, 0.5); lower=0, upper=2)
# Simulate Lotka-Volterra model.
p = [α, β, γ, δ]
# make sure tspan is large enough
newprob = remake(prob; p = p, tspan = [0; minimum(times)])
predicted = solve(newprob, Tsit5();
# seems like I cannot adjust tspan here?
# my particular solution needs t = 0 to be included
saveat = [0; times],
# simulate full system but save only second state
save_idxs = 2)
# Observations
data ~ MvNormal(predicted.u, σ^2 * I)
return nothing
end
model = fitlv(odedata[2,:], times, prob)
# Sample 3 independent chains with forward-mode automatic differentiation (the default).
chain = sample(model, NUTS(), MCMCSerial(), 1000, 3; progress=false)
plot(chain)
plot(; legend=false)
posterior_samples = sample(chain[[:α, :β, :γ, :δ]], 300; replace=false)
for p in eachrow(Array(posterior_samples))
sol_p = solve(prob, Tsit5(); p=p, saveat=times)
plot!(sol_p; alpha=0.1, color="#BBBBBB")
end
# Plot simulation and noisy observations.
plot!(sol; color=[1 2], linewidth=1)
scatter!(sol.t, odedata'; color=[1 2])
# create the same model but using the new conditioning | data syntax
@model function fitlv2(times, prob)
# Prior distributions.
σ ~ InverseGamma(2, 3)
α ~ truncated(Normal(1.5, 0.5); lower=0.5, upper=2.5)
β ~ truncated(Normal(1.2, 0.5); lower=0, upper=2)
γ ~ truncated(Normal(3.0, 0.5); lower=1, upper=4)
δ ~ truncated(Normal(1.0, 0.5); lower=0, upper=2)
# Simulate Lotka-Volterra model.
p = [α, β, γ, δ]
# I think this is needed to be able to sample from the prior
if !ismissing(times) && !ismissing(prob)
# make sure tspan is large enough
newprob = remake(prob; p = p, tspan = [0; minimum(times)])
predicted = solve(newprob, Tsit5();
# seems like I cannot adjust tspan here?
# my particular solution needs t = 0 to be included
saveat = [0; times],
# simulate full system but save only second state
save_idxs = 2)
# Observations
data ~ MvNormal(predicted.u, σ^2 * I)
end
return nothing
end
model2_prior = fitlv2(missing, prob)
prior2 = sample(model2_prior, NUTS(), MCMCSerial(), 1000, 3; progress=false)
model2 = fitlv2(times, prob) | (data = odedata[2,:], )
chain2 = sample(model2, NUTS(), MCMCSerial(), 1000, 3; progress=false)
# finally, our new addition to the model:
# what if there is some uncertainty on times?
# (in this case this matters because there is some external forcing function that depends on t)
@model function fitlvtau(times, prob)
# Prior distributions.
σ ~ InverseGamma(2, 3)
α ~ truncated(Normal(1.5, 0.5); lower=0.5, upper=2.5)
β ~ truncated(Normal(1.2, 0.5); lower=0, upper=2)
γ ~ truncated(Normal(3.0, 0.5); lower=1, upper=4)
δ ~ truncated(Normal(1.0, 0.5); lower=0, upper=2)
τ ~ Uniform(0, 0.3) # tau is the uncertainty on the time of the sample
# I also tried the below but I'm confused here...
# τ ~ filldist(truncated(Normal(0, 0.3); lower = 0, upper = 1),
# length(data))
if !ismissing(times)
# set up our ages with uncertainty added
# T = [0; times .+ τ] # first attempt
# ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Ta
# T = times .+ τ # what if it's the [0;] part?
# ERROR: MethodError: no method matching zero(::Base.TwicePrecision{Float64})
# then I tried setting it up like this:
# for i, t in enumerate(times)
# T[i] ~ Normal(t, τ^2)
# end
# and then switched to MvN for speed...
T ~ MvNormal(times, τ * I)
# which fails with
# ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 12})
# Simulate Lotka-Volterra model.
p = [α, β, γ, δ]
# make sure tspan is large enough
newprob = remake(prob; p = p, tspan = [0; minimum(times)])
predicted = solve(newprob, Tsit5();
saveat = T,
save_idxs = 2)
# real example would now interpolate some external value
# to timesteps of sol.t and then do some additional calculations here
# Observations.
data ~ MvNormal(predicted.u, σ * I)
end
return nothing
end
modtau_prior = fitlvtau(missing, prob)
chaintau_prior_miss = sample(modtau_prior, NUTS(), MCMCSerial(), 1000, 3)
# this works
modtau_prior = fitlvtau(times, prob)
chaintau_prior = sample(modtau_prior, NUTS(), MCMCSerial(), 1000, 3)
# fails with the dual cache problem
# ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DynamicPPL.DynamicPPLTag, Float64}, Float64, 12})
modtau = fitlvtau(times, prob) | (data = odedata[2,:], )
chaintau = sample(modtau, NUTS(), MCMCSerial(), 1000, 3)
# fails with this error:
# ERROR: DimensionMismatch: inconsistent array dimensions
# mess-around area
@model function fitlvtau(times, prob)
# Prior distributions.
σ ~ InverseGamma(2, 3)
α ~ truncated(Normal(1.5, 0.5); lower=0.5, upper=2.5)
β ~ truncated(Normal(1.2, 0.5); lower=0, upper=2)
γ ~ truncated(Normal(3.0, 0.5); lower=1, upper=4)
δ ~ truncated(Normal(1.0, 0.5); lower=0, upper=2)
τ ~ Uniform(0, 0.3) # tau is the uncertainty on the time of the sample
# I also tried the below but I'm confused here...
# τ ~ filldist(truncated(Normal(0, 0.3); lower = 0, upper = 1),
# length(data))
if !ismissing(times)
# set up our ages with uncertainty added
# first I tried simply times .+ τ but this also didn't work
# then I tried setting it up like this:
# for i, t in enumerate(times)
# T[i] ~ Normal(t, τ^2)
# end
# and then switched to MvN for speed...
# T ~ MvNormal(times, τ * I)
# Simulate Lotka-Volterra model.
p = [α, β, γ, δ]
# make sure tspan is large enough
predicted = solve(prob, Tsit5();
saveat = T, # <- this is where the error points
save_idxs = 2)
# real example would now interpolate some external value
# to timesteps of sol.t and then do some additional calculations here
# Observations.
data ~ MvNormal(predicted.u, σ * I)
end
return nothing
end