2
$\begingroup$

I'd like an automatic way to find the "best" smoothing factor s for a spline fit to a given set of data points. Here's a sample visualization of some data and the fit splines for various s values:

enter image description here

In this case, clearly s=2 (and at a lesser degree s=1) is not a good fit. On the other hand s=0.5 fits the data almost as good as s=0.1 but with less than half the number of knots and thus is less susceptible to overfitting. So my question is, what's a robust method to determine the "optimal", or at least a good enough s to fit the data?

$\endgroup$
7
  • 1
    $\begingroup$ You might consider cross-validation approaches. $\endgroup$
    – Glen_b
    Apr 12, 2018 at 5:59
  • $\begingroup$ Are you working in any particular programming language? Also, what type of spline are you working with? $\endgroup$
    – Jim
    Apr 12, 2018 at 6:16
  • $\begingroup$ I'm using scipy's UnivariateSpline $\endgroup$
    – gsakkis
    Apr 13, 2018 at 6:49
  • $\begingroup$ If you iterate over a range of smoothing factors, each time calculating R-squared (R2) as "R2 = 1.0 - (absolute_error_variance / dependent_data_variance)" (this would run quickly) and then plot smoothing factor vs. R-squared, you should be able to see where the approximate optimal smoothing factor is on that plot. This is easy to program and might be worth a try as it is simple to do. For a numpy array named X, the variance is X.var() so the R-squared calculation becomes trivial with numpy. $\endgroup$ Apr 14, 2018 at 15:06
  • $\begingroup$ Wow, the content in the scipy doc is incorrect. This is certainly not a smoothing spline. Could you post the data? I would like to provide at least a clean R solution. $\endgroup$
    – Jim
    Apr 14, 2018 at 19:04

2 Answers 2

1
$\begingroup$

What is a smoothing spline?

The Wikipedia article on smoothing splines does a good job in explaining that. To recap, given a set of data points, $\{ (x_i, y_i)_{i=1}^n \}$, a smoothing spline is a solution to the interpolation problem:

$$\underset{f}{\arg\min} \sum_{i=1}^n (y_i - f(x_i))^2 + \lambda \int_{x_{(1)}}^{x_{(n)}} f''(x)^2 dx,$$

with $f$ constrained to be piecewise cubic between different $x_i$. The first part measures the goodness of fit of such an $f$ to the observed data. The second part is a penalty term for the wiggliness (non-smoothness) of $f$.

Leaving it to us to find a good trade-off between fit and smoothness by means of $\lambda$.


Smoothing splines in R

Luckily R has the splines package that does the heavy lifting for us.

library(splines)

mydata <- read.csv(...)

myspline <- smooth.spline(x = mydata$x, y = mydata$y
                          , lambda = 8e-9 # optim 8.332658e-11
                          , cv = TRUE) 

xgrid <- sort(union(mydata$x
             , seq(from = min(mydata$x), to = max(mydata$x), by = 1))
             , decreasing = FALSE)

yhat_xgrid <- predict(myspline, x = xgrid)$y

plot(x = mydata$x, y = mydata$y, log = "x", ylim = c(0,1)
     , xlab = "x (log-scale)", ylab = "y"
     , col  = "lightblue", pch = 19)
lines(x = xgrid, y = yhat_xgrid, type = "l", col = "darkorange")
grid()
legend(...)

And we obtain this lovely plot.

Smoothing spline $\lambda = 8 \cdot 10^{-9}$.

The optimal values for $\lambda$ are $\hat{\lambda}^*_{\text{LOO}} = 8.33 \cdot 10^{-11}$ and $\hat{\lambda}^*_{\text{GCV}} = 5.81 \cdot 10^{-13}$. I like the one plotted: $\hat{\lambda}^*_{\text{Jim}} = 8 \cdot 10^{-9}$.

$\endgroup$
1
  • 1
    $\begingroup$ In stackoverflow.com/q/55363639/10698244 this post is referred to and it is asked for a smoothing spline implementation in Python. As it seems, there are different thoughts whether UnivariateSpline is a smoothing spline or not. Perhaps you could participate in the discussion there? $\endgroup$ Mar 28, 2019 at 20:17
0
$\begingroup$

In SciPy s-factor selection depends on the data scale. I use normalized s-factor in SplineCloud to avoid this ambiguity.

Basically, SplineCloud provides an interactive version of a smoothing splines method but extends it to the parametric splines. You can reuse obtained curves in your code using API and a client library for Python.

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.