API

A Spl object.

Fields

• H: an RObject generated by splines::bs()
• β: the coefficients for the B-spline.

aug(λ::AbstractFloat)

Augment λ with 8 different functions.

build_model(x::AbstractVector{T}; <keyword arguments>)

Construct design matrix and other internal variables for smoothing spline.

Arguments

• all_knots = false: whether to use all knots. If false, use the same rule as in R's smooth.spline.
• prop_nknots = 1.0: a proportion for using fewer knots. Suppose the number of knots is nknots, then the final number of knots is prop_nknots * nknots. Currently, it is only effective when all_knots = false.
• ε = 6.06e-6: a small number added to the diagonal of matrix Ω to ensure it is positive definite.

Returns

• B-spline design matrix B at x for cubic splines
• L: cholesky decomposition of Ω = LL'
• J: number of basis functions, which does not change for cubic splines, so it is only intended for smoothing splines

the above four are shared with the method for cubic splines, but for smoothing splines, it also returns

• mx, rx, idx, idx0: only for smoothing splines

build_model(x::AbstractVector{T}, J::Int; <keyword arguments>)

Construct design matrix and other internal variables for cubic spline with J basis functions.

Returns

• B: B-spline design matrix B at x for cubic splines
• rB: raw RObject of B

check_CI(; <keyword arguments>)

Conduct repeated experiments to check the overlap of confidence bands (default, check_acc = false) or accuracy of fitting curves (check_acc = true) between MLP generator and OPT solution.

Arguments

• n = 100: sample size
• σ = 0.1: noise level
• f::Function = exp: the truth curve
• seed = 1234: random seed for the simulated data
• check_acc = false: check overlap of confidence bands (default: false) or accuracy of fitting curves (true)
• nepoch0 = 5: number of epoch in the first step to fit the curve
• nepoch = 50: number of epoch in the second step to obtain the confidence band
• niter_per_epoch = 100: number of iterations in each epoch
• η0 = 1e-4: learning rate in step 1
• η = 1e-4: learning rate in step 2 (NOTE: lr did not make much difference, unify these two)
• K0 = 32: Monte Carlo size for averaging λ in step 2
• K = 32: Monte Carlo size for averaging λ in step 1 and for averaging y in step 2. (NOTE: unify these two Monte Carlo size)
• nB = 2000: number of bootstrap replications
• nrep = 5: number of repeated experiments
• fig = true: whether to plot
• figfolder = ~: folder for saving the figures if fig = true
• λs = exp.(range(-8, -2, length = 10)): region of continuous λ
• nhidden = 1000: number of hidden layers
• depth = 2: depth of MLP
• demo = false: whether to save internal results for demo purpose
• model_file = nothing: if not nothing, load the model from the file.
• gpu_id = 0: specify the id of GPU, -1 for CPU.
• prop_nknots = 0.2: proportion of number of knots in B-spline basis.
• backend = "flux": train MLP generator with Flux or PyTorch

ci_mono_ss_mlp(x::AbstractVector{T}, y::AbstractVector{T}, λs::AbstractVector{T}; )

Fit data x, y at each λs with confidence bands.

Arguments

• prop_nknots = 0.2: proportion of number of knots
• backend = "flux": flux or pytorch
• model_file: path for saving trained model
• nepoch0 = 3: number of epoch in training step 1
• nepoch = 3: number of epoch in training step 2
• niter_per_epoch = 100: number of iterations in each epoch
• M = 10: Monte Carlo size
• nhidden = 100: number of hidden units
• disable_progressbar = false: set true if generating documentation
• device = :cpu: train using :cpu or :gpu
• sort_in_nn = true: (only for backend = "flux") whether put sort in MLP
• eval_in_batch = false: (only for backend = "flux") Currently, Flux does not support sort in batch mode. A workaround with customized Zygote.batch_sort needs further verifications.

conf_band_width(CIs::AbstractMatrix)

Calculate width of confidence bands.

coverage_prob(CIs::AbstractMatrix, y0::AbstractVector)

Calculate coverage probability given n x 2 CI matrix CIs and true vector y0 of size n.

cv_mono_ss(x::AbstractVector{T}, y::AbstractVector{T}, λs::AbstractVector{T})

Cross-validation for monotone fitting with smoothing spline on y ~ x among parameters λs.

div_into_folds(N::Int; K = 10, seed = 1234)

Equally divide 1:N into K folds with random seed seed. If seed is negative, it is a non-random division, where the i-th fold would be the i-th equidistant range.

eval_penalty(model::Spl{T}, x::AbstractVector{T})

Evaluate the penalty matrix by R's fda::eval.penalty. To make sure the corresponding design matrix contructed by fda::eval.basis is the same as model.H, it asserts the norm difference should be smaller than sqrt(eps()).

fit(X, y, paras, method)

paras is either the number of basis functions, or the sequence of interior knots. Return a Spl object.

n = 100
x = rand(n) * 2 .- 1
y = x .^3 + randn(n) * 0.01
res = fit(x, y, 10, "monotone")

gen_data(n, σ, f::Union{Function, String}; xmin = -1, xmax = 1, k = 10)

Generate n data points (xi, yi) from curve f with noise level σ, i.e., yi = f(xi) + N(0, σ^2).

It returns four vectors, x, y, x0, y0, where

• x, y: pair points of length n.
• x0, y0: true curve without noise, represented by k*n points.

jaccard_index(a::AbstractVector, b::AbstractVector)

Calculate Jaccard Index for two confidence intervals a and b

jaccard_index(a::AbstractMatrix, b::AbstractMatrix)

Calculate Jaccard Index for two confidence intervals a[i, :] and b[i, :]

load_model(n::Int, J::Int, nhidden::Int, model_file::String; dim_lam = 8, gpu_id = 3)

Load trained model from model_file.

mono_cs(x::AbstractVector, y::AbstractVector, J::Int = 4; increasing::Bool = true)

Monotone splines with cubic splines.

mono_ss(x::AbstractVector, y::AbstractVector, λ = 1.0; prop_nknots = 1.0)

Monotone splines with smoothing splines, return a MonotoneSS object.

mono_ss(B::AbstractMatrix, y::AbstractVector, L::AbstractMatrix, J::Int, λ::AbstractFloat)

Monotone Fitting with Smoothing Splines given design matrix B and cholesky-decomposed matrix L.

Returns

• βhat: estimated coefficient
• yhat: fitted values
• (optional) B and L

mono_ss_mlp(x::AbstractVector, y::AbstractVector; λl, λu)

Fit monotone smoothing spline by training a MLP generator.

Arguments

• prop_nknots = 0.2: proportion of number of knots
• backend = flux: use flux or pytorch
• device = :cpu: use :cpu or :gpu
• nhidden = 100: number of hidden units
• disable_progressbar = false: disable progressbar (useful in Documenter.jl)

py_train_G_lambda(y::AbstractVector, B::AbstractMatrix, L::AbstractMatrix; <keyword arguments>)

Wrapper for training MLP generator using PyTorch.

Arguments

• η0, η: learning rate
• K0, K: Monte Carlo size
• nepoch0, nepoch: number of epoch
• nhidden, depth: size of MLP
• λl, λu: range of λ
• use_torchsort = false: torch.sort (default: false) or torchsort.soft_sort (true)
• sort_reg_strength = 0.1: tuning parameter when use_torchsort = true.
• model_file: path for saving trained model
• gpu_id = 0: use specified GPU
• niter_per_epoch = 100: number of iterations in each epoch
• disable_tqdm = false: set true when generating documentation
• λs_opt_train = nothing: if not nothing and instead a vector of λ (dim: N_t) used in the classical ML training step, then return the beta loss of GpBS procedure
• λs_opt_val = nothing: if not nothing and instead a vector of λ (dim: N_v) used in the classical ML validation step (integrative predictive loss), then return the beta loss of GpBS procedure
• βs_opt_train = nothing: if not nothing and instead a matrix of β (dim: N_t x J) used in the classical ML training step, then return the beta loss of GpBS procedure
• βs_opt_val = nothing: if not nothing and instead a matrix of β (dim: N_v x J) used in the classical ML validation step, then return the beta loss of GpBS procedure

smooth_spline(x::AbstractVector, y::AbstractVector, xnew::AbstractVector)

Perform smoothing spline on (x, y), and make predictions on xnew.

Returns: yhat, ynewhat,....

train_Gyλ(rawy::AbstractVector, rawB::AbstractMatrix, rawL::AbstractMatrix, model_file::String)

Train MLP generator G(y, λ) for λ ∈ [λl, λu] and y ~ N(f, σ²)

train_Gλ(rawy::AbstractVector, rawB::AbstractMatrix, rawL::AbstractMatrix; λl, λu)

Train MLP generator G(λ) for λ ∈ [λl, λu].

rcopy(s::Spl)

Convert RObject s.H as a Julia matrix, and s.β keeps the same.

predict(model::Spl{T}, xs::AbstractVector{T})
predict(X::Vector{Float64}, y::Vector{Float64}, J::Int, Xnew::AbstractVector{Float64}, ynew::AbstractVector{Float64}
predict(X::Vector{Float64}, y::Vector{Float64}, J::Int, Xnew::Vector{Float64}, ynew::Vector{Float64}, σ::Vector{Float64}

Make prediction based on fitted Spl on new points xs. If Xnew is provided, then also returns the prediction error ‖yhat - ynew‖_2^2.