# Bivariate isotonic regression. Bivariate isotonic regression with respect to simple (increasing) linear ordering on both variables. 1.1 ```r biviso(y, w = NULL, eps = NULL, eps2 = 1e-9, ncycle = 50000, fatal = TRUE, warn = TRUE) ``` ## Arguments - `y`: The matrix of observations to be isotonized. It must of course have at least two rows and at least two columns. - `w`: 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`. - `eps`: 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)`. - `eps2`: Criterion used to determine whether isotonicity is violated , whence whether (further) application of the pool adjacent violators procedure is required. - `ncycle`: 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.) - `fatal`: 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 . - `warn`: 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`)? ## Error Messages 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. ## Details See the paper by Bril et al., (**References**) and the references cited therein for details. ## Returns 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. ## References 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. ## Author(s) Rolf Turner rolfturner@posteo.net ## WARNING 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 Also `pava()` `pava.sa()` `ufit()` ## Examples ```r 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) ```

biviso function