Iso (version 0.0-21)

ufit: Unimodal isotonic regression.

Description

A "divide and conquer" algorithm is applied to calculate the isotonic regression of a set of data, for a unimodal order. If the mode of the unimodal order is not specified, then the optimal (in terms of minimizing the error sum of squares) unimodal fit is calculated.

Usage

ufit(y, lmode=NULL, imode=NULL, x=NULL, w=NULL, lc=TRUE, rc=TRUE,
        type=c("raw","stepfun","both"))

Value

If type=="raw" then the value is a list with components:

x

The argument x if this is specified, otherwise the default value.

y

The fitted values.

mode

The value of the location of the mode as determined by lmode or imode if one of these was specified. Otherwise it is the value of the location of the mode which was found to minimize the error sum of squares.

mse

The mean squared error.

If type=="both" then a component h which is the step function representation of the isotonic regression is added to the foregoing list.

If type=="stepfun" then only the step function representation h is returned.

Arguments

y

Vector of data whose isotonic regression is to be calculated.

lmode

Numeric scalar specifiing the location of the mode. It must be one of the values of x (see below) otherwise an error is thrown.

imode

Integer scalar specifying the index, amongst the values of x (see below) of the location of the mode. It must be one of the indices from 1 to n, where n is the length of y, otherwise an error is thrown.

It is an error to specify both lmode and imode.

Note that if neither lmode nor imode is specified then the function performs an exhaustive search among all possible mode locations for the optimal (in terms of minimising the error sum of squares) location.

x

A somewhat notional vector of \(x\) values corresponding to the data vector y; the value of the mode must be given, or will be determined in terms of these \(x\) values. Conceptually the model is y = m(x) + E, where m() is a unimodal function with mode at lmode, and where E is random "error". If x is not specified, it defaults to an equi-spaced sequence of length n on [0,1] (where n is the length of y).

w

Optional vector of weights to be used for calculating a weighted isotonic regression; if w is not specified then all weights are taken to equal 1.

lc

Logical scalar; should the isotonization be left continuous? If lc==FALSE then the value of the isotonization just before the mode is set to NA, which causes line plots to have a jump discontinuity at (just to the left of) the mode. The default is lc=TRUE.

rc

Logical scalar; should the isotonization be right continuous? If rc==FALSE then the value of the isotonization just after the mode is set to NA, which causes line plots to have a jump discontinuity at (just to the right of) the mode. The default is rc=TRUE.

type

String specifying the type of the output; see Value. May be abbreviated.

Author

Rolf Turner

Details

This function dynamically loads fortran subroutines "pava", "ufit" and "unimode" to do the actual work.

References

Mureika, R. A., Turner, T. R. and Wollan, P. C. (1992). An algorithm for unimodal isotonic regression, with application to locating a maximum. University of New Brunswick Department of Mathematics and Statistics Technical Report Number 92 -- 4.

Robertson, T., Wright, F. T. and Dykstra, R. L. (1988). Order Restricted Statistical Inference. Wiley, New York.

Shi, Ning-Zhong. (1988) A test of homogeneity for umbrella alternatives and tables of the level probabilities. Commun. Statist. --- Theory Meth. vol. 17, pp. 657 -- 670.

Turner, T. R., and Wollan, P. C. (1997) Locating a maximum using isotonic regression. Computational Statistics and Data Analysis vol. 25, pp. 305 -- 320.

See Also

pava() biviso()

Examples

y <- c(0,1,2,3,3,2)
f1 <- ufit(y,lmode=0.4) # The third entry of the default
                        # value of x = c(0.0,0.2,0.4,0.6,0.8,1.0).
f2 <- ufit(y,imode=3)   # Identical to f1.
f3 <- ufit(y,lmode=3,x=1:6)   # Effectively the same as f1 and f2.
                              # But is different  in appearance.
f4 <- ufit(y,imode=3,x=1:6)   # Identical to f3.

if (FALSE) {
    ufit(y,lmode=3)     # Throws an error.
    ufit(y,imode=7)     # Throws an error.
}

x <- c(0.00,0.34,0.67,1.00,1.34,1.67,2.00,2.50,3.00,3.50,4.00,4.50,
       5.00,5.50,6.00,8.00,12.00,16.00,24.00)
y <- c(0.0,61.9,183.3,173.7,250.6,238.1,292.6,293.8,268.0,285.9,258.8,
       297.4,217.3,226.4,170.1,74.2,59.8,4.1,6.1)
z <- ufit(y,x=x,type="b")
plot(x,y)
lines(z,col="red")
plot(z$h,do.points=FALSE,col.hor="blue",col.vert="blue",add=TRUE)
abline(v=z$mode,col="green",lty=2)