home *** CD-ROM | disk | FTP | other *** search
/ Big Green CD 8 / BGCD_8_Dev.iso / NEXTSTEP / UNIX / Educational / R-0.49-MI / R-0.49-I / library / splines < prev    next >
Encoding:
Text File  |  1997-09-14  |  8.4 KB  |  299 lines
  1. "bs"<- function(x, df = NULL, knots = NULL, degree = 3, intercept = F, Boundary.knots
  2.      = range(x))
  3. {
  4.     nx <- names(x)
  5.     x <- as.vector(x)
  6.     nax <- is.na(x)
  7.     if(nas <- any(nax))
  8.         x <- x[!nax]
  9.     if(!missing(Boundary.knots)) {
  10.         Boundary.knots <- sort(Boundary.knots)
  11.         outside <- (ol <- x < Boundary.knots[1]) | (or <- x > 
  12.             Boundary.knots[2])
  13.     }
  14.     else outside <- rep(F, length = length(x))
  15.     if(!missing(df) && missing(knots)) {
  16.         nIknots <- df - (degree + 1) + (1 - intercept)
  17.         if(nIknots < 0) {
  18.             nIknots <- 0
  19.             warning(paste("df was too small; have used ", (degree + 
  20.                 1) - (1 - intercept)))
  21.         }
  22.         if(nIknots > 0) {
  23.             knots <- seq(from = 0, to = 1, length = nIknots + 2)[ - 
  24.                 c(1, nIknots + 2)]
  25.             knots <- quantile(x[!outside], knots)
  26.         }
  27.         else knots <- NULL
  28.     }
  29.     Aknots <- sort(c(rep(Boundary.knots, degree + 1), knots))
  30.     if(any(outside)) {
  31.         warning("Some x values beyond boundary knots may cause ill-conditioned bases"
  32.             )
  33.         derivs <- seq(from = 0, to = degree)    
  34.     # tricky way to do factorial
  35.         scalef <- c(1, exp(cumsum(log(seq(degree)))))
  36.         basis <- array(0, c(length(x), length(Aknots) - degree - 1))
  37.         if(any(ol)) {
  38.             k.pivot <- Boundary.knots[1]
  39.             xl <- cbind(1, outer(x[ol] - k.pivot, 1:degree, "^"))
  40.             tt <- spline.des(Aknots, rep(k.pivot, degree + 1), 
  41.                 degree + 1, derivs)$design
  42.             basis[ol,  ] <- xl %*% (tt/scalef)
  43.         }
  44.         if(any(or)) {
  45.             k.pivot <- Boundary.knots[2]
  46.             xr <- cbind(1, outer(x[or] - k.pivot, 1:degree, "^"))
  47.             tt <- spline.des(Aknots, rep(k.pivot, degree + 1), 
  48.                 degree + 1, derivs)$design
  49.             basis[or,  ] <- xr %*% (tt/scalef)
  50.         }
  51.         if(any(inside <- !outside))
  52.             basis[inside,  ] <- spline.des(Aknots, x[inside], 
  53.                 degree + 1)$design
  54.     }
  55.     else basis <- spline.des(Aknots, x, degree + 1)$design
  56.     if(!intercept)
  57.         basis <- basis[, -1]
  58.     n.col <- ncol(basis)
  59.     if(nas) {
  60.         nmat <- matrix(NA, length(nax), n.col)
  61.         nmat[!nax,  ] <- basis
  62.         basis <- nmat
  63.     }
  64.     dimnames(basis) <- list(nx, 1:n.col)
  65.     a <- list(degree = degree, knots = knots, Boundary.knots = 
  66.         Boundary.knots, intercept = intercept, class = c("bs", "basis")
  67.         )
  68.     attributes(basis) <- c(attributes(basis), a)
  69.     basis
  70. }
  71. "ns"<-function(x, df = NULL, knots = NULL, intercept = F, Boundary.knots = range(x))
  72. {
  73.     nx <- names(x)
  74.     x <- as.vector(x)
  75.     nax <- is.na(x)
  76.     if(nas <- any(nax))
  77.         x <- x[!nax]
  78.     if(!missing(Boundary.knots)) {
  79.         Boundary.knots <- sort(Boundary.knots)
  80.         outside <- (ol <- x < Boundary.knots[1]) | (or <- x > 
  81.             Boundary.knots[2])
  82.     }
  83.     else outside <- rep(F, length = length(x))
  84.     if(!missing(df) && missing(knots)) {
  85. # df = number(interior knots) + 1 + intercept
  86.         nIknots <- df - 1 - intercept
  87.         if(nIknots < 0) {
  88.             nIknots <- 0
  89.             warning(paste("df was too small; have used ", 1 + 
  90.                 intercept))
  91.         }
  92.         if(nIknots > 0) {
  93.             knots <- seq(from = 0, to = 1, length = nIknots + 2)[ - 
  94.                 c(1, nIknots + 2)]
  95.             knots <- quantile(x[!outside], knots)
  96.         }
  97.         else knots <- NULL
  98.     }
  99.     Aknots <- sort(c(rep(Boundary.knots, 4), knots))
  100.     if(any(outside)) {
  101.         basis <- array(0, c(length(x), nIknots + 4))
  102.         if(any(ol)) {
  103.             k.pivot <- Boundary.knots[1]
  104.             xl <- cbind(1, x[ol] - k.pivot)
  105.             tt <- spline.des(Aknots, rep(k.pivot, 2), 4, c(0, 1))$
  106.                 design
  107.             basis[ol,  ] <- xl %*% tt
  108.         }
  109.         if(any(or)) {
  110.             k.pivot <- Boundary.knots[2]
  111.             xr <- cbind(1, x[or] - k.pivot)
  112.             tt <- spline.des(Aknots, rep(k.pivot, 2), 4, c(0, 1))$
  113.                 design
  114.             basis[or,  ] <- xr %*% tt
  115.         }
  116.         if(any(inside <- !outside))
  117.             basis[inside,  ] <- spline.des(Aknots, x[inside], 4)$
  118.                 design
  119.     }
  120.     else basis <- spline.des(Aknots, x, 4)$design
  121.     const <- spline.des(Aknots, Boundary.knots, 4, c(2, 2))$design
  122.     if(!intercept) {
  123.         const <- const[, -1]
  124.         basis <- basis[, -1]
  125.     }
  126.     qr.const <- qr(t(const))
  127.     basis <- as.matrix((t(qr.qty(qr.const, t(basis))))[,  - (1:2)])
  128.     n.col <- ncol(basis)
  129.     if(nas) {
  130.         nmat <- matrix(NA, length(nax), n.col)
  131.         nmat[!nax,  ] <- basis
  132.         basis <- nmat
  133.     }
  134.     dimnames(basis) <- list(nx, 1:n.col)
  135.     a <- list(degree = 4, knots = knots, Boundary.knots = Boundary.knots, 
  136.         intercept = intercept, class = c("ns", "basis"))
  137.     attributes(basis) <- c(attributes(basis), a)
  138.     basis
  139. }
  140. "predict.bs"<-function(object, newx, ...)
  141. {
  142.     if(missing(newx))
  143.         return(object)
  144.     a <- c(list(x = newx), attributes(object)[c("degree", "knots", 
  145.         "Boundary.knots", "intercept")])
  146.     do.call("bs", a)
  147. }
  148. "predict.ns"<-function(object, newx, ...)
  149. {
  150.     if(missing(newx))
  151.         return(object)
  152.     a <- c(list(x = newx), attributes(object)[c("degree", "knots", 
  153.         "Boundary.knots", "intercept")])
  154.     do.call("ns", a)
  155. }
  156. "interp.spline" <- function(x, y, ord = 4)
  157. {
  158.     x <- as.numeric(x)
  159.     lenx <- length(x)
  160.     y <- as.numeric(y)
  161.     leny <- length(y)
  162.     if(leny!=lenx)
  163.         stop("Lengths of x and y must be equal")
  164.     ind <- order(x)
  165.     x <- x[ind]
  166.     y <- y[ind]
  167.     ordm1 <- ord - 1
  168.     knots <- c(x[1:ordm1] + x[1] - x[ord], x, x[lenx + (1:ordm1) - ordm1] +
  169.         x[lenx] - x[lenx - ordm1])
  170.     derivs <- c(2, integer(lenx), 2)
  171.     x <- c(x[1], x, x[lenx])
  172.     sys.mat <- spline.des(knots, x, ord = ord, derivs = derivs)$design
  173.     coeff <- solve(sys.mat, c(0, y, 0))
  174.     list(knots = knots, coeff = coeff, order = ord)
  175. }
  176. "periodic.spline" <- function(x, y, knots, period = 2 * pi, ord = 4)
  177. {
  178.     x <- as.vector(x, "double")
  179.     y <- as.vector(y, "double")
  180.     lenx <- length(x)
  181.     if(lenx!=length(y))
  182.         stop("Lengths of x and y must match")
  183.     ind <- order(x)
  184.     x <- x[ind]
  185.     y <- y[ind]
  186.     if(any((x[-1] - x[ - lenx]) <= 0))
  187.         stop("Values of x must be strictly increasing")
  188.     if(!missing(knots))
  189.         period <- knots[length(knots) + 1 - ord] - knots[1]
  190.     else knots <- c(x[(lenx - (ord - 2)):lenx] - period, x, x[1:ord] + 
  191.             period)
  192.     if((x[lenx] - x[1]) >= period)
  193.         stop("The range of x values exceeds one period")
  194.     coeff.mat <- spline.des(knots, x, ord)$design
  195.     sys.mat <- coeff.mat[, (1:lenx)]
  196.     sys.mat[, 1:(ord - 1)] <- sys.mat[, 1:(ord - 1)] + coeff.mat[, lenx +
  197.         (1:(ord - 1))]
  198.     qrstr <- qr(sys.mat)
  199.     coeff <- qr.coef(qrstr, y)
  200.     coeff <- c(coeff, coeff[1:(ord - 1)])
  201.     list(knots = knots, coeff = coeff, order = ord, period = period, m = 
  202.         coeff.mat)
  203. }
  204. "spline.value" <- function(bspstr, x = seq(knots[ord], knots[ncoeff + 1], len = lenx), deriv = 0,
  205.     lenx = length(x))
  206. {
  207. ##    library.dynam("splines", "splines.o")
  208.     if(is.null(bspstr$knots))
  209.         stop("bspstr must be a B-spline structure")
  210.     if(missing(x) && missing(lenx))
  211.         stop("either x or lenx must be given to spline.value")
  212.     knots <- bspstr$knots
  213.     ord <- bspstr$order
  214.     ncoeff <- length(bspstr$coeff)
  215.     x.orig <- x
  216.     if(!is.null(period <- bspstr$period)) {
  217.         ind <- x < knots[ord]
  218.         if(any(ind))
  219.             x[ind] <- x[ind] + period * (1 + (knots[ord] - x[ind]) %/%
  220.                 period)
  221.         ind <- x > knots[ncoeff + 1]
  222.         if(any(ind))
  223.             x[ind] <- x[ind] - period * (1 + (x[ind] - knots[ncoeff +
  224.                 1]) %/% period)
  225.     }
  226.     xind <- order(x)
  227.     x <- x[xind]
  228.     rind <- order(xind)
  229.     z <- .C("spline_value",
  230.         as.double(bspstr$knots),
  231.         as.double(bspstr$coeff),
  232.         as.integer(length(bspstr$coeff)),
  233.         as.integer(bspstr$order),
  234.         as.double(x),
  235.         as.integer(lenx),
  236.         as.integer(deriv),
  237.         y = double(lenx))
  238.     list(x = x.orig, y = z$y[rind])
  239. }
  240. "spline.des" <- function(knots, x, ord = 4, derivs = integer(nx))
  241. {
  242. ##    library.dynam("splines", "splines.o")
  243.     knots <- sort(as.vector(knots))
  244.     x <- as.vector(x)
  245.     nk <- length(knots)
  246.     nx <- length(x)
  247.     ind <- order(x)
  248.     sortx <- x[ind]
  249.     ind <- order(ind)
  250.     if(sortx[1] < knots[ord] || sortx[nx] > knots[nk + 1 - ord])
  251.         stop(paste("The x data must be in the range", knots[ord], "to",
  252.             knots[nk + 1 - ord]))
  253.     if(length(derivs)!=nx)
  254.         stop("length of derivs must match length of x")
  255.     ncoef <- nk - ord
  256.     temp <- .C("spline_basis",
  257.         as.double(knots),
  258.         as.integer(ncoef),
  259.         as.integer(ord),
  260.         as.double(sortx),
  261.         as.integer(derivs),
  262.         as.integer(nx),
  263.         design = array(0, c(ord, nx)),
  264.         offsets = integer(nx))
  265.     design <- array(0, c(nx, ncoef))
  266.     d.ind <- array(c(rep(1:nx, rep(ord, nx)), outer(1:ord, temp$offsets,
  267.         "+")), c(nx * ord, 2))
  268.     design[d.ind] <- temp$design
  269.     list(knots = knots, order = ord, derivs = derivs, design = design[ind,  ])
  270. }
  271. "linear.interp" <- function(x, y, x0 = seq(x[1], x[lenx], len = lenx0), lenx0 = length(x0))
  272. {
  273. ##    library.dynam("splines","splines.o")
  274.     x <- as.numeric(x)
  275.     y <- as.numeric(y)
  276.     lenx <- length(x)
  277.     if(length(y)!=lenx)
  278.         stop("lengths of x and y must be the same")
  279.     xind <- order(x)
  280.     x <- x[xind]
  281.     y <- y[xind]
  282.     if(missing(x0) && missing(lenx0))
  283.         stop("either x0 or lenx0 must be given")
  284.     x0 <- as.numeric(x0)
  285.     x0ind <- order(x0)
  286.     x0ord <- x0[x0ind]
  287.     nvals <- length(x0)
  288.     if(x0ord[1] < x[1] || x0ord[nvals] > x[lenx])
  289.         stop(paste("x0 values must be in the range", x[1], "to", x[
  290.             lenx]))
  291.     y <- .C("lin_interp",
  292.         x,
  293.         y,
  294.         x0,
  295.         y = double(nvals),
  296.         as.integer(nvals))$y
  297.     list(x = x0, y = y[order(x0ind)])
  298. }
  299.