Main module.

Exports

IndexType = Union{Vector{<:Real}, Vector{<:Dates.TimeType}}

Permitted Julia types for Gaussian process time points. Real numbers are ingested directly, treated as time points.. Instances of Dates.TimeType are converted to numeric time points by using Dates.datetime2unix.

struct GPModel

A GPModel contains covariance kernel structures and parameters for modeling data.

Fields

• pf_state::Gen.ParticleFilterState: Internal particle set.
• config::GP.GPConfig=GP.GPConfig(): User-specific customization, refer to GP.GPConfig.
• ds::IndexType: Observed time points.
• y::Vector{<:Real}: Observed time series values.
• ds_transform::Transforms.LinearTransform: Transformation of time to direct space.
• y_transform::Transforms.LinearTransform: Transformation of observations to direct space.

Constructors

model = GPModel(
    ds::IndexType,
    y::Vector{<:Real};
    n_particles::Integer=8,
    config::GP.GPConfig=GP.GPConfig())

See also

To perform learning given the data, refer to

• AutoGP.fit_smc!
• AutoGP.fit_mcmc!
• AutoGP.fit_greedy!

add_data!(model::GPModel, ds::IndexType, y::Vector{<:Real})

Incorporate new observations (ds, y) into model.

covariance_kernels(model::GPModel)

Return Gaussian process covariance kernels in model.

effective_sample_size(model::GPModel)

Return effective sample size (ESS) of weighted particle collection.

function fit_greedy!(
        model::GPModel;
        max_depth::Integer=model.config.max_depth,
        verbose::Bool=false,
        check::Bool=false,
        callback_fn::Function=(; kwargs...) -> nothing)

Infer the structure and parameters of an appropriate Gaussian process covariance kernel for modeling the observed data. Inference is performed using greedy search, as described in Algorithm 2 of Kim and Teh, 2018. It is an error if max_depth is not a finite positive number.

A callback_fn can be provided to monitor the search progress at each stage. Its signature must contain a single varargs specifier, which will be populated with keys :model, :step, :elapsed at each step of the greedy search.

fit_mcmc!(
    model::GPModel;
    n_mcmc::Integer,
    n_hmc::Integer,
    biased::Bool=false,
    verbose::Bool=false,
    check::Bool=false,
    callback_fn::Function=(; kwargs...) -> nothing)

Perform n_mcmc steps of involutive MCMC on the structure, with n_hmc steps of Hamiltonian Monte Carlo sampling on the parameters per accepted involutive MCMC move.

A callback_fn can be provided to monitor the progress each MCMC step for which at least one particle (i.e, chain) accepted a transition. Its signature must contain a single varargs specifier, which will be populated with keys :model, :step, :elapsed.

fit_smc!(
    model::GPModel;
    schedule::Vector{<:Integer},
    n_mcmc::Int,
    n_hmc::Int,
    shuffle::Bool=true,
    biased::Bool=false,
    adaptive_resampling::Bool=true,
    adaptive_rejuvenation::Bool=false,
    hmc_config::Dict=Dict(),
    verbose::Bool=false,
    check::Bool=false,
    callback_fn::Function=(; kwargs...) -> nothing)

Infer the structure and parameters of an appropriate Gaussian process covariance kernel for modeling the observed data. Inference is performed using sequential Monte Carlo.

Arguments

• model::GPModel: Instance of the GPModel to use.
• schedule::Vector{<:Integer}: Schedule for incorporating data for SMC, refer to Schedule.
• n_mcmc::Int: Number of involutive MCMC rejuvenation steps.
• n_hmc::Int: Number of HMC steps per accepted involutive MCMC step.
• biased::Bool: Whether to bias the proposal to produce "short" structures.
• shuffle::Bool=true: Whether to shuffle indexes ds or incorporate data in the given order.
• adaptive_resampling::Bool=true: If true resamples based on ESS threshold, else at each step.
• adaptive_rejuvenation::Bool=false: If true rejuvenates only if resampled, else at each step.
• hmc_config::Dict: Configuration for HMC inference on numeric parameters. Allowable keys are:
  • n_exit::Integer=1: Number of successive rejections after which HMC loop is terminated.
  • L_param::Integer=10 Number of leapfrog steps for kernel parameters.
  • L_noise::Integer=10: Number of leapfrog steps for noise parameter.
  • eps_param::Float64=0.02: Step size for kernel parameters.
  • eps_noise::Float64=0.02: Step size for noise parameter.
• verbose::Bool=false: Report progress to stdout.
• check::Bool=false: Perform dynamic correctness checks during inference.
• config::GP.GPConfig=GP.GPConfig(): User-specific customization, refer to GP.GPConfig.
• callback_fn: A callback for monitoring inference, must be generated by AutoGP.Callbacks.make_smc_callback.

log_marginal_likelihood_estimate(model::GPModel)

Return estimate of marginal likelihood of data (in log space).

maybe_resample!(model::GPModel, ess_threshold::Real)

Resample the particle collection in model if ESS is below ess_threshold. Setting ess_threshold = AutoGP.num_particles(model) + 1 will ensure that resampling always takes place, since the ESS is upper bounded by the number of particles.

mcmc_parameters!(model::GPModel, n_hmc::Integer; verbose::Bool=false, check::Bool=false)

Perform n_hmc steps of Hamiltonian Monte Carlo sampling on the parameters.

mcmc_structure!(model::GPModel, n_mcmc::Integer, n_hmc::Integer;
    biased::Bool=false, verbose::Bool=false, check::Bool=false)

Perform n_mcmc steps of involutive MCMC on the structure, with n_hmc steps of Hamiltonian Monte Carlo sampling on the parameters per accepted involutive MCMC move.

num_particles(model::GPModel)

Return the number of particles.

observation_noise_variances(model::GPModel)

Return list of observation noise variances for each particle in model.

particle_weights(model::GPModel)

Return vector of normalized particle weights.

predictions = predict(
    model::GPModel,
    ds::IndexType;
    quantiles::Vector{Float64}=Float64[],
    noise_pred::Union{Nothing,Float64}=nothing)

Return predictions for new index point ds, and optionally quantiles corresponding to the provided quantiles (numbers between 0 and 1, inclusive). By default, the noise_pred of the new data is equal to the inferred noise of the observed data within each particle in model; using noise_pred=0. returns the posterior distribution over the noiseless function values.

The returned DataFrames.DataFrame has columns ["ds", "particle", "weight", "y_mean"], as well as any additional columns for the requested quantiles.

Example

julia> ds = [Dates.Date(2020,1,1), Dates.Date(2020,1,2)] # Dates to query
julia> GPModel.predict(model, ds; quantiles=[.025, 0.975])
16×6 DataFrame
 Row │ ds          particle  weight       y_0.025    y_0.975    y_mean
     │ Date        Int64     Float64      Float64    Float64    Float64
─────┼────────────────────────────────────────────────────────────────────
   1 │ 2020-01-01         1  4.97761e-22  -13510.0   14070.6      280.299
   2 │ 2020-01-02         1  4.97761e-22  -13511.0   14071.6      280.299
   3 │ 2020-01-01         2  0.279887       4504.73   8211.43    6358.08
   4 │ 2020-01-02         2  0.279887       4448.06   8154.3     6301.18
   5 │ 2020-01-01         3  0.0748059    -43638.6   65083.0    10722.2
   6 │ 2020-01-02         3  0.0748059    -43662.0   65074.7    10706.4
   7 │ 2020-01-01         4  0.60809      -17582.2   30762.4     6590.06
   8 │ 2020-01-02         4  0.60809      -17588.0   30771.5     6591.78

dist = predict_mvn(model::GPModel, ds::IndexType; noise_pred::Union{Nothing,Float64}=nothing)

Return an instance of Distributions.MixtureModel representing the overall posterior predictive distribution for data at index points ds. By default, the noise_pred of the new data is equal to the inferred noise of the observed data within each particle in model; using noise_pred=0. returns the posterior distribution over the noiseless function values.

The returned dist has precisely num_particles(model) components, each of type Distributions.MvNormal, with weights particle_weights(model). These objects can be retrieved using Distributions.components and Distributions.probs, respectively.

function predict_proba(model::GPModel, ds::IndexType, y::Vector{<:Real})

Compute predictive probability of data y at time points ds under model.

Example

julia> ds = [Dates.Date(2020,1,1), Dates.Date(2020,1,2)] # Dates to query
julia> y = [0.1, .0.5] # values to query
julia> GPModel.predict(model, ds, y)
7×3 DataFrame
 Row │ particle  weight       logp
     │ Int64     Float64      Float64
─────┼─────────────────────────────────
   1 │        1  0.0287155    -64.2388
   2 │        2  0.0437349    -59.7672
   3 │        3  0.576247     -62.6499
   4 │        4  0.00164846   -59.5311
   5 │        5  0.215255     -61.066
   6 │        6  0.134198     -64.4041
   7 │        7  0.000201078  -68.462

seed!(seed)

Set the random seed of the global random number generator.