Adding uncertainty on timestep of Bayesian ODE solve through `saveat`

Specific Domains Probabilistic programming
, ,
2

If i understand correctly, you are basically trying to fit a model that incorporates some type of measurement error in the time measurements? If that’s the case and you merely are trying to fit that parameter this I was able to sample from:

using Turing
using DifferentialEquations
using Distributions
using Random
using LinearAlgebra

# Set a seed for reproducibility
Random.seed!(14)

# Define Lotka-Volterra model
function lotka_volterra(du, u, p, t)
   α, β, γ, δ = p
   x, y = u
   du[1] = (α - β * y) * x # prey
   du[2] = (δ * x - γ) * y # predator
end

# Define initial-value problem
u0 = [1.0, 1.0]
p_true = [1.5, 1.0, 3.0, 1.0]
tspan = (0.0, 10.0)
prob = ODEProblem(lotka_volterra, u0, tspan, p_true)

# Generate synthetic data
times = 0.0:0.1:10.0
sol = solve(prob, Tsit5(); saveat=times)
odedata = Array(sol) + 0.1 * randn(size(Array(sol)))

@model function fitlv_time_uncertainty(times, prob, data)
    # Prior distributions
    σ ~ InverseGamma(2, 3)
    α ~ truncated(Normal(1.5, 0.5), 0.5, 2.5)
    β ~ truncated(Normal(1.2, 0.5), 0, 2)
    γ ~ truncated(Normal(3.0, 0.5), 1, 4)
    δ ~ truncated(Normal(1.0, 0.5), 0, 2)
    τ ~ Uniform(0, 0.3)

    # Time uncertainty
    T = times .+ rand(Normal(0, τ), length(times))

    # Simulate Lotka-Volterra model
    p = [α, β, γ, δ]
    predicted = solve(remake(prob, p=p, tspan=(0.0, maximum(T))), Tsit5(); saveat=T)

   # Observations
   for i in 1:length(times)
        data[:, i] ~ MvNormal(predicted(T[i]), σ^2 * I)
   end
end

# Fit the model
model = fitlv_time_uncertainty(times, prob, odedata)
chain = sample(model, NUTS(), 1000)

I essentially just assumed you were sampling time descrepancies from a 0 mean normal with some standard deviation parameter tau. Hope this is useful.

1 Like