| Title: | 'Constrained Quantile Regression with Cubic B-Splines' |
|---|---|
| Description: | Quantile regression with cubic B-splines under monotonicity and convexity constraints using the Karlin-Studden SOCP formulation. The method is described in Abbes (2026) <doi:10.5281/zenodo.17427913>. This R implementation is intended for demonstration and prototyping; all B-spline and polynomial functions have been rewritten for consistency. A faster version written in 'Python' is available at <https://github.com/alexandreabbes/Constrained-Quantile-Regression-with-cubic-splines>. |
| Authors: | Alexandre Abbes [aut, cre] |
| Maintainer: | Alexandre Abbes <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1.0 |
| Built: | 2026-06-23 16:32:03 UTC |
| Source: | https://github.com/cran/BsplineQuantReg |
Applies Karlin-Studden SOCP constraints to ensure positivity of a quadratic polynomial on the interval [0,1].
apply_karlin_constraints(p2, p1, p0, z0, verbose = FALSE)apply_karlin_constraints(p2, p1, p0, z0, verbose = FALSE)
p2 |
Coefficient of u^2 |
p1 |
Coefficient of u |
p0 |
Constant term |
z0 |
Auxiliary SOCP variable |
verbose |
boolean FALSE (default) or TRUE. |
List of CVXR constraints
Computes the values of all B-spline basis functions at given points.
bs_direct(Basis, xvalues)bs_direct(Basis, xvalues)
Basis |
Object returned by |
xvalues |
Vector of evaluation points |
Matrix of basis function values (n_splines x length(xvalues))
Build B-spline basis in piecewise polynomial form Computes local polynomial coefficients for each B-spline basis function on each interval. Polynomials are expressed in the canonical basis Uses De Boor's recursion formula.
Bspline_base(sn, degree = 3, der = 0, verbose = FALSE)Bspline_base(sn, degree = 3, der = 0, verbose = FALSE)
sn |
Extended knot vector (including endpoint repetitions) This means if t0..tkn it the set of knots then sn should be given as a vector with "degree" times t_0 and t_kn at the begining and the ends. its length is number of intervals+1+2*degree. |
degree |
B-spline degree (default = 3 for cubic) |
der |
Derivative order (0 = original basis) |
verbose |
boolean FALSE (default) or TRUE. |
A list containing:
base |
Coefficients in the local bases, in the form of an 3-d array [j,nu,coeff], j : the number of the spline in the basis, nu: the number of the interval in the extended notation, coeff : the coefficients in decreasing order convention on the local bases (t-s_nu)^l, l=3..1 base[j,,] is a matrix of piecewise polynomial function compatible with the pp-form. |
base0 |
Coefficients in canonical basis (centered at 0) (1, t, t^2, t^3) centered at the interval origin. |
knots |
Extended knot vector |
int_knots |
Internal knots (effective partition including ends) |
degree |
Spline degree |
n_splines |
Number of basis functions |
deriv_order |
Applied derivative order |
sn <- c(0,0,0,0,1,2,3,4,5,5,5,5) basis <- Bspline_base(sn, degree=3) x=(0:(5*100))/100 y=bs_direct(basis,x) matplot(x,t(y)) #or simple: #view_basis(basis)sn <- c(0,0,0,0,1,2,3,4,5,5,5,5) basis <- Bspline_base(sn, degree=3) x=(0:(5*100))/100 y=bs_direct(basis,x) matplot(x,t(y)) #or simple: #view_basis(basis)
Computes the basis of order der derivatives of a B-spline basis.
Bspline_deriv(bspline, der = 2, verbose = FALSE)Bspline_deriv(bspline, der = 2, verbose = FALSE)
bspline |
Object returned by |
der |
Derivative order |
verbose |
boolean FALSE (default) or TRUE. |
A list similar to Bspline_base for the derivative basis
Converts a B-spline basis to normalized first and second derivative coefficients on each interval.
bspline_to_deriv_coeffs_pp(tn, degree = 3, xvalues = 0, verbose = FALSE)bspline_to_deriv_coeffs_pp(tn, degree = 3, xvalues = 0, verbose = FALSE)
tn |
Knot vector (effective partition, not extended) |
degree |
Spline degree (default = 3) |
xvalues |
Evaluation points for design matrix (0 = no evaluation) |
verbose |
boolean FALSE (default) or TRUE. |
A list containing:
d0 |
Design matrix (if xvalues provided) |
d1 |
First derivative coefficients [a3, a2, a1] for each interval |
d2 |
Second derivative values at knots |
Converts a polynomial expressed in the basis (t-a)^k to its representation in the basis (t-b)^k using Taylor's formula. P(t) = sum c_k (t-a)^k P(t) = sum c"_k (t-b)^k
change_polynomial_base_taylor(coeffs_a, a, b)change_polynomial_base_taylor(coeffs_a, a, b)
coeffs_a |
Coefficients in basis centered at a (decreasing powers) |
a |
Original expansion point |
b |
New expansion point |
Coefficients in basis centered at b
Check if package is beta version
is_beta()is_beta()
TRUE if beta version
Get package version
package_version()package_version()
Current version string
Evaluates a polynomial at one or more points.
poly_eval(p, xvalues)poly_eval(p, xvalues)
p |
Coefficient vector (decreasing powers) |
xvalues |
vector at which to evaluate the polynomial |
Vector of same length as xvalues : polynomial values at the points xvalues
# P(x) = 1 + x + x^2 poly_eval(c(1, 1, 1), c(0, 1, 2)) # returns c(1, 3, 7)# P(x) = 1 + x + x^2 poly_eval(c(1, 1, 1), c(0, 1, 2)) # returns c(1, 3, 7)
Adds two polynomials represented by coefficients in decreasing power order.
polyadd(p1, p2, verbose = FALSE)polyadd(p1, p2, verbose = FALSE)
p1 |
First polynomial (coefficient vector) |
p2 |
Second polynomial (coefficient vector) |
verbose |
boolean FALSE (default) or TRUE. |
Coefficient vector of the sum
polyadd(c(1, 1), c(1, -1)) # returns c(2, 0)polyadd(c(1, 1), c(1, -1)) # returns c(2, 0)
der of a polynomial.Polynomial derivative
Computes the derivative of order der of a polynomial.
polyderiv(p, der = 1)polyderiv(p, der = 1)
p |
Coefficient vector (decreasing powers) |
der |
Derivative order (default = 1) |
Coefficients of the derivative polynomial
# P(x) = x^2 -> P'(x) = 2x polyderiv(c(1, 0, 0), 1) # returns c(2, 0)# P(x) = x^2 -> P'(x) = 2x polyderiv(c(1, 0, 0), 1) # returns c(2, 0)
Multiplies two polynomials represented by coefficients in decreasing power order.
polymul(p1, p2, ord = 0, verbose = FALSE)polymul(p1, p2, ord = 0, verbose = FALSE)
p1 |
First polynomial (coefficient vector, decreasing powers) |
p2 |
Second polynomial (coefficient vector, decreasing powers) |
ord |
Unused (compatibility parameter) |
verbose |
boolean FALSE (default) or TRUE. |
Coefficient vector of the product polynomial (decreasing powers)
# (1 + x) * (1 + x) = 1 + 2x + x^2 polymul(c(1, 1), c(1, 1)) # returns c(1, 2, 1)# (1 + x) * (1 + x) = 1 + 2x + x^2 polymul(c(1, 1), c(1, 1)) # returns c(1, 2, 1)
Removes leading zeros from a polynomial coefficient vector.
reduce_pol(p, verbose = FALSE)reduce_pol(p, verbose = FALSE)
p |
Polynomial coefficient vector(coef in decreasing order) |
verbose |
boolean FALSE (default) or TRUE. |
Reduced vector (without leading zeros)
reduce_pol(c(0,0, 1, 1))reduce_pol(c(0,0, 1, 1))
Computes derivative values of a B-spline at knots (efficient because it directly uses polynomial coefficients).
Spline_der_knots(Bspline, der = 1)Spline_der_knots(Bspline, der = 1)
Bspline |
Object returned by |
der |
Derivative order (default = 1) |
Matrix of derivative values (n_splines x n_knots)
Evaluates a spline (linear combination of B-splines) at given points.
spline_eval(Bspline, xvalues)spline_eval(Bspline, xvalues)
Bspline |
Spline object (list with coefficients on the Bspline basis, degree, extended knots) |
xvalues |
Vector of evaluation points |
vector of same length as xvalues, with Spline values at the requested points
{ # Create and evaluate a spline sn <- c(0,0,0,0,1,2,3,4,5,5,5,5) basis <- Bspline_base(sn, degree=3) basis$coefficients <- runif(basis$n_splines) y <- spline_eval(basis, seq(0,5,length=100))}{ # Create and evaluate a spline sn <- c(0,0,0,0,1,2,3,4,5,5,5,5) basis <- Bspline_base(sn, degree=3) basis$coefficients <- runif(basis$n_splines) y <- spline_eval(basis, seq(0,5,length=100))}
Performs quantile regression using cubic B-splines, with optional monotonicity constraints (via Karlin-Studden) and convexity constraints.
SplineConstQuantRegBs3( xtab, ytab, knots, tau, monot = 0, convcons = 0, solver = "CLARABEL", weight = NULL, verbose = FALSE )SplineConstQuantRegBs3( xtab, ytab, knots, tau, monot = 0, convcons = 0, solver = "CLARABEL", weight = NULL, verbose = FALSE )
xtab |
Predictor vector (x) |
ytab |
Response vector (y) |
knots |
Knot vector or number of knots (quantiles are then used) |
tau |
Quantile (between 0 and 1) |
monot |
Monotonicity constraint vector per interval: 1 = increasing, -1 = decreasing, 0 = unconstrained. If scalar, repeated. |
convcons |
Convexity constraint vector per knot: 1 = convex, -1 = concave, 0 = unconstrained. If scalar, repeated. |
solver |
CVXR solver to use (default = "CLARABEL") |
weight |
Observation weights (default = 1 for all) |
verbose |
boolean FALSE (default) or TRUE. |
A list containing:
coefficients |
B-spline coefficients (including y mean) |
degree |
Spline degree (always 3) |
knots |
Knot vector used |
int_knots |
Same as knots (compatibility) |
Abbes, A. (2025). Quantile regression with cubic polynomial splines under shape constraints with applications . Zenodo. doi:10.5281/zenodo.16999784
de Boor, C. (1978). A Practical Guide to Splines. Springer-Verlag. doi:10.1007/978-1-4612-6333-3
Karlin, S., & Studden, W. J. (1966). Tchebycheff Systems: With Applications in Analysis and Statistics. Interscience.
Koenker, R., & Bassett, G. (1978). Regression Quantiles. Econometrica, 46(1), 33-50. doi:10.2307/1913643
Koenker, R. (2025). quantreg: Quantile Regression. R package version 5.99. https://CRAN.R-project.org/package=quantreg
Ng, P., & Maechler, M. (2024). cobs: Constrained B-Splines. R package version 1.3-8. https://CRAN.R-project.org/package=cobs
Related R packages:
quantreg - Quantile regression with linear programming
cobs - Constrained B-sines (linear and quadratic only)
Other implementations:
MATLAB/Python versions: https://github.com/alexandreabbes/Constrained-Quantile-Regression-with-cubic-splines
#optional set.seed(42) x <- seq(0, 1, length=100) y <- 2*x + sin(6*pi*x)/2 + rnorm(100, 0, 0.05) knots <- quantile(x, probs=seq(0,1,length.out=10)) # Median quantile regression without constraints fit <- SplineConstQuantRegBs3(x, y, knots, tau=0.5) # With increasing monotonicity constraint fit_monot <- SplineConstQuantRegBs3(x, y, knots, tau=0.5, monot=1) # With convexity constraint fit_convex <- SplineConstQuantRegBs3(x, y, knots, tau=0.5, convcons=1)#optional set.seed(42) x <- seq(0, 1, length=100) y <- 2*x + sin(6*pi*x)/2 + rnorm(100, 0, 0.05) knots <- quantile(x, probs=seq(0,1,length.out=10)) # Median quantile regression without constraints fit <- SplineConstQuantRegBs3(x, y, knots, tau=0.5) # With increasing monotonicity constraint fit_monot <- SplineConstQuantRegBs3(x, y, knots, tau=0.5, monot=1) # With convexity constraint fit_convex <- SplineConstQuantRegBs3(x, y, knots, tau=0.5, convcons=1)
Runs quantile regression tests with and without constraints, and displays results. Demo function.
test_karlin_simple(verbose = FALSE, seed = NULL)test_karlin_simple(verbose = FALSE, seed = NULL)
verbose |
boolean FALSE (default) or TRUE. |
seed |
(default=NULL) value for the random generator. |
No return value, produces plots.
test_karlin_simple()test_karlin_simple()
Plots all functions of a B-spline basis.
view_basis(Bspline, xvalues = 0)view_basis(Bspline, xvalues = 0)
Bspline |
Object returned by |
xvalues |
Vector of evaluation points for plotting (by default 100 points are computed in the knot range) |
No return value, called for side effects (generates a plot)