Description Usage Arguments Details Value Author(s) References See Also Examples
See details below.
1 2 3 4 | BoundedIsoMeanTwo(g1, w1, g2, w2, K1 = 1000, K2 = 400,
delta = 10^(-4), errorPrec = 10, output = TRUE)
BoundedAntiMeanTwo(g1, w1, g2, w2, K1 = 1000, K2 = 400,
delta = 10^(-4), errorPrec = 10, output = TRUE)
|
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. |
errorPrec |
Computation of stopping criterion is expensive. Therefore, the stopping criterion is
only evaluated at every |
output |
Should intermediate results be output? |
We consider the problem of estimating two isotonic (antitonic) regression curves g_1^\circ and g_2^\circ under the constraint that g_1^\circ ≤ g_2^\circ. Given two sets of n data points y_1, …, y_n and z_1, …, z_n that are observed at (the same) deterministic design points x_1, …, x_n with weights w_{1,i} and w_{2,i}, respectively, the estimates are obtained by minimizing the Least Squares criterion
L_2(a, b) = ∑_{i=1}^n (y_i - a_i)^2 w_{1,i} + ∑_{i=1}^n (z_i - b_i)^2 w_{2,i}
over the class of pairs of vectors (a, b) such that a and b are isotonic (antitonic) and a_i ≤ b_i for all i = {1, …, n}. The estimates are computed with a projected subgradient algorithm where the projection is calculated using a suitable version of the pool-adjacent-violaters algorithm (PAVA).
The algorithm is implemented for antitonic curves in the function BoundedAntiMeanTwo
.
The function BoundedIsoMeanTwo
solves the same problem for isotonic curves, by simply invoking
BoundedAntiMeanTwo
and suitably flipping some of the arguments.
g1 |
The estimated function \hat g_1^\circ. |
g2 |
The estimated function \hat g_2^\circ. |
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. |
Fadoua Balabdaoui fadoua@ceremade.dauphine.fr
http://www.ceremade.dauphine.fr/~fadoua
Kaspar Rufibach (maintainer) kaspar.rufibach@gmail.com
http://www.kasparrufibach.ch
Filippo Santambrogio filippo.santambrogio@math.u-psud.fr
http://www.math.u-psud.fr/~santambr/
Balabdaoui, F., Rufibach, K., Santambrogio, F. (2009). Least squares estimation of two ordered monotone regression curves. Preprint.
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
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 | ## ========================================================
## 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
## --------------------------------------------------------
shor <- BoundedAntiMeanTwo(g1, w1, g2, w2, errorPrec = 20,
delta = 10^(-10))
## corresponding isotonic problem
shor2 <- BoundedIsoMeanTwo(-g2, w2, -g1, w1, errorPrec = 20,
delta = 10^(-10))
## the following vectors are equal
shor$g1 - -shor2$g2
shor$g2 - -shor2$g1
## --------------------------------------------------------
## for comparison, compute estimates via cyclical projection
## algorithm due to Dykstra (1983) (isotonic problem)
## --------------------------------------------------------
dykstra1 <- BoundedIsoMeanTwoDykstra(-g2, w2, -g1, w1,
delta = 10^(-10))
## the following vectors are equal
shor2$g1 - dykstra1$g1
shor2$g2 - dykstra1$g2
## --------------------------------------------------------
## Checking of solution
## --------------------------------------------------------
# This compares the first component of shor$g1 with a^*_1:
c(shor$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(shor$g1, shor$g2, g1, g2)), xlab = expression(x),
ylab = "measurements and estimates")
points(x, g2, col = 3)
lines(x, shor$g1 + 0.01, col = 2, type = 's', lwd = 2)
lines(x, shor$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)
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 via cyclical projection algorithm due to
## Dysktra (1983)
## --------------------------------------------------------
dykstra1 <- BoundedIsoMeanTwoDykstra(g1, w1, g2, w2,
delta = 10^-10, output = TRUE)
## --------------------------------------------------------
## compute smoothed versions
## --------------------------------------------------------
g1_mon <- dykstra1$g1
g2_mon <- dykstra1$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("upper 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.