This section analyzes the polarization hole data using the monotone splines techniques.

using MonotoneSplines
using Plots
using DelimitedFiles

First of all, we load the data.

current_folder = @__DIR__
data = readdlm(joinpath(current_folder, "ph.dat"));
x = data[:, 1]
y = data[:, 2]
x0 = range(minimum(x), maximum(x), length=500)

-9.14043389977092:0.018317502805152146:0.0

Then we can check how the data looks like

scatter(x, y, label = "")

Monotone Cubic Splines

Perform the monotone cubic splines with different number of basis functions $J=4, 10$

fit_mcs4 = mono_cs(x, y, 4, increasing = false)
plot!(x0, predict(fit_mcs4, x0), label = "J = 4", legend = :bottomleft)

fit_mcs10 = mono_cs(x, y, 10, increasing = false)
plot!(x0, predict(fit_mcs10, x0), label = "J = 10")

Monotone Smoothing Splines

Perform smoothing splines

yhat_ss, yhatnew_ss, _, λ = MonotoneSplines.smooth_spline(x, y, x0);

use the same $\lambda$,

fit_mss = mono_ss(x, y, λ, increasing  = false)

MonotoneSplines.MonotoneSS(-9.14043389977092, 9.14043389977092, [0.0, 1.454834583480431e-5, 0.04025396544546997, 0.040275046566675785, 0.05310974245421181, 0.07158510316354796, 0.07160005755490678, 0.08108010267692153, 0.08109742411895293, 0.08881246768312859  …  0.7734524042117233, 0.7735584982043695, 0.7899565110302229, 0.8038873718450078, 0.8040938195476999, 0.8438188519580327, 0.8438778678860593, 0.8646148315616505, 0.9414804833829471, 1.0], [0.6734819458337752 0.32650220992558177 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.9991961351507619 … 0.0 0.0], [0.9999999999999999 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 1.0], [-206209.0105682694 206209.0105682694 … 0.0 0.0; 0.0 0.0 … -51.26494840399594 51.26494840399594], [6.244923202518329e7 0.0 … 0.0 0.0; 0.0 1.1034563248616778e6 … 0.0 0.0; … ; 0.0 0.0 … 0.015373619727786093 0.0; -6.242667422657813e7 0.0 … 0.009318550356941242 -0.40422156229033135], [-0.09662692156915031, -0.09662692188778219, -0.09662692239856684, -0.0966269230243472, -0.09662692375552343, -0.09662692462374925, -0.09662692557369894, -0.09662692655791523, -0.09662692764774626, -0.0966269287382746  …  -0.09663421980402108, -0.09663422157187751, -0.09663422325006109, -0.09663422489365407, -0.09663422653413008, -0.09663422815044044, -0.09663422973921397, -0.09663423132452188, -0.09663423290349896, -0.09663423447887037], [-0.09662692167319747, -0.09662692448119378, -0.09662693006182466, -0.09662694001389979, -0.09662694676916647, -0.09662695823409825, -0.09662696846703415, -0.0966269694776969, -0.09662696947811114, -0.09662696847095167  …  -0.09662696846508866, -0.0966269694772154, -0.09662696947761472, -0.09662696846856444, -0.09662695830685132, -0.09662694731624359, -0.0966269400265264, -0.09662693006409406, -0.09662692448211467, -0.09662692188819288])

then plot it

plot!(x0, yhatnew_ss, ls = :dot, label = "Smoothing Spline (λ = $(round(λ, sigdigits = 3)))")
plot!(x0, predict(fit_mss, x0), ls = :solid, label = "Monotone Smoothing Spline (λ = $(round(λ, sigdigits = 3)))")

Monotone smoothing splines with cross-validation

Alternatively, we can find the optimal tuning parameter $\lambda$ by cross-validation,

λs = exp.(-10:0.2:1)
errs, B, L, J = cv_mono_ss(x, y, λs, increasing = false)
λopt = λs[argmin(errs)]

1.2214027581601699

Fit with the optimal tuning parameter

fit_mss2 = mono_ss(x, y, λopt, increasing = false)
plot!(x0, predict(fit_mss2, x0), label = "Monotone Smoothing Spline (λ = $(round(λopt, sigdigits = 3)))")

where the cross-validation error curve is as follows,

scatter(log.(λs), errs, label = "")