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 <= i="" <="n$" and="" $m[0]="" :="m[n+1]" for="" concavity.="" convexity="" (convex=TRUE), replace $>=$ by $

Usage

conreg(x, y = NULL, w = NULL, convex = FALSE, tol = 1e-07, maxit = c(200, 20), adjTol = TRUE, verbose = FALSE)

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.

optional vector of weights of the same length as x; defaults to all 1.

convex

logical indicating if convex or concave regression is desired.

tol

convergence tolerance; do not make this too small!

maxit

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

adjTol

logical indicating if the convergence test tolerance is to be adjusted (increased) in some cases.

verbose

logical indicating if knot placement iterations should be “reported”.

Value

an object of class conreg which is basically a list with componentsNote that there are several methods defined for conreg objects, see predict.conreg. Notably print and plot; also predict, residuals, fitted, knots.

Details

The algorithm is an active-set method, needing some numerical tolerance because of rounding errors in computation of finite difference rations.

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)
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)
title("Concave Regression in R")

## 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))
r.3. <- conreg(1:3,c(4:5,8))
stopifnot(resid(r.1) == 0,
          resid(r.2) == 0,
          resid(r.3) == 0,
          all.equal(fitted(r.3.),
                    c(11,17,23)/3, tol=1e-12))