BoundedAntiMean, BoundedIsoMean {OrdMonReg} | R Documentation |
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).
BoundedIsoMean(y, w, a = NA, b = NA) BoundedAntiMean(y, w, a = NA, b = NA)
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. |
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.
The bounded isotonic (antitonic) estimate (hat f(x))_{x in X}.
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
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.
The functions BoundedAntiMeanTwo
and BoundedIsoMeanTwo
for the problem of
estimating two ordered antitonic (isotonic) regression
functions. The function BoundedIsoMean
depends on the function MA
.
## -------------------------------------------------------- ## 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)