Conditional Branching of Parameter in Turing.jl

Specific Domains Modelling & Simulations
, ,
3

Thank you so much for your advice.
Your activities have helped me a lot.

I’ve tried many things, but it didn’t work.

For example, in the following code, I removed γ from the parameter and wrote a conditional branch in the model, but it doesn’t seem to be working properly.

I think I’ve made a fundamental mistake.
I would appreciate it if you could show me a concrete method.

using Turing, Distributions, DifferentialEquations
using MCMCChains, Plots, StatsPlots
using Random
Random.seed!(12);

function lotka_volterra(du,u,p,t)
    x, y = u
    α, β, δ = p
    γ = 3.0
    du[1] = dx = (α - β*y)x
    du[2] = dy = (δ*x - γ)y
end
p = [1.5, 1.1, 1.0]
u0 = [1.0,1.0]
prob = ODEProblem(lotka_volterra,u0,(0.0,10.0),p)
sol = solve(prob,Tsit5())
plot(sol)

odedata1 = Array(solve(prob,Tsit5(),saveat=0.1))
odedata2 = odedata1 .+ rand()
odedata3 = odedata1 .+ rand()
odedata = zeros(Float64, 2, 101, 3)
odedata[:,:,1] = odedata1
odedata[:,:,2] = odedata2
odedata[:,:,3] = odedata3

Turing.setadbackend(:forwarddiff)

@model function fitlv(data)
    σ ~ 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)

    var1 = [80, 40, 70]
    γ = zeros(Float64, length(var1))

    for l in 1:length(var1)
        if var1[l] > 60
            γ[l] = 5.0
        else
            γ[l] = 2.0
        end
    end

    p = [α,β,δ]
    prob = ODEProblem(lotka_volterra,u0,(0.0,10.0),p)
    predicted = solve(prob,Tsit5(),saveat=0.1)

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

model = fitlv(odedata)
chain = sample(model, NUTS(.65),1000)
plot(chain)