OrdMonReg (version 1.0.3)

BoundedAntiMean, BoundedIsoMean: Compute least square estimate of an iso- or antitonic 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 not given in Barlow et al. (1972). A short note detailing this proof is available from the authors of Balabdaoui et al. (2009) on request.

Usage

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

Value

The bounded isotonic (antitonic) estimate $(\hat g^\circ)_{i=1}^n$.

Details

The bounded isotonic regression problem is given by: For $x_1 \le \ldots \le x_n$ let $y_i, i = 1, \ldots, n$ be measurements of some quantity at the $x_i$'s, with true mean function $g^\circ(x)$. The goal is to estimate $g^\circ$ using least squares, i.e. to minimize $$L(a) = \sum_{i=1}^n w_i(y_i - a_i)^2$$ over the class of vectors $a$ that are isotonic and satisfy $$a_{L, i} \le a_i \le a_{U, i} \ \ \mathrm{for} \ \ \mathrm{all} \ \ i = 1, \ldots, n$$ and two fixed isotonic vectors $a_L$ and $a_U$. 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.

References

Balabdaoui, F., Rufibach, K., Santambrogio, F. (2009). Least squares estimation of two ordered monotone 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')
    
## --------------------------------------------------------
## '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))