API Reference

MH sample

Fields:

• theta: parameters
• x_t: state at current time step $t=0$
• x_init: initial state of the training trajectory, i.e., state at $t=-T$; this is used to adapt the proposal distribution in the staged sampler

ODE_MH(u_t::Function, t_m::AbstractVector{<:AbstractFloat}, y::AbstractMatrix{<:AbstractFloat}, t_span::Tuple{Float64,Float64}, K::Int, K_b::Int, k_d::Int, f_theta!::Function, g_theta!::Function, log_pdf_w_theta::Function, log_pdf_theta::Function, theta_0::AbstractVector{<:AbstractFloat}, log_pdf_x_init::Function, x_init_0::AbstractVector{<:AbstractFloat}, propose_z::Function, log_proposal_ratio_z::Function; ODE_solver=RK4(), ODE_solver_opts=(dt=0.1, adaptive=false), print_progress=true)

Run a Metropolis-Hastings (MH) sampler with ODE-integrated latent trajectory to obtain samples $\{\theta, x(t=0)\}^{[1:K]}$ from the joint parameter and state posterior distribution $p(\theta, x(t=0) \mid \mathbb{D}=\{u(t), y_{1:M}\})$.

Arguments

• u_t: input trajectory; function of time $t$
• t_m: time points of the output measurements, must lie within t_span
• y: output measurements
• t_span: timespan of the training trajectory
• K: number of models/scenarios to be sampled
• K_b: length of the burn in period
• k_d: number of models/scenarios to be skipped to decrease correlation (thinning)
• f_theta!: dynamics function parametrized by theta (mutating); has inputs $(\dot{x}, \theta, x, u, t)$
• g_theta!: measurement function parametrized by theta (mutating); has inputs $(g, \theta, x, u, t)$
• log_pdf_w_theta: function that returns the logarithm of the probability density function of the measurement noise parametrized by theta; has inputs $(\theta, w)$
• log_pdf_theta: function that returns the logarithm of the probability density function of theta (prior); has input $(\theta)$
• theta_0: $\theta$ used to initialize the MH sampler
• log_pdf_x_init: function that returns the logarithm of probability density function of $x(t=-T)$ (prior); has input $x(t=-T)$
• x_init_0: $x(t=-T)$ used to initialize the MH sampler
• propose_z: function that proposes new parameters $z' = [\theta'; x'(t=-T)]$ (proposal distribution); has input ($z = [\theta; x(t=-T)$])
• log_proposal_ratio_z: function that computes the logarithm of the ratio of proposal densities for $z$; has input arguments $(z, z')$
• ODE_solver: ODE solver algorithm to use (we use RK4 as default as this is also used in the optimal control formulation)
• ODE_solver_opts: ODE solver setting
• print_progress: if true, the progress is printed (default: true)

Returns

• MH_samples: MH samples
• acceptance_ratio: acceptance ratio of the MH sampler
• runtime: runtime of the sampling process

staged_ODE_MH(u_t::Function, t_m::AbstractVector{<:AbstractFloat}, y::AbstractMatrix{<:AbstractFloat}, t_span::Tuple{Float64,Float64}, K::Int, K_b::Int, k_d::Int, f_theta!::Function, g_theta!::Function, log_pdf_w_theta::Function, log_pdf_theta::Function, theta_0::AbstractVector{<:AbstractFloat}, log_pdf_x_init::Function, x_init_0::AbstractVector{<:AbstractFloat}, proposal_z_cov_0::AbstractMatrix{<:AbstractFloat}, M_chunk::Int, K_stage::Int, alpha::Union{AbstractFloat,AbstractVector{<:AbstractFloat}}; regularizer::Float64=1e-8, ODE_solver=RK4(), ODE_solver_opts=(dt=0.1, adaptive=false), print_progress=true)

Run a Metropolis-Hastings (MH) sampler with ODE-integrated latent trajectory with incremental data and adaptive proposal to obtain samples $\{\theta, x(t=0)\}^{[1:K]}$ from the joint parameter and state posterior distribution $p(\theta, x(t=0) \mid \mathbb{D}=\{u(t), y_{1:M}\})$. The number of data points used in the likelihood computation is gradually increased by a fixed chunk size. At each stage, the MH sampler is run on the current data subset, and the proposal distribution is adapted based on the empirical covariance of the collected samples.

Arguments

• u_t: input trajectory; function of time $t$
• t_m: time points of the output measurements, must lie within t_span
• y: output measurements
• t_span: timespan of the training trajectory
• K: number of models/scenarios to be sampled
• K_b: length of the burn in period
• k_d: number of models/scenarios to be skipped to decrease correlation (thinning)
• f_theta!: dynamics function parametrized by theta (mutating); has inputs $(\dot{x}, \theta, x, u, t)$
• g_theta!: measurement function parametrized by theta (mutating); has inputs $(g, \theta, x, u, t)$
• log_pdf_w_theta: function that returns the logarithm of the probability density function of the measurement noise parametrized by theta; has inputs $(\theta, w)$
• log_pdf_theta: function that returns the logarithm of the probability density function of theta (prior); has input $(\theta)$
• theta_0: $\theta$ used to initialize the MH sampler
• log_pdf_x_init: function that returns the logarithm of probability density function of $x(t=-T)$ (prior); has input $x(t=-T)$
• x_init_0: $x(t=-T)$ used to initialize the MH sampler
• proposal_z_cov_0: initial covariance of the multivariate normal proposal distribution for $z = [\theta; x(t=-T)]$ (e.g., covariance of the prior)
• M_chunk: number of data points added at each stage
• K_stage: number of samples per stage
• alpha: proposal scaling factor (scalar or vector; if a vector its i‐th element is used at stage i)
• regularizer: small constant added to the diagonal of the proposal covariance matrix to ensure positive definiteness (default: 1e-8)
• ODE_solver: ODE solver algorithm to use (we use RK4 as default as this is also used in the optimal control formulation)
• ODE_solver_opts: ODE solver setting
• print_progress: if true, the progress is printed (default: true)

Returns

• MH_samples: final samples from full-data posterior
• acceptance_ratio: vector containing the acceptance ratio of each stage
• runtime: total runtime of the staged sampling process

compute_autocorrelation(MH_samples::Vector{MH_sample}; max_lag::Int=100)

Compute the autocorrelation function (ACF) of the MH samples.

Arguments

• MH_samples: MH samples
• max_lag: maximum lag at which to calculate the ACF

Returns

• autocorrelation: matrix containing the ACF for each variable

compute_ess(MH_samples::Vector{MH_sample}; max_lag::Int=100)

Compute the effective sample size (ESS) for each parameter and initial state.

Arguments

• MH_samples: MH samples
• max_lag: maximum lag for autocorrelation estimation

Returns

• ess: vector of ESS estimates for all variables

compute_gelman_rubin(MH_chains::Vector{Vector{MH_sample}})

Compute the Gelman–Rubin statistic R̂ for each parameter and initial state from a vector of PMCMC chains. The Gelman–Rubin statistic R̂ quantifies convergence by comparing within-chain to between-chain variance. R̂ close to 1 (typically R̂ < 1.05) indicates good convergence across chains.

Arguments

• MH_chains: vector of chains, where each chain is a vector of MH samples

Returns

• R_hat: vector of R̂ values, one for each variable

solve_MH_OCP(MH_samples::Vector{MH_sample}, n_u::Int, f_theta::Function, g_theta::Function, H::Float64, N::Int, c::Function, c_f::Function, h_scenario::Function, h_u::Function; U_init=nothing, MH_samples_pre_solve=nothing, K_warmup=0, solver_opts=nothing, rk4_step_size=0.1, print_progress=true)

Solve the continous time scenario optimal control problem of the following form:

$\min_{u(\cdot)} J_{sc}(u(\cdot)) := \frac{1}{K} \sum_{k=1}^{K} [ c_f(x^{[k]}(H)) + \int_0^H c(u(t), x^{[k]}(t), t) dt],$

subject to:

\[\begin{aligned} \forall t \in [0,H], \forall k \in \mathbb{N}_{\leq K}, \\ x^{[k]}(t) = \Phi(t; \theta^{[k]}, x^{[k]}(0), u(\cdot)), \\ h_{scenario}(u(t), x^{[k]}(t), t) \leq 0. \end{aligned}\]

Here, $\Phi(t; \theta^{[k]}, x^{[k]}(0), u(\cdot))$ denotes the solution at time $t$ of the ODE $\dot{x}(t) = f_{\theta}(x(t), u(t))$ with initial condition $x(0) = x^{[k]}(0)$ and parameter $\theta$ under the input trajectory $u(\cdot)$.

Arguments

• MH_samples: MH samples
• n_u: number of control inputs
• f_theta: dynamics function parametrized by theta (non mutating); has inputs $(\theta, x, u, t)$
• g_theta: measurement function parametrized by theta (non mutating); has inputs $(\theta, x, u, t)$
• H: horizon of the OCP (as time)
• N: number of discretization steps in the horizon
• c: function with input arguments $(u(t), x^{[k]}(t), t)$ that returns the running cost of scenario k at time t
• c_f: function with input argument $x^{[k]}(H)$ that returns the terminal cost of scenario k
• h_scenario: function with input arguments $(u(t), x^{[k]}(t), t)$ that returns the constraint vector for scenario k at time t; a feasible solution must satisfy $h_{\mathrm{scenario}} \leq 0$ for all scenarios and all discretization points
• h_u: function with input arguments $(u(t), t)$ that returns the constraint vector for the control inputs; a feasible solution satisfy $h_u \leq 0$ at all discretization points
• U_init: initial guess for the input trajectory - either a n_u x N array, a function with input argument t, or nothing (default: nothing)
• MH_samples_pre_solve: if provided, an initial guess for the input trajectory is obtained by solving an OCP with the samples in MH_samples_pre_solve only
• K_warmup: if K_warmup > 0 and MH_samples_pre_solve is provided, an initial guess for the the input trajectory is obtained in a two stage process: first, an OCP with only K_warmup samples from MH_samples_pre_solve is solved and then an OCP with all samples in MH_samples_pre_solve
• solver_opts: SolverOptions struct containing options of the solver
• rk4_step_size: step size of the RK4 integrator used to simulate the system dynamics
• print_progress: if set to true, the progress is printed

Returns

• U_opt: piecewise constant optimal input trajectory, array of dimension n_u x N
• X_opt: states at the discretization points for all scenarios, array of dimension n_x x N x K
• t_grid: time grid of the discretization points, array of dimension N + 1
• J_sc_opt: optimal cost $J_{sc}$
• solve_successful: true if the optimization was successful, false otherwise
• iterations: number of iterations of the solver
• runtime: runtime of the optimization