conreg: Convex / Concave Regression

Description

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 < \ldots < 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}) \ge (m_{i+1} - m_i)/(x_{i+1} - x_i),$$ for $1 \le i \le n$ and $m_0 := m_{n+1} := -\infty$, for concavity. For convexity (convex=TRUE), replace $\ge$ by $\le$ and $-\infty$ by $+\infty$.

Usage

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)

Value

an object of class conreg which is basically a list with components

x: sorted (and possibly aggregated) abscissa values x.
y: corresponding y values.
w: corresponding weights, only for "Duembgen06_R".
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 call to conreg() used.
iter: integer (vector of length one or two) with the number of iterations used (in the outer and inner loop for "Duembgen06_R").

Note that there are several methods defined for conreg objects, see predict.conreg or methods(class = "conreg").

Notably print and plot; also

predict, residuals, fitted,

knots.

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.

Arguments

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 xy.coords.

w

for method "Duembgen06_R" only: optional vector of weights of the same length as x; defaults to all 1.

convex

logical indicating if convex or concave regression is desired.

method

a character string indicating the method used,

"Duembgen06_R": is an active set method written by Lutz Duembgen (University of Berne, CH) in Matlab in 2006 and translated to R by Martin Maechler.

"SR"

is an R interface to the C code of a Support Reduction algorithm written by Piet Groeneboom (TU Delft, NL) and donated to the cobs package in July 2012.

tol

convergence tolerance(s); do not make this too small!

maxit

maximal number of (outer and inner) iterations of knot selection.

adjTol

(for "Duembgen06_R" only:) logical indicating if the convergence test tolerance is to be adjusted (increased) in some cases.

verbose

logical or integer indicating if (and how much) knot placement and fitting iterations should be “reported”.

Author

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.

Details

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.

Examples

Run this code


## 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)))

Run the code above in your browser using DataLab