Multivariate Distributions & Parameter Inference

willsharpless · March 24, 2021, 9:33pm

Hi all,

I was trying to adapt the example from this tutorial on parameter inference
to a generalized Lotka-Volterra model with n equations when I ran into some issues. In order to write multivariate distributions into the Turing model (and avoid the much discussed (Turing type error: expected Float64, got ForwardDiff.Dual)) I followed the instructions outlined at the performance tips note, however, in this workflow this pushed the error to the subsequent “predicted = …” line when evaluating with the sampled values because my ODE function is not fit to handle the ForwardDiff.Dual{…} type.

How should I restructure this to get around the type errors?

Here is the code:

##
using DifferentialEquations
using Plots
using Random
using Turing
# Turing.turnprogress(false)
## Generate Model and Synthetic Data
# Random.seed!(14);

function gLV!(dx,x,p,t)
    x[x .< 1e-8] .= 0
    mu, A, np = p
    dx .= x.*(mu + A*x)

    # no cc catch
    if any(x .> 100.0)
        dx .= 0
    end
end

function mult_noise!(dx,x,p,t)
    mu, A, np = p
    dx .= x.*np
end

mu = [0.4; 0.8; 0.2].*0.1
A = [-1.5 0.5 -0.2; -0.7 -1.3 0.1; 0.8 -0.2 -0.9]
np = [0.01; 0.01; 0.01] # species noise percentages

x0 = [0.1; 0.1; 0.1].*0.1
tspan = (0.0, 144.0)

# Definite Model with Read-Out Noise
np.=0

prob1 = ODEProblem(gLV!, x0, tspan, (mu, A, np))
sol = solve(prob1, Tsit5(), saveat = 24)
syndata = Array(sol) + 0.003*randn(size(Array(sol)))
plot(sol, alpha=0.3)
scatter!(sol.t, syndata')
title!("Definite")


## MCMC Fittin
Turing.setadbackend(:forwarddiff)

@model function fitglv(data, prob1, ::Type{T} = Float64) where {T}
    σ ~ InverseGamma(2, 3)

    n = size(prob1.u0,1)
    mu = Vector{T}(undef, n)
    A = Matrix{T}(undef, n, n)

    for i = 1:n
        mu[i] ~ truncated(Normal(0.4, 0.2), 0, 1.5)
    end

    for i = 1:n
        for j = 1:n
            if i == j
                A[i,j] ~ truncated(Normal(-0.8,1),-3,0)
            else
                A[i,j] ~ truncated(Normal(-0.3,1),-3,3)
            end
        end
    end

    np .= 0
    p = (mu, A, np)
    prob = remake(prob1, p=p)
    predicted = solve(prob, Tsit5(), saveat=24)

    for i = 1:length(predicted)
        data[:,i] ~ MvNormal(predicted[i], σ)
    end
end

model = fitglv(syndata, prob1)
chain = mapreduce(c -> sample(model, NUTS(.65), 1000), chainscat, 1:3)

and the full error I am gettting:

TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing,Float64,7}
setindex!(::Array{Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.VarInfo{NamedTuple{(:σ, :mu, :A),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:σ,Tuple{}},Int64},Array{InverseGamma{Float64},1},Array{DynamicPPL.VarName{:σ,Tuple{}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:mu,Tuple{Tuple{Int64}}},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:mu,Tuple{Tuple{Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:A,Tuple{Tuple{Int64,Int64}}},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:A,Tuple{Tuple{Int64,Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64},DynamicPPL.Model{var"#103#104",(:data, :prob1, :T),(:T,),(),Tuple{Array{Float64,2},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Array{Float64,1},Array{Float64,2},Array{Float64,1}},ODEFunction{true,typeof(gLV!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,typeof(SciMLBase.DEFAULT_OBSERVED),Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},SciMLBase.StandardODEProblem},Type{Float64}},Tuple{Type{Float64}}},DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40},(),AdvancedHMC.DiagEuclideanMetric}},DynamicPPL.DefaultContext},Float64},Float64,7}, ::Int64) at array.jl:847
macro expansion at broadcast.jl:932 [inlined]
macro expansion at simdloop.jl:77 [inlined]
copyto! at broadcast.jl:931 [inlined]
copyto! at broadcast.jl:886 [inlined]
materialize! at broadcast.jl:848 [inlined]
materialize! at broadcast.jl:845 [inlined]
gLV!(::Array{Float64,1}, ::Array{Float64,1}, ::Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.VarInfo{NamedTuple{(:σ, :mu, :A),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:σ,Tuple{}},Int64},Array{InverseGamma{Float64},1},Array{DynamicPPL.VarName{:σ,Tuple{}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:mu,Tuple{Tuple{Int64}}},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:mu,Tuple{Tuple{Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:A,Tuple{Tuple{Int64,Int64}}},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:A,Tuple{Tuple{Int64,Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64},DynamicPPL.Model{var"#103#104",(:data, :prob1, :T),(:T,),(),Tuple{Array{Float64,2},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Array{Float64,1},Array{Float64,2},Array{Float64,1}},ODEFunction{true,typeof(gLV!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,typeof(SciMLBase.DEFAULT_OBSERVED),Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},SciMLBase.StandardODEProblem},Type{Float64}},Tuple{Type{Float64}}},DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40},(),AdvancedHMC.DiagEuclideanMetric}},DynamicPPL.DefaultContext},Float64},Float64,7},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.VarInfo{NamedTuple{(:σ, :mu, :A),Tuple{DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:σ,Tuple{}},Int64},Array{InverseGamma{Float64},1},Array{DynamicPPL.VarName{:σ,Tuple{}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:mu,Tuple{Tuple{Int64}}},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:mu,Tuple{Tuple{Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}},DynamicPPL.Metadata{Dict{DynamicPPL.VarName{:A,Tuple{Tuple{Int64,Int64}}},Int64},Array{Truncated{Normal{Float64},Continuous,Float64},1},Array{DynamicPPL.VarName{:A,Tuple{Tuple{Int64,Int64}}},1},Array{Float64,1},Array{Set{DynamicPPL.Selector},1}}}},Float64},DynamicPPL.Model{var"#103#104",(:data, :prob1, :T),(:T,),(),Tuple{Array{Float64,2},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Array{Float64,1},Array{Float64,2},Array{Float64,1}},ODEFunction{true,typeof(gLV!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,typeof(SciMLBase.DEFAULT_OBSERVED),Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},SciMLBase.StandardODEProblem},Type{Float64}},Tuple{Type{Float64}}},DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40},(),AdvancedHMC.DiagEuclideanMetric}},DynamicPPL.DefaultContext},Float64},Float64,7},2},Array{Float64,1}}, ::Float64) at PostBin_Analysis_v1.1.jl:33
ODEFunction at scimlfunctions.jl:334 [inlined]
initialize!(::Ordinar...

willsharpless · March 24, 2021, 9:52pm

Note, I tried using filldist() to see if it would work instead of the {T} fix, but my problem with this is that I need a different distribution for the diagonal of the matrix. My overwriting with CartesianIndex is not working:

Error:
MethodError: no method matching axes(::DistributionsAD.MatrixOfUnivariate{Continuous,Truncated{Normal{Float64},Continuous,Float64},FillArrays.Fill{Truncated{Normal{Float64},Continuous,Float64},2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}}}, ::Int64)