How to infer diffusion parameters based on its discrete-time observations?
The most convenient solution is to use the package DiffusionMCMC.jl. Below, we explain how to do it without this package, using solely DiffusionDefinition.jl, ObservationSchemes.jl and GuidedProposals.jl.
Define some transition kernel for the parameter updates. We will use a simple random walker
# define a simple transition kernel
function customkernel(θ, s::Symbol, scale=0.1)
θ° = deepcopy(θ)
θ°[s] += 2.0*scale*(rand()-0.5)
θ°
end
The main routine
# Perform inference on a single parameter for the data in the `recording`, using
# Guided Proposals with the auxiliary law `AuxLaw`.
function simple_inference(AuxLaw, recording, dt, θ; ρ=0.5, num_steps=10^4, ϵ=0.3)
# -------------------------------------------------------------------------#
# Initializations #
# -------------------------------------------------------------------------#
# time-grids for the forward-simulation of trajectories #
tts = OBS.setup_time_grids(recording, dt) #
# laws of guided proposals #
PP = build_guid_prop(AuxLaw, recording, tts) #
# laws of guided proposals for parameter proposals #
PP° = deepcopy(PP) #
#
# starting point #
# NOTE `rand` for `KnownStartingPt` simply returns the starting position #
y1 = rand(recording.x0_prior) #
# initialize the `accepted` trajectory #
XX, WW, Wnr = rand(PP, y1) #
# initialize the containers for the `proposal` trajectory #
XX°, WW° = trajectory(PP) #
#
ll = loglikhd(PP, XX) #
paths = [] #
imp_num_accpt = 0 #
param_num_accpt = 0 #
pname = collect(keys(θ))[1] #
θθ = Float64[θ[pname],] #
# -------------------------------------------------------------------------#
# MCMC
for i in 1:num_steps
# impute a path
_, ll° = rand!(PP, XX°, WW°, WW, ρ, Val(:ll), y1; Wnr=Wnr)
# Metropolis–Hastings accept/reject step
if rand() < exp(ll°-ll)
XX, WW, XX°, WW° = XX°, WW°, XX, WW
ll = ll°
imp_num_accpt += 1
end
# update parameter s
θ° = customkernel(θ, pname, ϵ)
DD.set_parameters!(PP°, θ°)
recompute_guiding_term!(PP°)
_, ll° = GP.solve_and_ll!(XX°, WW, PP°, y1)
if rand() < exp(ll°-ll) # uniform updates have no contribution to ll
XX, PP, θ, XX°, PP°, θ° = XX°, PP°, θ°, XX, PP, θ
ll = ll°
param_num_accpt += 1
end
append!(θθ, [θ[pname]])
# progress message
if i % 100 == 0
println(
"$i. ll=$ll, $pname=$(θ[pname]), imp accpt rate: ",
"$(imp_num_accpt/100), updt accpt rate: $(param_num_accpt/100)"
)
imp_num_accpt = param_num_accpt = 0
end
# save intermediate path for plotting
i % 400 == 0 && append!(paths, [deepcopy(XX)])
end
paths, θθ
end
Example
We can now infer a single parameter of FitzHugh–Nagumo model. Define which parameters will not change:
DD.const_parameter_names(::Type{<:FitzHughNagumo}) = (:ϵ, :γ, :β, :σ)
DD.const_parameter_names(::Type{<:FitzHughNagumoAux}) = (:ϵ, :γ, :β, :σ, :t0, :T, :vT, :xT)
and then, run the inference
paths, θθ = simple_inference(
FitzHughNagumoAux, recording, 0.001, Dict(:s=>0.0); ρ=0.96, num_steps=10^4, ϵ=0.3
)
It takes about 25sec on my laptop...
The chain (truth used to generate the data is $0.8$)
And the paths: