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.038395854934803304, 0.03839375928415281, 0.03838826738987779, 0.0383809710334689, 0.03837032302482713, 0.03835874092192127, 0.038343576755696655, 0.03832796392734008, 0.03830780462522126, 0.038275376136259295  …  0.04137759437681002, 0.04164021277916051, 0.041913206221465035, 0.042196674158723774, 0.042490657941377304, 0.0427951479690365, 0.0431100944414775, 0.04343538190253497, 0.043770785937548115, 0.04411610106206591], [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], [7049.07849620016 0.0 … 0.0 0.0; -9701.411896739612 2348.196637868324 … 0.0 0.0; … ; 0.0 0.0 … 5.534829241098749 0.0; 0.0 0.0 … -5.518834323116495 0.011477206321819491], 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.0288973266267285
 -1.0169967436017602
 -0.9933147704827087
 -0.872131550938475
 -0.7482172252734769
 -0.6003933892529125
 -0.5091857299915032
 -0.4314585858082689
 -0.3829045954368462
 -0.3418844599722128
  ⋮
  0.7071221238215786
  0.7630132296751112
  0.8023649437876629
  0.8337982230387274
  0.8622202511290816
  0.9051040121500407
  0.9533887247878259
  0.9941865936505095
  1.0025945371846579

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.001106585471596121

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)