This section shows how to conduct monotone fitting with monotone splines with the existing optimization toolbox. The smoothing parameter can be tuned by cross-validation. We also compare the monotone splines with the popular smoothing splines (R's implementation smooth.spline).

using MonotoneSplines
using Plots
using Random
using RCall

Cubic Curve

n = 100
seed = 1234
σ = 0.2
Random.seed!(seed)
x = rand(n) * 2 .- 1
y = x .^3 + randn(n) * σ
λs = exp.(range(-10, -4, length = 100));

Perform cross-validation for monotone fitting with smoothing splines,

@time errs, B, L, J = MonotoneSplines.cv_mono_ss(x, y, λs)

([0.03836494968310219, 0.038356054988443196, 0.038346860326416564, 0.03833810480828669, 0.03833030336201492, 0.0383181815212615, 0.038306086526776, 0.03829277370690645, 0.03827569988291665, 0.03825256780043086  …  0.041377653962167346, 0.04164026736341345, 0.041913261668632816, 0.04219672525984907, 0.042490701739591435, 0.04279519591471907, 0.04311015485047943, 0.04343543081198881, 0.0437708461692481, 0.04411615119567598], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [3722.3861277938377 0.0 … 0.0 0.0; -5484.443813834549 2256.9009327110152 … 0.0 0.0; … ; 0.0 0.0 … 5.503905144035493 0.0; 0.0 0.0 … -5.4879608135793845 0.011543876496072695], 64)

Then plot the CV curve

scatter(log.(λs), errs, title = "seed = $seed")

Then we can choose λ which minimized the CV error.

idx = argmin(errs)
λopt = λs[idx]

0.0008846636600765526

Fit with λopt

βhat, yhat = MonotoneSplines.mono_ss(B, y, L, J, λopt);

Alternatively,

res = MonotoneSplines.mono_ss(x, y, λopt);
yhat = res.fitted
βhat = res.β

64-element Vector{Float64}:
 -1.0288757376477982
 -1.010651571190564
 -0.9869589289078865
 -0.8657114782087676
 -0.7478717956212413
 -0.6000033126726187
 -0.5064227831953195
 -0.42894948406762573
 -0.36508688261887273
 -0.3056529887751692
  ⋮
  0.7140294931741021
  0.784167204977876
  0.8180581915475537
  0.8494013724090633
  0.8659588322917697
  0.9073886912733284
  0.9556781925759968
  0.9941873855077844
  1.0025953675523138

Plot it

scatter(x, y)
scatter!(x, yhat)

We can also compare it with smooth.spline,

spl = R"smooth.spline($x, $y)"

RCall.RObject{RCall.VecSxp}
Call:
smooth.spline(x = `#JL`$x, y = `#JL`$y)

Smoothing Parameter  spar= 0.8358156  lambda= 0.001106585 (13 iterations)
Equivalent Degrees of Freedom (Df): 7.05542
Penalized Criterion (RSS): 3.174281
GCV: 0.03674491

it also determine λ by cross-validation,

λ = rcopy(R"$spl$lambda")

0.0011065854744521935

we can plot its fitting values together,

yhat_ss = rcopy(R"predict($spl, $x)$y")
scatter!(x, yhat_ss)

For ease of demonstrating other examples, we wrap up the above procedures as a function

function demo_mono_ss(x, y, λs)
    errs, B, L, J = MonotoneSplines.cv_mono_ss(x, y, λs)
    fig1 = plot(log.(λs), errs, xlab = "λ", ylab = "CV error", legend=false)
    λopt = λs[argmin(errs)]
    λ_mono_ss = [round(λopt, sigdigits = 4), round(log(λopt), sigdigits=4)]
    yhat = MonotoneSplines.mono_ss(x, y, λopt).fitted
    fig2 = scatter(x, y, label = "obs.")
    scatter!(fig2, x, yhat, label = "mono_ss (λ = $(λ_mono_ss[1]), logλ = $(λ_mono_ss[2]))")
    # ss
    spl = R"smooth.spline($x, $y)"
    λ = rcopy(R"$spl$lambda")
    λ_ss = [round(λ, sigdigits = 4), round(log(λ), sigdigits=4)]
    yhat_ss = rcopy(R"predict($spl, $x)$y")
    scatter!(fig2, x, yhat_ss, label = "ss (λ = $(λ_ss[1]), logλ = $(λ_ss[2]))")
    return plot(fig1, fig2, size = (1240, 420))
end

demo_mono_ss (generic function with 1 method)

Growth Curve

λs = exp.(range(-10, 0, length = 100));

σ = 3.0

σ = 3.0
Random.seed!(seed)
x, y, x0, y0 = MonotoneSplines.gen_data(n, σ, z->1/(1-0.42log(z)), xmin = 0, xmax = 10)
scatter(x, y)
scatter!(x0, y0)

demo_mono_ss(x, y, λs)

σ = 2.0

σ = 2.0
Random.seed!(seed)
x, y, x0, y0 = MonotoneSplines.gen_data(n, σ, z->1/(1-0.42log(z)), xmin = 0, xmax = 10)
scatter(x, y)
scatter!(x0, y0)

demo_mono_ss(x, y, λs)

σ = 0.5

σ = 0.5
Random.seed!(seed)
x, y, x0, y0 = MonotoneSplines.gen_data(n, σ, z->1/(1-0.42log(z)), xmin = 0, xmax = 10)
scatter(x, y)
scatter!(x0, y0)

demo_mono_ss(x, y, λs)

Logistic Curve

λs = exp.(range(-10, 0, length = 100));

σ = 0.2

σ = 0.2
Random.seed!(seed)
x, y, x0, y0 = MonotoneSplines.gen_data(n, σ, z->exp(z)/(1+exp(z)), xmin = -5, xmax = 5)
scatter(x, y)
scatter!(x0, y0)

demo_mono_ss(x, y, λs)

σ = 1.0

σ = 1.0
Random.seed!(seed)
x, y, x0, y0 = MonotoneSplines.gen_data(n, σ, z->exp(z)/(1+exp(z)), xmin = -5, xmax = 5)
scatter(x, y)
scatter!(x0, y0)

demo_mono_ss(x, y, λs)