Convex / Concave Regression
Compute a univariate concave or convex regression, i.e., for given vectors, x,y,w in R^n, where x has to be strictly sorted (x_1 < x_2 < … < x_n), compute an n-vector m minimizing the weighted sum of squares sum(i=1..n; w_i * (y_i - m_i)^2) under the constraints
(m[i] - m[i-1])/(x[i] - x[i-1]) >= (m[i+1] - m[i])/(x[i+1] - x[i]),
for 1 <= i <= n and
m[0] := m[n+1] := -Inf,
for concavity.
For convexity (convex=TRUE), replace >= by
<= and -Inf by +Inf.
conreg(x, y = NULL, w = NULL, convex = FALSE,
method = c("Duembgen06_R", "SR"),
tol = c(1e-10, 1e-7), maxit = c(500, 20),
adjTol = TRUE, verbose = FALSE)x, y |
numeric vectors giving the values of the predictor and
response variable. Alternatively a single “plotting”
structure (two-column matrix / y-values only / list, etc) can be
specified: see |
w |
for |
convex |
logical indicating if convex or concave regression is desired. |
method |
a character string indicating the method used,
|
tol |
convergence tolerance(s); do not make this too small! |
maxit |
maximal number of (outer and inner) iterations of knot selection. |
adjTol |
(for |
verbose |
logical or integer indicating if (and how much) knot placement and fitting iterations should be “reported”. |
Both algorithms need some numerical tolerances because of rounding
errors in computation of finite difference ratios.
The active-set "Duembgen06_R" method notably has two different
such tolerances which were both 1e-7= 10^{7} up to March
2016.
The two default tolerances (and the exact convergence checks) may change in the future, possibly to become more adaptive.
an object of class conreg which is basically a list with components
x |
sorted (and possibly aggregated) abscissa values |
y |
corresponding y values. |
w |
corresponding weights, only for |
yf |
corresponding fitted values. |
convex |
logical indicating if a convex or a concave fit has been computed. |
iKnots |
integer vector giving indices of the knots, i.e. locations where the fitted curve has kinks. Formally, these are exactly those indices where the constraint is fulfilled strictly, i.e., those i where (m[i] - m[i-1])/(x[i] - x[i-1]) > (m[i+1] - m[i])/(x[i+1] - x[i]). |
call |
the |
iter |
integer (vector of length one or two) with the number of
iterations used (in the outer and inner loop for |
Note that there are several methods defined for conreg objects,
see predict.conreg or methods(class = "conreg").
Also, interpSplineCon() to construct a more smooth
(cubic) spline, and isIsplineCon() which checks
if the int is strictly concave or convex the same as the
conreg() result from which it was constructed.
Lutz Duembgen programmed the original Matlab code in July 2006; Martin Maechler ported it to R, tested, catch infinite loops, added more options, improved tolerance, etc; from April 27, 2007.
isoreg for isotone (monotone) regression;
CRAN packages ftnonpar, cobs, logcondens.
## Generated data :
N <- 100
f <- function(X) 4*X*(1 - X)
xx <- seq(0,1, length=501)# for plotting true f()
set.seed(1)# -> conreg does not give convex cubic
x <- sort(runif(N))
y <- f(x) + 0.2 * rnorm(N)
plot(x,y, cex = 0.6)
lines(xx, f(xx), col = "blue", lty=2)
rc <- conreg(x,y)
lines(rc, col = 2, force.iSpl = TRUE)
# 'force.iSpl': force the drawing of the cubic spline through the kinks
title("Concave Regression in R")
y2 <- y
## Trivial cases work too:
(r.1 <- conreg(1,7))
(r.2 <- conreg(1:2,7:6))
(r.3 <- conreg(1:3,c(4:5,1)))Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.