如何在 Julia 中调整贝叶斯 ODE 的超参数?

How to tune the hyperparameters of a Bayesian ODE fit in Julia?

我一直在尝试使用不同的 ODE 方程来复制 https://diffeqflux.sciml.ai/dev/examples/BayesianNODE_NUTS/,但我收到的结果没有不确定性量化,是不是因为我的初始值 u0 更高:

你能告诉我哪里出了问题吗?

   using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Plots, AdvancedHMC, MCMCChains
using JLD, StatsPlots

function Arps!(du,u,p,t)
    y= u[1]
    #x, y = u
   # Di,b,n,tau = p
    n,tau = p
    #du[1]=dx=-(x * Di * x^b)
    du[1]=dy=-(n *((t^n)/tau) * y/t)
end
tspan=(1.0,50.0)
tsteps = 1:1:50
u0 = [16382.9]
p=[0.48,15.92]
prob_trueode = ODEProblem(Arps!,u0,tspan,p)
ode_data = Array(solve(prob_trueode, Tsit5(), saveat = tsteps))
ode_data =ode_data[1,:]

dudt= FastChain(FastDense(1, 30, tanh),
                  FastDense(30, 1))

prob_neuralode = NeuralODE(dudt, tspan, Tsit5(), saveat = tsteps)

function predict_neuralode(p)
    Array(prob_neuralode(u0, p))
end

function loss_neuralode(p)
    pred = predict_neuralode(p)
    loss = sum(abs2, ode_data .- pred)
    return loss, pred
end

l(θ) = -sum(abs2, ode_data .- predict_neuralode(θ)) - sum(θ .* θ)


function dldθ(θ)
    x,lambda = Flux.Zygote.pullback(l,θ)
    grad = first(lambda(1))
    return x, grad
end

metric  = DiagEuclideanMetric(length(prob_neuralode.p))

h = Hamiltonian(metric, l, dldθ)


integrator = Leapfrog(find_good_stepsize(h, Float64.(prob_neuralode.p)))


prop = AdvancedHMC.NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)

adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.45, prop.integrator))

samples, stats = sample(h, prop, Float64.(prob_neuralode.p), 500, adaptor, 500; progress=true)
losses = map(x-> x[1],[loss_neuralode(samples[i]) for i in 1:length(samples)])

################### RETRODICTED PLOTS: TIME SERIES #################
pl = scatter(tsteps, ode_data, color = :red, label = "Data: Var1", xlabel = "t", title = "Spiral Neural ODE")
for k in 1:300
    resol = predict_neuralode(samples[100:end][rand(1:400)])
    plot!(tsteps,resol[1,:], alpha=0.04, color = :red, label = "")
end

idx = findmin(losses)[2]
prediction = predict_neuralode(samples[idx])
plot!(tsteps,prediction[1,:], color = :black, w = 2, label = "")

最可能的原因是后验样本的损失函数幅度过高,因此后验样本结果超出范围并且在您的图上不可见。

这可以通过 (a) 添加神经 ODE 输出的比例因子并确保损失函数不是从非常高的幅度开始或 (b) 增加神经网络中的层数来解决架构/改变激活函数。

通过在Neural ODE中添加比例因子,我得到了如下图所示的良好结果: