Fused regression, graphs and trend filtering

using HighDimensionalModels, Random, LinearAlgebra, Statistics
import KernelSmoothing
rng = Xoshiro(80)

Random.Xoshiro(0x33845dd4fd4c3729, 0x88e4542925b33b71, 0x51a654e24e1c4041, 0xde60537f28d2ffcf, 0x7b134b3612b8296a)

Fused regression on ordered coefficients

X = randn(rng, 80, 8)
beta = [0.,0,1,1,1,-.5,-.5,0]
y = X*beta + .3randn(rng,80)
levels = [.2,.08,.03,.01]
path = fused_lasso_path(X,y,levels;l1_ratio=.25)
@assert all(s -> s.converged, path.solutions)
println("Fused degrees of freedom: ", dof(path))
cv = cv_fused_lasso(X,y,levels;l1_ratio=.25,nfolds=4,rng=Xoshiro(1),threads=false)
println("CV level: ", cv.λmin)
D = graph_difference(8, [(j,j+1) for j in 1:7])
gpath = GeneralizedLassoPath(X,y,D,levels)
println("Graph BIC: ", bic(gpath))

Fused degrees of freedom: [4, 4, 6, 7]
CV level: 0.03
Graph BIC: [-87.21501940627665, -155.8201515267873, -171.86802981862806, -167.05795618798766]

Compare trend filtering and local linear regression on the same data

x = collect(1.:80)
truth = @. .02*x + .03*max(x-25,0) - .06*max(x-55,0)
signal = truth + .1randn(rng,80)
tpath = trend_filter_path(signal,[.03,.01,.003,.001];order=1)
valid = findall(isfinite,bic(tpath))
j = valid[argmin(bic(tpath)[valid])]
tf = coef(tpath,tpath.λpath[j])
lp = KernelSmoothing.local_polynomial(x,signal,x;bandwidth=8.,degree=1)
println("Trend filtering MSE on this signal: ", mean(abs2,tf-truth))
println("Local linear MSE on this signal: ", mean(abs2,lp-truth))
@assert all(isfinite,tf) && all(isfinite,lp)

Trend filtering MSE on this signal: 0.0012088422916674251
Local linear MSE on this signal: 0.001345486724327882