BoundedAntiMean, BoundedIsoMean {OrdMonReg}R Documentation

Compute least square estimate of an anti- or isotonic function, bounded below and above by fixed functions

Description

This function computes the bounded least squares isotonic regression estimate, where the bounds are two functions such that the estimate is above the lower and below the upper function. To find the solution, we use the pool-adjacent-violaters algorithm for a suitable set function M, as discussed in Balabdaoui et al. (2009). The problem was initially posed in Barlow et al. (1972), including a remark (on p. 57) that the PAVA can be used to solve it. However, a formal proof is only given in Balabdaoui et al. (2009).

Usage

BoundedIsoMean(y, w, a = NA, b = NA)
BoundedAntiMean(y, w, a = NA, b = NA)

Arguments

y Vector in R^n of measurements.
w Vector in R^n of weights.
a Vector in R^n that gives lower bound.
b Vector in R^n that gives upper bound.

Details

The bounded isotonic regression problem is given by: For X = {x_1 <= ... <= x_n} let g(x_i), i = 1, ..., n be measurements of some quantity, with true mean f(x). The goal is to estimate f using least squares, i.e. to minimize

L(f) = sum_{x in X} w(x)(f(x) - g(x))^2

over all functions f that are isotonic and satisfy

a(x) <= f(x) <= b(x) mathrm{for} mathrm{all} x in X

and two fixed isotonic functions a and b. This problem can be solved using a suitable modification of the pool-adjacent-violaters algorithm, see Barlow et al. (1972, p. 57) and Balabdaoui et al. (2009).

The function BoundedAntiMean solves the same problem for antitonic curves, by simply invoking BoundedIsoMean flipping some of the arguments.

Value

The bounded isotonic (antitonic) estimate (hat f(x))_{x in X}.

Author(s)

Fadoua Balabdaoui fadoua@ceremade.dauphine.fr
http://www.ceremade.dauphine.fr/~fadoua

Kaspar Rufibach (maintainer) kaspar.rufibach@ifspm.uzh.ch
http://www.biostat.uzh.ch/aboutus/people/rufibach.html

Filippo Santambrogio filippo@ceremade.dauphine.fr
http://www.ceremade.dauphine.fr/~filippo

References

Balabdaoui, F., Rufibach, K., Santambrogio, F. (2009). Least squares estimation of two ordered antitonic regression curves. Preprint.

Barlow, R. E., Bartholomew, D. J., Bremner, J. M., Brunk, H. D. (1972). Statistical inference under order restrictions. The theory and application of isotonic regression. John Wiley and Sons, London - New York - Sydney.

See Also

The functions BoundedAntiMeanTwo and BoundedIsoMeanTwo for the problem of estimating two ordered antitonic (isotonic) regression functions. The function BoundedIsoMean depends on the function MA.

Examples

## --------------------------------------------------------
## generate data
## --------------------------------------------------------
set.seed(23041977)
n <- 35
x <- 1:n / n
f0 <- - 3 * x + 5
g0 <- 1 / (x + 0.5) ^ 2 + 1 
g <- g0 + 3 * rnorm(n)

## --------------------------------------------------------
## compute estimate
## --------------------------------------------------------
g_est <- BoundedAntiMean(g, w = rep(1 / n, n), a = -rep(Inf, n), b = f0)

## --------------------------------------------------------
## plot observations and estimate
## --------------------------------------------------------
par(mar = c(4.5, 4, 3, 0.5))
plot(0, 0, type = 'n', main = "Observations, upper bound and estimate 
    for bounded antitonic regression", xlim = c(0, max(x)), ylim = 
    range(c(f0, g)), xlab = expression(x), ylab = "observations and estimate")
points(x, g, col = 1)
lines(x, g0, col = 1, lwd = 2, lty = 2)
lines(x, f0, col = 2, lwd = 2, lty = 2)
lines(x, g_est, col = 3, lwd = 2)
legend("bottomleft", c("truth", "data", "upper bound", "estimate"), 
    lty = c(1, 0, 1, 1), lwd = c(2, 1, 2, 2), pch = c(NA, 1, NA, NA), 
    col = c(1, 1:3), bty = 'n')
    
## Not run: 
## --------------------------------------------------------
## 'BoundedIsoMean' is a generalization of 'isoMean' in the 
## package 'logcondens'
## --------------------------------------------------------
library(logcondens)
n <- 50
y <- sort(runif(n, 0, 1)) ^ 2 + rnorm(n, 0, 0.2)

isoMean(y, w = rep(1 / n, n))
BoundedIsoMean(y, w = rep(1 / n, n), a = -rep(Inf, n), b = rep(Inf, n))
## End(Not run)    

[Package OrdMonReg version 1.0.0 Index]