Bivariate isotonic regression with respect to simple (increasing) linear ordering on both variables.
biviso(y, w = NULL, eps = NULL, eps2 = 1e-9, ncycle = 50000,
fatal = TRUE, warn = TRUE)A matrix of the same dimensions as y containing the
corresponding isotonic values. It has an attribute icycle
equal to the number of cycles required to achieve convergence of the algorithm.
The matrix of observations to be isotonized. It must of course have at least two rows and at least two columns.
A matrix of weights, greater than or equal to zero, of the same
dimension as y. If left NULL then w is
created as a matrix all of whose entries are equal to 1.
Convergence criterion. The algorithm is deemed to have converged
if each entry of the output matrix, after the completion of the
current iteration, does not differ by more than eps from
the corresponding entry of the matrix after the completion of the
previous iteration. If this argument is not supplied it defaults
to sqrt(.Machine$double.eps).
Criterion used to determine whether isotonicity is “violated”, whence whether (further) application of the “pool adjacent violators” procedure is required.
The maximum number of cycles of the iteration procedure. Must be
at least 2 (otherwise an error is given). If the procedure has not
converged after ncycle iterations then an error is given.
(See below.)
Logical scalar. Should the function stop if the subroutine
returns an error code other than 0 or 4? If fatal is
FALSE then output is returned by the function even if
there was a “serious” fault. One can set fatal=FALSE
to inspect the values of the objective matrix at various interim
stages prior to convergence. See Examples.
Logical scalar. Should a warning be produced if the subroutine
returns a value of ifault equal to 4 (or to any other
non-zero value when fatal has been set to FALSE)?
The subroutine comprising Algorithm AS 206 produces an error
code ifault with values from 0 to 6
The meaning of these codes is as follows:
0: No error.
1: Convergence was not attained in ncycle cycles.
2: At least one entry of w was negative.
3: Either nrow(y) or ncol(y) was less than 2.
4: A near-zero weight less than delta=0.00001 was
replaced by delta.
5: Convergence was not attained and a non-zero weight
was replaced by delta.
6: All entries of w were less than delta.
If ifault==4 a warning is given. All of the other non-zero
values of ifault result in an error being given.
Rolf Turner rolfturner@posteo.net
This function appears not to achieve exact isotonicity, at least not quite. For instance one can do:
set.seed(42)
u <- matrix(runif(400),20,20)
iu <- biviso(u)
any(apply(iu,2,is.unsorted))
and get TRUE. It turns out that columns 13, 14, and 16 of
iu have exceptions to isotonicity. E.g. six of the values
of diff(iu[,13]) are less than zero. However only one of
these is less than sqrt(.Machine$double.eps), and then only
“marginally” smaller.
So some of these negative values are “numerically different”
from zero, but not by much. The largest in magnitude in this
example, from column 16, is -2.217624e-08 --- which is
probably not of “practical importance”.
Note also that this example occurs in a very artificial context in which there is no actual isotonic structure underlying the data.
See the paper by Bril et al., (References) and the references cited therein for details.
Bril, Gordon; Dykstra, Richard; Pillers Carolyn, and Robertson, Tim ; Isotonic regression in two independent variables; Algorithm AS 206; JRSSC (Applied Statistics), vol. 33, no. 3, pp. 352-357, 1984.
pava() pava.sa() ufit()
x <- 1:20
y <- 1:10
xy <- outer(x,y,function(a,b){a+b+0.5*a*b}) + rnorm(200)
ixy <- biviso(xy)
set.seed(42)
u <- matrix(runif(400),20,20)
v <- biviso(u)
progress <- list()
for(n in 1:9) progress[[n]] <- biviso(u,ncycle=50*n,fatal=FALSE,warn=FALSE)
Run the code above in your browser using DataLab