AtomisticQoIs

AtomisticQoIs

Author: Joanna Zou <jjzou@mit.edu> Version: 0.1.0 Year: 2024 Notes: A Julia library for computing expectation-valued quantities of interest (QoIs) from molecular dynamics simulation.

@kwdef mutable struct AutocorrelationQoI <: QoI

Defines quantities for computing two-point invariant statistics, in particular: mean and variance of the autocorrelation of the function H over the Gibbs measure p with lag t, e. g. Ep[H(xt)H(x0)] and Varp[H(xt)H(x0)]

Arguments

• H :: Function : function to evaluate autocorrelation
• lag :: Real : autocorrelation lag
• s :: Function : score (drift) function

struct CircularDomain <: HittingDomain

Defines circular domain for N-D state space.

Arguments

• center :: Vector : coordinate of center point (N-D vector)
• radius :: Real : circle radius

Divergence

A struct of abstract type Divergence contains properties needed to compute a divergence metric between measures of the stochastic process.

struct GaussQuadrature <: QuadIntegrator

Contains parameters for 1D or 2D Gauss quadrature integration.

Arguments

• ξ :: Union{Vector{<:Real}, Vector{<:Vector}} : quadrature points
• w :: Union{Vector{<:Real}, Vector{<:Vector}} : quadrature weights
• d :: Real : dimension

GibbsIntegrator

A struct of abstract type GibbsIntegrator computes the expectation of a function h(x, θ) with respect to the invariant measure p(x, θ).

@kwdef mutable struct GibbsQoI <: QoI

Defines quantities for computing one-point invariant statistics, in particular: mean and variance of the observable function h over the Gibbs measure p, e. g. Ep[h(θ)] and Varp[h(θ)]

Arguments

• h :: Function : function to evaluate expectation, receiving inputs (x,θ)
• p :: Gibbs : Gibbs (invariant) distribution with θ == nothing (undefined parameters)
• ∇h :: Union{Function, Nothing} = nothing : gradient of h (if nothing, gradient is computed using ForwardDiff)

@kwdef mutable struct GreenKuboQoI <: QoI

Defines quantities for computing a Green Kubo integral, in particular: the integrated autocorrelation function up to time lag t, e. g. ∫₀ᵗ E_p[H(xs)H(x0)] ds

# Arguments

• H :: Function : function to evaluate autocorrelation
• lag :: Real : autocorrelation lag
• s :: Function : score (drift) function

struct HMC <: Sampler

Defines the struct with parameters of the Hamiltonian Monte Carlo algorithm for sampling.

Arguments

• L :: Int : time length of integrating Hamiltonian's equations
• ϵ :: Real : step size

HittingDomain

Abstract type defining domain in state space for computing the first hitting time

@kwdef mutable struct HittingTimeQoI <: QoI

Defines quantities for computing statistics of the first hitting time, in particular: mean and variance of the first hitting time τ over the transient path measure P up to time τ, e. g. EP[τ(θ)] and VarP[τ(θ)]

Arguments

• x0 :: Union{Real, Vector{<:Real}} : initial state of the stochastic process
• D :: HittingDomain : hitting domain
• s :: Function : score (drift) function of the stochastic process

struct ISMC <: ISIntegrator

Contains parameters for Importance Sampling (IS) integration using MC samples from a biasing distribution. This method is implemented when the biasing distribution g can be analytically sampled using rand().

Arguments

• g :: Distribution : biasing distribution
• n :: Int : number of samples

struct ISMCMC <: ISIntegrator

Contains parameters for Importance Sampling (IS) integration using MCMC samples from a biasing distribution. This method is implemented when the biasing distribution g cannot be analytically sampled.

Arguments

• g :: Distribution : biasing distribution
• n :: Int : number of samples
• sampler :: Sampler : type of sampler (see Sampler)
• ρ0 :: Distribution : prior distribution of the state

struct ISMixSamples <: ISIntegrator

Contains parameters for Importance Sampling (IS) integration using samples from a mixture biasing distribution. The mixture weights are computed based on a kernel distance metric of the parameter with respect to parameters of the mixture component distributions. This method is implemented with the user providing samples from each component mixture distribution.

Arguments

• g :: MixtureModel : mixture biasing distribution
• n :: Int : number of samples
• knl :: Kernel : kernel function to compute mixture weights
• xsamp :: Vector : Vector of sample sets from each component distribution
• normint :: GibbsIntegrator : Integrator for the approximating the normalizing constant of each component distribution

struct ISSamples <: ISIntegrator

Contains pre-defined Monte Carlo samples from the biasing distribution to evaluate in integration. This method is implemented with the user providing samples from the biasing distribution.

Arguments

• g :: Distribution : biasing distribution
• xsamp :: Vector : fixed set of samples
• normint :: Union{GibbsIntegrator, Nothing} : integrator for computing normalizing constant of biasing distribution (required for mixture models)

struct IntervalDomain <: HittingDomain

Defines bounds of the interval domain for 1D state space.

Arguments

• bounds :: Vector{<:Real} : lower and upper bounds [ll, ul]

struct MALA <: Sampler

Defines the struct with parameters of the Metropolis-Adjusted Langevin Algorithm (MALA) for sampling.

Arguments

• L :: Int : time length between evaluating accept/reject criterion
• ϵ :: Real : step size

struct MCMC <: MCIntegrator

Contains parameters for Monte Carlo (MC) integration using MCMC samples. This method is implemented when the distribution cannot be analytically sampled.

Arguments

• n :: Int : number of samples
• sampler :: Sampler : type of sampler (see Sampler)
• ρ0 :: Union{Distribution, Real, Vector{<:Real}} : initial state or prior distribution of the state

struct MCPaths <: PathIntegrator

Contains parameters for generating path realizations to compute Monte Carlo averages over path measures.

Arguments

• M :: Int : number of paths
• n :: Int : max number of samples per path
• sampler :: Sampler : type of sampler (see Sampler)
• ρ0 :: Union{Distribution, Real, Vector{<:Real}} : initial state or prior distribution of the state

struct MCSamples <: MCIntegrator

Contains pre-defined Monte Carlo samples to evaluate in integration. This method is implemented with the user providing samples from the distribution.

Arguments

• xsamp :: Vector : fixed set of samples

struct MonteCarlo <: MCIntegrator

Contains parameters for Monte Carlo (MC) integration using MC samples. This method may only be implemented when the distribution can be analytically sampled using rand().

Arguments

• n :: Int : number of samples

struct NUTS <: Sampler

Defines the struct with parameters of the No-U-Turn-Sampler algorithm for sampling.

Arguments

• ϵ :: Union{<:Real, Nothing} : step size

PathIntegrator

A struct of abstract type PathIntegrator computes the expectation of a functional h with respect to the path measure P.

QoI

Abstract type defining observable quantities of interest (QoI) from a stochastic system.

struct RectangularDomain <: HittingDomain

Defines rectangular domain for N-D state space.

Arguments

• bounds :: Vector : 2-dim. vector of uniform lower and upper bounds, or N-dim. vector of lower and upper bounds in each dimension
• dim :: Real : dimension of state space (N)

function EffSampleSize(q::GibbsQoI, xsamp::Vector)

Computes the effective sample size of a set of nMC samples from an MCMC chain.

Arguments

• q :: GibbsQoI : Gibbs QoI struct
• xsamp :: Vector : nMC-vector of x-samples

Outputs

• ESS :: Float64 : effective sample size (where 1 < ESS < nMC)

function EffSampleSize(hsamp::Vector)

Computes the effective sample size of a set of nMC samples from an MCMC chain.

Arguments

• hsamp :: Vector : nMC-vector of d-dimensional function evaluation h

Outputs

• ESS :: Float64 : effective sample size (where 1 < ESS < nMC)

function ISWeightDiagnostic(wsamp::Vector)

Returns a diagnostic value for importance sampling weights using an unnormalized biasing distribution. The importance sampling estimate is reliable when the diagnostic is less than 5.

Arguments

• wsamp :: Vector : nMC-vector of 1-dim. weights

Outputs

• diagnostic :: Float64 : diagnostic value (want < 5)

function ISWeightVariance(wsamp::Vector)

Returns the variance of importance sampling weights.

Arguments

• wsamp :: Vector : nMC-vector of weights

Outputs

• var(wsamp) :: Float64 : variance of weights

function MCSE(q::GibbsQoI, xsamp::Vector)

Estimates the Monte Carlo standard error (se = std(h(x))/√n) using an empirical estimate of standard deviation.

Arguments

• q :: GibbsQoI : Gibbs QoI struct
• xsamp :: Vector : vector of samples

Outputs

• se :: Float64 : Monte Carlo standard error

function MCSEbm(q::GibbsQoI, xsamp::Vector)

Estimates the Monte Carlo standard error (se = std(h(x))/√n) using the batch means (BM) method.

Arguments

• q :: GibbsQoI : Gibbs QoI struct
• xsamp :: Vector : vector of samples

Outputs

• se :: Float64 : Monte Carlo standard error

function MCSEbm(q::GibbsQoI, xsamp::Vector)

Estimates the Monte Carlo standard error (se = std(h(x))/√n) using the overlapping batch means (OBM) method.

Arguments

• q :: GibbsQoI : Gibbs QoI struct
• xsamp :: Vector : vector of samples

Outputs

• se :: Float64 : Monte Carlo standard error

function assign_param(qoi::GibbsQoI, θ::Vector)

Helper function modifying a GibbsQoI for a specific value of the parameter θ

function computegradqoi(θ::Union{Real, Vector{<:Real}}, qoi::GibbsQoI, integrator::GibbsIntegrator; gradh::Union{Function, Nothing}=nothing)

Computes the gradient of the QoI with respect to the parameters θ, e. g. ∇θ E_p[h(θ)].

Arguments

• θ :: Union{Real, Vector{<:Real}} : parameters to evaluate expectation
• qoi :: QoI : QoI object containing summand h and Gibbs measure p
• integrator :: GibbsIntegrator : struct specifying method of integration
• gradh :: Union{Function, Nothing} : gradient of qoi.h; if nothing, compute with ForwardDiff

Outputs

• res :: NamedTuple : contains output quantities (see compute_qoi for more info)

function compute_metrics(θ::Vector, h::Vector, w::Vector)

Computes importance sampling diagnostic metrics, including variance of IS weights ("wvar"), ESS of IS weights ("wESS"), and diagnostic measure ("wdiag").

Arguments

• θ :: Vector : vector of parameters
• h :: Vector : evaluations of summand function h(x) at samples
• w :: Vector : evaluations of importance sampling weights at samples

Outputs

• metrics :: Dict{String, Vector} : dictionary of importance sampling diagnostics

function expectationisstable(xsamp::Vector, ϕ::Function, f::Gibbs, g::Distribution; normint=nothing)

Helper function for computing stable importance sampling weights (using log formulation).

function gaussquad(nquad::Integer, gaussbasis::Function, args...; limits::Vector=[-1,1])

Computes 1D Gauss quadrature points with a change of domain to limits=[ll, ul].

Arguments

• nquad::Integer : number of quadrature points
• gaussbasis::Function : polynomial basis function from FastGaussQuadrature (gausslegendre, gaussjacobi, etc.)
• args... : additional arguments for the gaussbasis() function
• limits::Vector : defines lower and upper limits for rescaling

Outputs

• ξ :: Vector{<:Real} : quadrature points ranging [ll, ul]
• w :: Vector{<:Real} : quadrature weights

function gaussquad_2D(nquad::Integer, gaussbasis::Function, args...; limits::Vector=[-1,1])

Computes 2D (symmetric) Gauss quadrature points with a change of domain to limits=[ll, ul].

Arguments

• nquad::Integer : number of quadrature points
• gaussbasis::Function : polynomial basis function from FastGaussQuadrature (gausslegendre, gaussjacobi, etc.)
• args... : additional arguments for the gaussbasis() function
• limits::Vector : defines lower and upper limits for rescaling

Outputs

• ξ :: Vector{<:Real} : quadrature points ranging [ll, ul]
• w :: Vector{<:Real} : quadrature weights

function hasapproxnormconst(d::Distribution)

Returns 'true' if distribution d requires an approximation of its normalizing constant.

Arguments

• d :: Distribution : distribution to check

Outputs

• flag :: Bool : true or false

function hasupdf(d::Distribution)

Returns 'true' if distribution d has an unnormalized pdf function (updf(), logupdf()).

Arguments

• d :: Distribution : distribution to check

Outputs

• flag :: Bool : true or false

function hits_domain(x::Union{Real, Vector{<:Real}, D::HittingDomain)

Returns bool of whether the point x falls within the hitting domain.

Arguments

• x::Union{Real, Vector{<:Real} : current state in the simulation trajectory
• D::HittingDomain : hitting time domain

function sample(lπ::Function, gradlπ::Function, sampler::NUTS, n::Int64, x0::Real)

Draw samples from target distribution π(x) using the No-U-Turn Sampler (NUTS) for scalar-valued support (x::Real).

Arguments

• lπ :: Function : log likelihood of target π
• gradlπ :: Function : gradient of log likelihood of π
• sampler :: HMC : Sampler struct specifying sampling algorithm
• n :: Int64 : number of samples
• x0 :: Real : initial state

Outputs

• samples :: Vector : vector of samples from π

function sample(lπ::Function, gradlπ::Function, sampler::NUTS, n::Int64, x0::Vector{<:Real})

Draw samples from target distribution π(x) using the No-U-Turn Sampler (NUTS) from the AdvancedHMC package. Assumes multi-dimensional support (x::Vector).

Arguments

• lπ :: Function : log likelihood of target π
• gradlπ :: Function : gradient of log likelihood of π
• sampler :: HMC : Sampler struct specifying sampling algorithm
• n :: Int64 : number of samples
• x0 :: Vector{<:Real} : initial state

Outputs

• samples::Vector{Vector} : vector of samples from π

function sample(lπ::Function, gradlπ::Function, sampler::Union{HMC,MALA}, n::Int64, x0::Real)

Draw samples from target distribution π(x) for scalar-valued support (x::Real).

Arguments

• lπ :: Function : log likelihood of target π
• gradlπ :: Function : gradient of log likelihood of π
• sampler :: Union{HMC,MALA} : Sampler struct specifying sampling algorithm
• n :: Int64 : number of samples
• x0 :: Real : initial state

Outputs

• samples :: Vector : vector of samples from π