smbsp: Smoothing B-spline of order n >= 0

Description

Optimize smoothing B-spline coefficients (smbsp) and knot positions (fitsmbsp) such that residual squared sum is minimized for all y columns.

Usage

smbsp(
  x,
  y,
  n = 3L,
  xki = NULL,
  nki = 1L,
  lieq = NULL,
  monotone = 0,
  positive = 0,
  mat = NULL,
  estSD = FALSE,
  tol = 1e-10
)
fitsmbsp(
  x,
  y,
  n = 3L,
  xki = NULL,
  nki = 1L,
  lieq = NULL,
  monotone = 0,
  positive = 0,
  control = list(),
  estSD = FALSE,
  tol = 1e-10
)

Value

Function, smoothing B-splines respecting optional constraints (generated by par2bsp()).

Arguments

x: Numeric vector, abscissa points
y: Numeric vector or matrix or data.frame, ordinate values to be smoothed (one set per column in case of matrix or data.frame)
n: Integer scalar, polynomial order of B-splines (3 by default)
xki: Numeric vector, strictly internal B-spline knots, i.e. lying strictly inside of x bounds. If NULL (by default), they are estimated with the help of iknots(). This vector is used as initial approximation during optimization process. Must be non decreasing if not NULL.
nki: Integer scalar, internal knot number (1 by default). When nki==0, it corresponds to polynomial regression. If xki is not NULL, this parameter is ignored.
lieq: List, equality constraints to respect by the smoothing spline, one list item per y column. By default (NULL), no constraint is imposed. Constraints are given as a 2-column matrix (xe, ye) where for each xe, an ye value is imposed. If a list item is NULL, no constraint is imposed on corresponding y column.
monotone: Numeric scalar or vector, if monotone > 0, resulting B-spline weights must be increasing; if monotone < 0, B-spline weights must be decreasing; if monotone == 0 (default), no constraint on monotonicity is imposed. If 'monotone' is a vector it must be of length ncol(y), in which case each component indicates the constraint for corresponding column of y.
positive: Numeric scalar, if positive > 0, resulting B-spline weights must be >= 0; if positive < 0, B-spline weights must be decreasing; if positive == 0 (default), no constraint on positivity is imposed. If 'positive' is a vector it must be of length ncol(y), in which case each component indicates the constraint for corresponding column of y.
mat: Numeric matrix of basis vectors, if NULL it is recalculated by bsc(). If provided, it is the responsibility of the user to ensure that this matrix be adequate to xki vector.
estSD: Logical scalar, if TRUE, indicates to calculate: SD of each y column, covariance matrix and SD of spline coefficients. All these values can be retrieved with bsppar() call (FALSE by default). These estimations are made under assumption that all y points have uncorrelated noise. Optional constraints are not taken into account of SD.
tol: Numerical scalar, relative tolerance for small singular values that should be considered as 0 if s[i] <= tol*s[1]. This parameter is ignored if estSD=FALSE (1.e-10 by default).
control: List, passed through to nlsic() call

Details

If constraints are set, we use nlsic::lsie_ln() to solve a least squares problem with equality constraints in least norm sens for each y column. Otherwise, nlsic::ls_ln_svd() is used for the whole y matrix. The solution of least squares problem is a vector of B-splines coefficients qw, one vector per y column. These vectors are used to define B-spline function which is returned as the result.

NB. When nki >= length(x)-n-1 (be it from direct setting or calculated from length(xki)), it corresponds to spline interpolation, i.e. the resulting spline will pass exactly by (x,y) points (well, up to numerical precision).

Border and external knots are fixed, only strictly internal knots can move during optimization. The optimization process is constrained to respect a minimal distance between knots as well as to bound them to x range. This is done to avoid knots getting unsorted during iterations and/or going outside of a meaningful range.

Examples

Run this code

  x=seq(0, 1, length.out=11)
  y=sin(pi*x)+rnorm(x, sd=0.1)
  # constraint B-spline to be 0 at the interval ends
  fsm=smbsp(x, y, nki=1, lieq=list(rbind(c(0, 0), c(1, 0))))
  # check parameters of found B-splines
  bsppar(fsm)
  plot(x, y) # original "measurements"
  # fine grained x
  xfine=seq(0, 1, length.out=101)
  lines(xfine, fsm(xfine)) # fitted B-splines
  lines(xfine, sin(pi*xfine), col="blue") # original function
  # visualize knot positions
  xk=bsppar(fsm)$xk
  points(xk, fsm(xk), pch="x", col="red")
 # fit broken line with linear B-splines
 x1=seq(0, 1, length.out=11)
 x2=seq(1, 3, length.out=21)
 x3=seq(3, 4, length.out=11)
 y1=x1+rnorm(x1, sd=0.1)
 y2=-2+3*x2+rnorm(x2, sd=0.1)
 y3=4+x3+rnorm(x3, sd=0.1)
 x=c(x1, x2, x3)
 y=c(y1, y2, y3)
 plot(x, y)
 f=fitsmbsp(x, y, n=1, nki=2)
 lines(x, f(x))