BoundedAntiMeanTwo, BoundedIsoMeanTwo {OrdMonReg}R Documentation

Compute solution to the problem of two ordered antitonic (or isotonic) curves

Description

See details below.

Usage

BoundedAntiMeanTwo(g1, w1, g2, w2, K1 = 1000, K2 = 400, 
    delta = 10^(-4), errorPrec = 10, output = TRUE)
BoundedIsoMeanTwo(g1, w1, g2, w2, K1 = 1000, K2 = 400, 
    delta = 10^(-4), errorPrec = 10, output = TRUE)

Arguments

g1 Vector in R^n, measurements of upper function.
w1 Vector in R^n, weights for upper function.
g2 Vector in R^n, measurements of lower function.
w2 Vector in R^n, weights for lower function.
K1 Upper bound on number of iterations.
K2 Number of iterations where step length is changed from the inverse of the norm of the subgradient to a diminishing function of the norm of the subgradient.
delta Upper bound on the error, defines stopping criterion. See Table 1 in Balabdaoui et al. (2009).
errorPrec Computation of stopping criterion is expensive. Therefore, the stopping criterion is only evaluated at every errorPrec-th iteration of the algorithm.
output Should intermediate results be output?

Details

We consider the problem of estimating two antitonic (isotonic) regression curves g_1^* and g_2^* under the constraint that g_1^* >= g_2^*. Given two sets of n data points g_1(x_1), ..., g_1(x_n) and g_2(x_1), ..., g_2(x_n) that are observed at (the same) deterministic points x_1, ..., x_n with weight functions w_1 and w_2, respectively, the estimates are obtained by minimizing the Least Squares criterion

L(f_1, f_2) = sum_{i=1}^n (g_1(x_i) - f_1(x_i))^2 w_1(x_i)+ sum_{i=1}^n (g_2(x_i) - f_2(x_i))^2 w_2(x_i)

over the class of pairs of functions (f_1, f_2) such that f_1 and f_2 are antitonic and f_1(x_i) >= f_2(x_i) for all i = {1, ..., n}. The estimates are computed with an projected subgradient algorithm where the projection is calculated using a suitable version of the pool-adjacent-violaters algorithm (PAVA).

The function BoundedIsoMeanTwo solves the same problem for isotonic curves, by simply invoking BoundedAntiMeanTwo and suitably flipping some of the arguments.

Value

g1 The estimated function g_1^*.
g2 The estimated function g_2^*.
L Value of the least squares criterion at the minimum.
error Value of error.
k Number of iterations performed.
tau Step length at final iteration.

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.

See Also

The functions BoundedAntiMean and BoundedIsoMean for the problem of estimating one antitonic (isotonic) regression function bounded above and below by fixed functions. The function BoundedAntiMeanTwo depends on the functions BoundedAntiMean, bstar_n, LSfunctional, and Subgradient.

Examples


## ========================================================
## The first example uses simulated data
## For the analysis of the mechIng dataset see below
## ========================================================

## --------------------------------------------------------
## initialization
## --------------------------------------------------------
set.seed(23041977)
n <- 100
x <- 1:n
g1 <- 1 / x^2 + 2
g1 <- g1 + 3 * rnorm(n)
g2 <- 1 / log(x+3) + 2
g2 <- g2 + 4 * rnorm(n)
w1 <- runif(n)
w1 <- w1 / sum(w1)
w2 <- runif(n)
w2 <- w2 / sum(w2)

## --------------------------------------------------------
## compute estimates
## --------------------------------------------------------
res <- BoundedAntiMeanTwo(g1, w1, g2, w2, errorPrec = 20)
    
## corresponding isotonic problem
res2 <- BoundedIsoMeanTwo(-g1, w1, -g2, w2, errorPrec = 20)
    
## the following vectors are equal
res$g1 - -res2$g2    
res$g2 - -res2$g1
        
## --------------------------------------------------------
## Checking of solution
## --------------------------------------------------------
# This compares the first component of res$g1 with a^*_1:
c(res$g1[1], astar_1(g1, w1, g2, w2))

## --------------------------------------------------------
## plot original functions and estimates
## --------------------------------------------------------
par(mfrow = c(1, 1), mar = c(4.5, 4, 3, 0.5))
plot(x, g1, col = 2, main = "Original observations and estimates in problem 
two ordered antitonic regression functions", xlim = c(0, max(x)), ylim = 
range(c(res$g1, res$g2, g1, g2)), xlab = expression(x), 
ylab = "measurements and estimates")
points(x, g2, col = 3)
lines(x, res$g1 + 0.01, col = 2, type = 's', lwd = 2)
lines(x, res$g2 - 0.01, col = 3, type = 's', lwd = 2)
legend("bottomleft", c(expression("upper estimated function g"[1]*"*"), 
    expression("lower estimated function g"[2]*"*")), lty = 1, col = 2:3, 
    lwd = 2, bty = "n")

## ========================================================
## Analysis of the mechIng dataset
## ========================================================

## --------------------------------------------------------
## input data
## --------------------------------------------------------
data(mechIng)

## for quick analysis only use randomly chosen 200 observations
set.seed(23041977)
mechIng <- mechIng[sort(sample(1:dim(mechIng)[1])[1:200]), ]

x <- mechIng$x
n <- length(x)
g1 <- mechIng$g1
g2 <- mechIng$g2
w1 <- rep(1, n)
w2 <- w1

## --------------------------------------------------------
## compute unordered estimates
## --------------------------------------------------------
g1_pava <- BoundedIsoMean(y = g1, w = w1, a = NA, b = NA)
g2_pava <- BoundedIsoMean(y = g2, w = w2, a = NA, b = NA)

## --------------------------------------------------------
## compute estimates
## --------------------------------------------------------
res1 <- BoundedIsoMeanTwo(g1, w1, g2, w2, K1 = 1000, K2 = 400, 
    delta = 10^-4, errorPrec = 20, output = TRUE)

## --------------------------------------------------------
## compute smoothed versions
## --------------------------------------------------------
g1_mon <- res1$g1
g2_mon <- res1$g2   

kernel <- function(x, X, h, Y){
    tmp <- dnorm((x - X) / h) 
    res <- sum(Y * tmp) / sum(tmp)
    return(res)
    }
h <- 0.1 * n^(-1/5)

g1_smooth <- rep(NA, n)
g2_smooth <- g1_smooth
for (i in 1:n){
    g1_smooth[i] <- kernel(x[i], X = x, h, g1_mon)
    g2_smooth[i] <- kernel(x[i], X = x, h, g2_mon)
}
            
## --------------------------------------------------------
## plot original functions and estimates
## --------------------------------------------------------
par(mfrow = c(2, 1), oma = c(0, 0, 2, 0), mar = c(4.5, 4, 2, 0.5), 
    cex.main = 0.8, las = 1) 

plot(0, 0, type = 'n', xlim = c(0, max(x)), ylim = 
    range(c(g1, g2, g1_mon, g2_mon)), xlab = "x", ylab = 
    "measurements and estimates", main = "ordered antitonic estimates")
points(x, g1, col = grey(0.3), pch = 20, cex = 0.8)
points(x, g2, col = grey(0.6), pch = 20, cex = 0.8)
lines(x, g1_mon + 0.1, col = 2, type = 's', lwd = 3)
lines(x, g2_mon - 0.1, col = 3, type = 's', lwd = 3)
legend(0.2, 10, c(expression("upper isotonic function g"[1]*"*"), 
    expression("lower isotonic function g"[2]*"*")), lty = 1, col = 2:3, 
    lwd = 3, bty = "n")

plot(0, 0, type = 'n', xlim = c(0, max(x)), ylim = 
    range(c(g1, g2, g1_mon, g2_mon)), xlab = "x", ylab = "measurements and 
    estimates", main = "smoothed ordered antitonic estimates")
points(x, g1, col = grey(0.3), pch = 20, cex = 0.8)
points(x, g2, col = grey(0.6), pch = 20, cex = 0.8)
lines(x, g1_smooth + 0.1, col = 2, type = 's', lwd = 3)
lines(x, g2_smooth - 0.1, col = 3, type = 's', lwd = 3)
legend(0.2, 10, c(expression("lower isotonic smoothed function "*tilde(g)[1]*"*"), 
    expression("lower isotonic smoothed function "*tilde(g)[2]*"*")), 
    lty = 1, col = 2:3, lwd = 3, bty = "n")

par(cex.main = 1)
title("Original observations and estimates in mechanical engineering example", 
    line = 0, outer = TRUE)

[Package OrdMonReg version 1.0.0 Index]