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.040255767815961296, 0.040275046566675785, 0.053120406912553514, 0.07158510316354796, 0.07160151251650909, 0.08108010267692153, 0.08880327201348134, 0.08881246768312859  …  0.7735049160513184, 0.7735584982043695, 0.7900313672860578, 0.8038873718450078, 0.8041381144717783, 0.8438188519580327, 0.8519292828788234, 0.8646148315616505, 0.9417901105035358, 1.0], [0.6734819458337752 0.32650221063496454 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.9991961711349227 … 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.53763434273872 51.53763434273872], [6.242667422805045e7 0.0 … 0.0 0.0; -6.24492310139133e7 428.8793536648739 … 0.0 0.0; … ; 0.0 0.0 … 71.43886038880825 0.0; 0.0 0.0 … -70.05267748833886 0.011036672741150104], [-0.00010540172415586472, -0.00010545870524991872, -0.00024925570161359484, -0.0003172203598952579, -0.00031793711533506633, -0.000318100045055394, -0.0007954049898739803, -0.000795468937069027, -0.0007955053865183091, -0.0007955411073331438  …  -0.47302855217086537, -0.5926980987370047, -0.5926980988756441, -0.5926980993424537, -1.001966027524817, -1.0019660276424656, -1.0019660280184943, -1.3849584959992836, -1.384958497054536, -1.3849584977221798], [-0.00010542260778072007, -0.00039813937290813833, -0.0007955739462636377, -0.0018508856542169736, -0.0027937719583265358, -0.0029818678046731302, -0.005703884111338646, -0.006235107368740723, -0.006235107375298541, -0.0057070183314291995  …  -0.00570232674677156, -0.006235107359004756, -0.006235107367300503, -0.005705108724883886, -0.002981868442153324, -0.0027937789471503167, -0.0018509087654751675, -0.0007955763520098499, -0.00039833929452060083, -0.0001055743044414837])

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)]

4.5399929762484854e-5

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 = "")