#!/bin/sh # This is a shell archive (produced by GNU sharutils 4.2). # To extract the files from this archive, save it to some FILE, remove # everything before the `!/bin/sh' line above, then type `sh FILE'. # # Made on 1998-04-25 07:06 PDT by . # Source directory was `/f7/trevor/MDA'. # # Existing files will *not* be overwritten unless `-c' is specified. # # This shar contains: # length mode name # ------ ---------- ------------------------------------------ # 592 -rw-r--r-- README # 48056 -rw-r--r-- dumpdata.mda # 60886 -rw-r--r-- helpfile.shar # 81669 -rw-r--r-- ratfor.shar # 1000 -rw-r--r-- makefile # save_IFS="${IFS}" IFS="${IFS}:" gettext_dir=FAILED locale_dir=FAILED first_param="$1" for dir in $PATH do if test "$gettext_dir" = FAILED && test -f $dir/gettext \ && ($dir/gettext --version >/dev/null 2>&1) then set `$dir/gettext --version 2>&1` if test "$3" = GNU then gettext_dir=$dir fi fi if test "$locale_dir" = FAILED && test -f $dir/shar \ && ($dir/shar --print-text-domain-dir >/dev/null 2>&1) then locale_dir=`$dir/shar --print-text-domain-dir` fi done IFS="$save_IFS" if test "$locale_dir" = FAILED || test "$gettext_dir" = FAILED then echo=echo else TEXTDOMAINDIR=$locale_dir export TEXTDOMAINDIR TEXTDOMAIN=sharutils export TEXTDOMAIN echo="$gettext_dir/gettext -s" fi touch -am 1231235999 $$.touch >/dev/null 2>&1 if test ! -f 1231235999 && test -f $$.touch; then shar_touch=touch else shar_touch=: echo $echo 'WARNING: not restoring timestamps. Consider getting and' $echo "installing GNU \`touch', distributed in GNU File Utilities..." echo fi rm -f 1231235999 $$.touch # if mkdir _sh20801; then $echo 'x -' 'creating lock directory' else $echo 'failed to create lock directory' exit 1 fi # ============= README ============== if test -f 'README' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'README' '(file already exists)' else $echo 'x -' extracting 'README' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'README' && #To install the mda software: #------------------------- mkdir .Data mkdir .Data/.Help sh ratfor.shar make cd .Data/.Help sh ../../helpfile.shar cd ../../ # pick one of the following #S < dumpdata.mda Splus < dumpdata.mda # -------------------------- # In S(plus) change the objects ".bruto.object" and ".mars.object" to point # to the appropriate places. # # Trevor Hastie and Rob Tibshirani. # # Send bug reports to: # Trevor Hastie trevor@playfair.stanford.edu # Phone: 415-725-2231 Fax: 415-725-8977 # Statistics Department, Sequoia Hall, Stanford University, Ca94305 SHAR_EOF $shar_touch -am 0601190695 'README' && chmod 0644 'README' || $echo 'restore of' 'README' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'README:' 'MD5 check failed' 6ec0db289d9222043e9ea1071a183938 README SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'README'`" test 592 -eq "$shar_count" || $echo 'README:' 'original size' '592,' 'current size' "$shar_count!" fi fi # ============= dumpdata.mda ============== if test -f 'dumpdata.mda' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'dumpdata.mda' '(file already exists)' else $echo 'x -' extracting 'dumpdata.mda' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'dumpdata.mda' && ".First.lib"<- function(lib.loc, section) { X path <- paste(lib.loc, section, "mda.so", sep = "/") X dyn.load.shared(path) X library(classif) } ".bruto.object"<- "/fs/trevor/docs/discr/BRUTO.o" ".mars.object"<- "/n/rice/usr/trevor/docs/discr/MARS.o" ".mda.object"<- "/usr/trevor/docs/mda/mda.o" "bruto"<- function(x, y, w = rep(1, n), wp = rep(1/np, np), dfmax, cost = 2, maxit.select X = 20, maxit.backfit = 20, thresh = 0.0001, trace.bruto = F, X start.linear = T, fit.object, ...) { X this.call <- match.call() X y <- as.matrix(y) X x <- as.matrix(x) X np <- ncol(y) X d <- dim(x) X nq <- d[2] X n <- d[1] X xnames <- dimnames(x)[[2]] X if(!length(xnames)) X xnames <- NULL X ynames <- dimnames(y)[[2]] X if(!length(ynames)) X ynames <- NULL X storage.mode(x) <- "double" X storage.mode(y) <- "double" X storage.mode(w) <- "double" X storage.mode(wp) <- "double" X storage.mode(cost) <- "double" X if(missing(fit.object)) { X nknotl <- function(n) X { X a1 <- log(50)/log(2) X a2 <- log(100)/log(2) X a3 <- log(140)/log(2) X a4 <- log(200)/log(2) # Cutoff Criteria X cx <- as.vector(cut(n, c(0, 50, 200, 800, 3200))) X if(is.na(cx)) X cx <- 5 X floor(switch(cx, X n, X 2^(a1 + ((a2 - a1) * (n - 50))/150), X 2^(a2 + ((a3 - a2) * (n - 200))/600), X 2^(a3 + ((a4 - a3) * (n - 800))/2400), X 200 + (n - 3200)^0.20000000000000001) + 6) X } X check.range <- apply(x, 2, var) X if(any(check.range < .Machine$single.eps)) X stop("A column of x is constant; do not include an intercept column" X ) X nkmax <- nknotl(n) - 4 X coef <- matrix(double(nkmax * np * nq), ncol = nq) X ybar <- apply(y * w, 2, sum)/sum(w) X if(start.linear && (nq * cost > n)) X start.linear <- F X if(start.linear) { X start.fit <- polyreg(x, y, w) X eta <- fitted(start.fit) X coef[seq(from = 2, by = 2, length = np), ] <- t( X start.fit$coef)[, -1] X type <- as.integer(rep(2, nq)) X df <- as.double(rep(1, nq)) X } X else { X eta <- outer(rep(1, n), ybar) X type <- integer(nq) X df <- double(nq) X } X nk <- integer(nq) X knot <- matrix(double((nkmax + 4) * nq), ncol = nq) X Match <- matrix(integer(n * nq), ncol = nq) X nef <- integer(nq) X lambda <- double(nq) X xrange <- matrix(double(2 * nq), 2, nq) X df <- double(nq) X if(missing(dfmax)) X dfmax <- (2 * nkmax)/3 X if(length(dfmax) != nq) X dfmax <- rep(dfmax, length = nq) X if(cost > 0) { X TD <- (n - sum(df))/cost X TT <- dfmax > TD X if(any(TT)) X dfmax[TT] <- TD X } X storage.mode(dfmax) <- "double" X } X else { X this.call <- fit.object$call X ybar <- fit.object$ybar X nkmax <- fit.object$nkmax X dfmax <- fit.object$dfmax X eta <- fit.object$fitted.values X if(is.null(eta)) X eta <- predict(fit.object, x) X nk <- fit.object$nk X knot <- fit.object$knot X Match <- fit.object$Match X nef <- fit.object$nef X lambda <- fit.object$lambda X coef <- fit.object$coef X type <- codes(fit.object$type) X xrange <- fit.object$xrange X maxit.select <- 0 X maxit.backfit <- fit.object$nit[2] X df <- fit.object$df X } X maxit <- as.integer(c(maxit.select, maxit.backfit)) X names(df) <- xnames X names(maxit) <- c("selection", "backfitting") X gcv.select <- if(maxit.select) matrix(double(maxit.select * nq), nq, X maxit.select) else double(1) X gcv.backfit <- if(maxit.backfit) matrix(double(maxit.backfit * nq), nq, X maxit.backfit) else double(1) X df.select <- if(maxit.select) matrix(double(maxit.select * nq), nq, X maxit.select) else double(1) X names(lambda) <- xnames X if(!is.loaded(Fortran.symbol("bruto"))) X dyn.load(.bruto.object) X fit <- .Fortran("bruto", X x, X as.integer(n), X as.integer(nq), X y, X as.integer(np), X w, X knot = knot, X nkmax = as.integer(nkmax), X nk = nk, X wp, X Match = Match, X nef = nef, X dfmax = dfmax, X cost = cost, X lambda = lambda, X df = df, X coef = coef, X type = type, X xrange = xrange, X gcv.select = gcv.select, X gcv.backfit = gcv.backfit, X df.select = df.select, X maxit = maxit, X nit = maxit, X fitted.values = eta, X residuals = y - eta, X as.double(thresh), X double((2 * np + 2) * ((n + 1) + 1) + (2 * np + 16) * (n + 1) + X 2 * (n + 1) + np), X trace.bruto)[c("knot", "nkmax", "nk", "Match", "nef", "dfmax", X "cost", "lambda", "df", "coef", "type", "xrange", "gcv.select", X "gcv.backfit", "df.select", "maxit", "nit", "fitted.values", X "residuals")] X if(TN <- fit$nit[1]) { X TT <- fit$gcv.select[, seq(TN), drop = F] X dimnames(TT) <- list(xnames, NULL) X } X else TT <- NULL X fit$gcv.select <- TT X if(TN) { X TT <- fit$df.select[, seq(TN), drop = F] X dimnames(TT) <- list(xnames, NULL) X } X else TT <- NULL X fit$df.select <- TT X if(TN <- fit$nit[2]) { X TT <- fit$gcv.backfit[, seq(TN), drop = F] X dimnames(TT) <- list(xnames, NULL) X } X else TT <- NULL X fit$gcv.backfit <- TT X TT <- factor(fit$type, levels = 1:3, labels = c("excluded", "linear", X "smooth")) X names(TT) <- xnames X fit$type <- TT X fit$ybar <- ybar X fit$call <- this.call X structure(fit, class = "bruto") } "coef.fda"<- function(object, type = c("canonical", "discriminant"), ...) { X type <- match.arg(type) X fit <- object$fit X Coefs <- fit$coef X if(is.null(Coefs)) X stop("No explicit coefficients in this formulation") X dimension <- object$dimension X Coefs <- Coefs %*% object$theta[, seq(dimension), drop = F] X lambda <- object$values X alpha <- sqrt(lambda[seq(dimension)]) X sqima <- sqrt(1 - lambda[seq(dimension)]) X Coefs <- scale(Coefs, F, sqima * alpha) X if(type == "discriminant") X Coefs <- Coefs %*% t(object$means) X Coefs } "confusion"<- function(object, ...) UseMethod("confusion") "confusion.default"<- function(predict, true, ...) { X if(inherits(predict, "data.frame")) X confusion.list(predict, true) X else { X jt <- table(predict, true) X jd <- dimnames(jt) X jn <- unlist(jd) X ju <- jn[duplicated(jn)] X j1 <- jd[[1]][!match(jd[[1]], ju, 0)] X j2 <- jd[[2]][!match(jd[[2]], ju, 0)] X jt <- jt[c(ju, j1), c(ju, j2), drop = F] X realjt <- jt[ju, ju, drop = F] X ntot <- sum(jt) X mismatch <- (ntot - sum(realjt))/ntot X structure(jt, error = (1 - sum(diag(realjt))/sum(realjt)), X mismatch = if(mismatch > 0) mismatch else NULL) X } } "confusion.fda"<- function(object, data, ...) { X if(missing(data)) X return(object$confusion) X Terms <- terms(object) X attr(Terms, "intercept") <- 0 X m <- model.frame(Terms, data) X x <- model.matrix(Terms, m) X g <- model.extract(m, response) X confusion.default(predict(object, x, ...), g) } "confusion.list"<- function(pred, truth) { X dd <- names(pred) X y <- seq(dd) X x <- attr(pred, "dimension") X if(!length(x)) X x <- seq(dd) X for(i in y) { X confi <- confusion(pred[, i], truth) X y[i] <- attr(confi, "error") X } X return(x, y) } "contr.fda"<- function(p = rep(1, d[1]), contrast.default = contr.helmert(length(p))) { X d <- dim(contrast.default) X sqp <- sqrt(p/sum(p)) X x <- cbind(1, contrast.default) * outer(sqp, rep(1, d[2] + 1)) X qx <- qr(x) X J <- qx$rank X qr.qy(qx, diag(d[1])[, seq(2, J)])/outer(sqp, rep(1, J - 1)) } "df.inv"<- function(d, df, lambda = 10^4, iterations = 10) { X df.diriv <- function(d, lambda) X - sum((d * lambda)/(d + lambda)^2) X current.df <- sum(d/(d + lambda)) X # cat("lambda,df", lambda, current.df, "\n") X if(abs((df - current.df)/df) < 0.0001 | iterations == 1) X return(lambda, df = current.df) X else { X lambda <- exp(log(lambda) - (current.df - df)/df.diriv(d, X lambda)) X Recall(d, df, lambda, iterations - 1) X } } "dumpdata.mda"<- c("bruto", "coef.fda", "confusion", "contr.fda", "fda", "mars", X "model.matrix.mars", "polybasis", "polyreg", "pplot", "predict.bruto", X "predict.fda", "predict.mars", "predict.polyreg", "print.fda", X "softmax", ".bruto.object", ".mars.object", "confusion.default", X "confusion.fda", "confusion.list", "fix.theta", "gen.ridge", X "kmeans.start", "llmult", "lvq.start", "mda", "mda.fit", "mda.means", X "mda.start", "mean.penalty", "plot.fda", "pplot4", "predict.mda", X "print.mda", "shrink", "shrink.mda", "transform.penalty") "fda"<- function(formula = formula(data), data = sys.parent(), weights, theta, X dimension = J - 1, eps = .Machine$single.eps, method = polyreg, X keep.fitted = (n * dimension < 1000), ...) { # # Flexible Discriminant Analysis # Function for fitting models described in # Hastie, Tibshirani and Buja, 1994, JASA # "Flexible Discriminant Analysis by Optimal Scoring" # Hastie, Buja and Tibshirani, 1995, Annals of Statistics # "Penalized Discriminant Analysis" # Modified 2/15/95 by T Hastie # X this.call <- match.call() X # This extracts the x and g from the formula or data frame # # -------< not for human consumption >-------- X m <- match.call(expand = F) X m[[1]] <- as.name("model.frame") X m <- m[match(names(m), c("", "formula", "data", "weights"), 0) > 0] X m <- eval(m, sys.parent()) X Terms <- attr(m, "terms") X g <- model.extract(m, response) X attr(Terms, "intercept") <- 0 X x <- model.matrix(Terms, m) X dd <- dim(x) X n <- dd[1] X weights <- model.extract(m, weights) X if(!length(weights)) X weights <- rep(1, n) X else if(any(weights < 0)) X stop("negative weights not allowed") X if(length(g) != n) X stop("g should have length nrow(x)") X fg <- factor(g) # # if some levels are missing, this gets rid of them # X prior <- table(fg) X prior <- prior/sum(prior) X cnames <- levels(fg) X g <- as.numeric(fg) X J <- length(cnames) # # # construct indicator matrix for response # # Initialization uses orthogonal contrasts for theta # Even if contrasts are supplied, they need to be orthogonalized X iswt <- F X if(missing(weights)) X dp <- table(g)/n X else { X weights <- (n * weights)/sum(weights) X dp <- tapply(weights, g, sum)/n X iswt <- T X } X if(missing(theta)) X theta <- contr.helmert(J) X theta <- contr.fda(dp, theta) X Theta <- theta[g, , drop = F] # #Theta is now an n x K matrix of responses for the (nonparametric) # regression, normalized wrt the data (i.e. theta normalized wrt dp) # where K is min(J-1, ncol of starting theta) # X fit <- method(x, Theta, weights, ...) X if(iswt) X Theta <- Theta * weights X ssm <- t(Theta) %*% fitted(fit)/n # # This n could be (n-1) to make it "unbiassed" # # now we need the svd of ssm (unweighted) # X ed <- svd(ssm, nu = 0) X thetan <- ed$v X lambda <- ed$d X # Note: the discriminant eigenvalues are a transformation # of the optimal scaling values X lambda[lambda > 1 - eps] <- 1 - eps # # If lambda is one we get errors, so we illiminate this possibility X discr.eigen <- lambda/(1 - lambda) X pe <- (100 * cumsum(discr.eigen))/sum(discr.eigen) X dimension <- min(dimension, sum(lambda > eps)) X if(dimension == 0) { X warning("degenerate problem; no discrimination") X return(structure(list(dimension = 0, fit = fit, call = X this.call), class = "fda")) X } X thetan <- thetan[, seq(dimension), drop = F] # X pe <- pe[seq(dimension)] #Now produce projected centroids # X alpha <- sqrt(lambda[seq(dimension)]) X sqima <- sqrt(1 - lambda[seq(dimension)]) # # package up results # X vnames <- paste("v", seq(dimension), sep = "") X means <- scale(theta %*% thetan, F, sqima/alpha) X dimnames(means) <- list(cnames, vnames) X names(lambda) <- c(vnames, rep("", length(lambda) - dimension)) # X names(pe) <- vnames X obj <- structure(list(percent.explained = pe, values = lambda, means = X means, theta.mod = thetan, dimension = dimension, prior = prior, X fit = fit, call = this.call, terms = Terms), class = "fda") X obj$confusion <- confusion(predict(obj), fg) X # get rid of the fitted values; these take up too much space X if(!keep.fitted) obj$fit$fitted.values <- NULL # X obj } "fix.theta"<- function(theta, Q) { X M <- t(theta) %*% Q %*% theta X eM <- eigen(M, sym = T) X scale(theta %*% eM$vectors, F, sqrt(eM$values)) } "gen.ridge"<- function(x, y, weights, lambda = 1, omega, df, ...) { # # ||Y-XB||^2 + lambda*B^T Omega Beta # = ||Y-X* B*||^2 + lambda ||B*||^2 # where X*=XU*, Omega=UDU^T, U*=UD^{-1/2} and B*=D^{1/2}U^T B # This is a simple ridge problem now # Let X* = UDV^T # then H = UDV^T(V^TDV + lambda I)^{-1} VDU^T # = UD(D^2 + lambda I)^{-1}DU^T # with trace sum(d_j^2/(lambda+d_j^2) # The coefficients B* are given by V(D^2 + lambda I)^{-1}DU^TY # and then we see that U*B*=B # X if(missing(df) && lambda <= .Machine$single.eps) return(polyreg(x, y)) X d <- dim(x) X mm <- apply(x, 2, mean) X x <- scale(x, mm, F) X simple <- if(missing(omega)) T else F X if(!simple) { X if(!all(match(c("values", "vectors"), names(omega), F))) X stop("You must supply an eigen-decomposed version of omega" X ) X vals <- pmax(.Machine$single.eps, sqrt(omega$values)) X basis <- scale(omega$vectors, F, vals) X x <- x %*% basis X } X svd.x <- svd(x) X dd <- svd.x$d X if(!missing(df)) { X while(sum(dd^2/(dd^2 + lambda)) > df) lambda <- lambda * 10 X junk <- df.inv(dd^2, df, lambda) X lambda <- junk$lambda X df <- junk$df X } X else df <- sum(dd^2/(dd^2 + lambda)) X y <- (t(t(y) %*% svd.x$u) * dd)/(dd^2 + lambda) X coef <- svd.x$v %*% y X fitted <- x %*% coef X if(!simple) X coef <- basis %*% coef X structure(list(fitted.values = fitted, coefficients = coef, df = df, X lambda = lambda, xmeans = mm), class = "gen.ridge") } "kmeans.start"<- function(x, g, subclasses) { X cnames <- levels(g <- factor(g)) X J <- length(cnames) X g <- as.numeric(g) X weights <- as.list(cnames) X names(weights) <- cnames X subclasses <- rep(subclasses, length = length(cnames)) # X R <- sum(subclasses) X cl <- rep(seq(J), subclasses) X cx <- x[seq(R), , drop = F] X if(!is.loaded(symbol.For("kmns"))) X dyn.load("mda.o") X for(j in seq(J)) { X nc <- subclasses[j] X which <- cl == j X if(nc <= 1) { X cx[which, ] <- apply(x[g == j, , drop = F], 2, mean) X wmj <- matrix(1, sum(g == j), 1) X } X else { X xx <- x[g == j, , drop = F] X start <- xx[sample(1:nrow(xx), size = nc), ] X TT <- kmeans(xx, start) X cx[which, ] <- TT$centers X wmj <- diag(nc)[TT$cluster, ] X } X dimnames(wmj) <- list(NULL, paste("s", seq(nc), sep = "")) X weights[[j]] <- wmj X } X list(x = cx, cl = factor(cl, labels = cnames), weights = weights) } "laplacian"<- function(size = 16, compose = F) { #build gamma #build gamma # Here I follow very closely the material on page 635 in JASA 1991 # of O'Sullivan's article on discretized Laplacian Smoothing X gmat <- matrix(0, size, size) X xx <- seq(size) X for(v in xx) X gmat[, v] <- sqrt(2/size) * cos(((xx - 0.5) * pi * (v - 1))/ X size) X gmat[, 1] <- gmat[, 1]/sqrt(2) X lvec <- - (2 * size^2) * (1 - cos(((xx - 1) * pi)/size)) X gmat <- kronecker(gmat, gmat) X lvec <- rep(lvec, rep(size, size)) + rep(lvec, size) X if(compose) X gmat %*% (lvec^2 * t(gmat)) X else list(vectors = gmat, values = lvec^2) } "llmult"<- function(p, g) { X index <- cbind(seq(along = g), as.numeric(g)) X p <- p[index] X -2 * sum(log(p[p > .Machine$single.eps])) } "lvq.start"<- function(x, g, subclasses) { X cnames <- levels(fg <- factor(g)) X J <- length(cnames) X g <- as.numeric(g) X weights <- as.list(cnames) X names(weights) <- cnames X subclasses <- rep(subclasses, length = length(cnames)) # X size <- sum(subclasses) X if(!is.loaded(symbol.For("olvq"))) X dyn.load("mda.o") X cb <- lvqinit(x, g, size = size) X TT <- olvq1(x, g, codebk = cb) X TT <- lvq3(x, g, codebk = TT) X cl <- as.numeric(TT$cl) X R <- length(cl) X cx <- TT$x X p <- ncol(cx) #Compute the weights based on within class assignments # X for(j in seq(J)) { X which <- cl == j X number <- sum(which) X if(number == 0) { X cx <- rbind(cx, apply(x[g == j, ], 2, mean)) X cl <- c(cl, j) X wmj <- matrix(1, sum(g == j), 1) X number <- 1 X } X else if(number == 1) X wmj <- matrix(1, sum(g == j), 1) X else { X jcx <- cx[which, ] X jcl <- seq(number) X jcluster <- lvqtest(list(x = jcx, cl = jcl), x[g == j, X ]) X needed <- unique(jcluster) X rcl <- rep(0, number) X rcl[needed] <- j X cl[which] <- rcl X wmj <- diag(number)[jcluster, needed, drop = F] X number <- length(needed) X } X dimnames(wmj) <- list(NULL, paste("s", seq(number), sep = "")) X weights[[j]] <- wmj X } X TT <- cl > 0 X list(x = cx[TT, , drop = F], cl = factor(cl[TT], labels = cnames), X weights = weights) } "make.dumpdata.mda"<- function() dump(dump.list, "dumpdata.mda") "mars"<- function(x, y, w = rep(1, nrow(x)), wp, degree = 1, nk = max(21, 2 * ncol(x) + X 1), penalty = 2, thresh = 0.001, prune = T, trace.mars = F, X forward.step = T, prevfit = NULL, ...) { X this.call <- match.call() X if((nk %% 2) != 1) X nk <- nk - 1 X x <- as.matrix(x) X np <- dim(x) X n <- np[1] X p <- np[2] X y <- as.matrix(y) X nclass <- ncol(y) X if(is.null(np)) { X np <- c(length(x), 1) X x <- as.matrix(x) X } X if(forward.step) { X interms <- 1 X lenb <- nk X bx <- matrix(rep(0, nrow(x) * nk), nrow = n) X res <- matrix(rep(0, nrow(x) * ncol(y)), nrow = n) X fullin <- rep(0, nk) X cuts <- NULL X factor <- NULL X } X else { X bx <- model.matrix.mars(prevfit, x, full = T) X interms <- ncol(bx) X lenb <- prevfit$lenb X o <- prevfit$all.terms X fullin <- rep(0, ncol(bx)) X fullin[o] <- 1 X res <- prevfit$res X factor <- prevfit$factor X cuts <- prevfit$cuts X if(missing(penalty)) X penalty <- prevfit$penalty X degree <- prevfit$degree X nk <- lenb X thresh <- prevfit$thresh X } X if(missing(penalty) & (degree > 1)) X penalty <- 3 X if(!missing(wp)) { X if(any(wp <= 0)) X stop("wp should all be positive") X wp <- sqrt(wp/sum(wp)) X y <- y * outer(rep(1, n), wp) X } X else wp <- NULL X tagx <- x X storage.mode(tagx) <- "integer" X for(j in 1:p) { X tagx[, j] <- order(x[, j]) X } X bestin <- rep(0, nk) X flag <- matrix(rep(0, nk * p), nrow = nk, ncol = p) X if(is.null(cuts)) X cuts <- matrix(rep(0, nk * p), nrow = nk, ncol = p) X if(is.null(factor)) { X dir <- matrix(rep(0, nk * p), nrow = nk, ncol = p) X } X else { X dir <- factor X } X alpha <- rep(0, nclass) X beta <- matrix(rep(0, nk * nclass), nrow = nk) X bestgcv <- 0 X storage.mode(y) <- "double" X storage.mode(x) <- "double" X storage.mode(bx) <- "double" X storage.mode(flag) <- "integer" X storage.mode(cuts) <- "double" X storage.mode(dir) <- "double" X storage.mode(res) <- "double" X storage.mode(beta) <- "double" X lenscrat <- 1 + n + 2 * n * nk + 4 * nk * nk + 3 * nk + 3 * nk * nclass + X 3 * nclass + 28 * n + 51 # X if(!is.loaded(Fortran.symbol("marss"))) X dyn.load(.mars.object) X junk <- .Fortran("marss", X as.integer(n), X as.integer(n), X as.integer(p), X as.integer(nclass), X as.matrix(y), X as.matrix(x), X as.double(w), X as.matrix(tagx), X as.integer(degree), X as.integer(nk), X as.double(penalty), X as.double(thresh), X as.logical(forward.step), X as.integer(interms), X as.logical(prune), X bx = as.matrix(bx), X fullin = as.integer(fullin), X lenb = as.integer(lenb), X bestgcv = as.double(bestgcv), X bestin = as.integer(bestin), X flag = as.matrix(flag), X cuts = as.matrix(cuts), X dir = as.matrix(dir), X res = as.matrix(res), X alpha = as.double(alpha), X beta = as.matrix(beta), X double(lenscrat), X integer(4 * nk), X trace.mars) X lenb <- junk$lenb X all.terms <- seq(lenb)[junk$fullin[1:lenb] == 1] X selected.terms <- seq(lenb)[junk$bestin[1:lenb] == 1] X coefficients <- junk$beta[seq(selected.terms), , drop = F] X residuals <- junk$res X fitted.values <- y - residuals X if(!is.null(wp)) { X TT <- outer(rep(1, n), wp) X residuals <- residuals/TT X fitted.values <- fitted.values/TT X coefficients <- coefficients/outer(rep(1, length(selected.terms X )), wp) X } X dir <- junk$dir[seq(lenb), , drop = F] X dimnames(dir) <- list(NULL, dimnames(x)[[2]]) X cutss <- junk$cuts[seq(lenb), , drop = F] X x <- junk$bx[, selected.terms, drop = F] X structure(list(call = this.call, all.terms = all.terms, selected.terms X = selected.terms, penalty = penalty, degree = degree, nk = nk, X thresh = thresh, gcv = junk$bestgcv, factor = dir, cuts = cutss, X residuals = residuals, fitted.values = fitted.values, lenb = X junk$lenb, coefficients = coefficients, x = x), class = "mars") } "mda"<- function(formula = formula(data), data = sys.parent(), subclasses = 3, sub.df X = NULL, tot.df = NULL, dimension = sum(subclasses) - 1, eps = .Machine$ X single.eps, iter = 5, weights = mda.start(x, g, subclasses, trace, ...), X method = polyreg, keep.fitted = (n * dimension < 1000), trace = F, ... X ) { # # Mixture Discriminant Analysis # Function for fitting models described in # Hastie and Tibshirani 1995 JRSSB # modified by Trevor on 2/15/95 # X this.call <- match.call() # # This extracts the x and g from the formula or data frame # # -------< not for human consumption >-------- X m <- match.call(expand = F) X m[[1]] <- as.name("model.frame") X m <- m[match(names(m), c("", "formula", "data"), 0) > 0] X m <- eval(m, sys.parent()) X Terms <- attr(m, "terms") X g <- model.extract(m, response) X attr(Terms, "intercept") <- 0 X x <- model.matrix(Terms, m) X dd <- dim(x) X n <- dd[1] X m <- dd[2] # # # Define a function needed later X rowmin <- function(mat) X { X ncc <- ncol(mat) X if(ncc == 1) X return(drop(mat)) X rowm <- pmin(mat[, 1], mat[, 2]) X if(ncc == 2) X return(rowm) X else { X for(j in seq(from = 3, to = ncc)) X rowm <- pmin(rowm, mat[, j]) X } X rowm X } # ---------------------------------------------- # X if(length(g) != n) stop("g should have length nrow(x)") # # turn g into a factor # if some levels are missing, this gets rid of them X fg <- factor(g) # # weights is a special beast. It is a list of matrices of probabilities # that are the subclass probabilites. The names correspond to the levels # of the (implicit) factor g # # One can supply an mda object itself as weights # X if(inherits(weights, "mda")) { X if(is.null(weights$weights)) X weights <- predict(weights, x, type = "weights", g = fg X ) X else weights <- weights$weights X } X subclasses <- sapply(weights, ncol) # # extract codes and class labels # I don't use codes since the class labels can get mixed up # X prior <- table(fg) X dim(prior) <- NULL X prior <- prior/sum(prior) X cnames <- levels(fg) X g <- as.numeric(fg) # X J <- length(cnames) # X Assign <- split(seq(sum(subclasses)), rep(seq(J), subclasses)) X names(Assign) <- cnames # # # see if shrinking was called for # X if(!is.null(tot.df)) { X if(tot.df >= sum(subclasses)) X tot.df <- NULL X } X if(!is.null(sub.df)) { X sub.df <- rep(sub.df, length = length(prior)) X sub.df <- pmin(sub.df, subclasses) X if(all(sub.df == subclasses)) X sub.df <- NULL X } # generate starting sub-class weight X for(counter in seq(iter)) { X fit <- mda.fit(x, g, weights, assign.theta = Assign, Rj = X subclasses, sub.df = sub.df, tot.df = tot.df, dimension X = dimension, eps = .Machine$single.eps, method = X method, trace = trace, ...) # # predict.fda works on a mda.fit object # X dmat <- predict.fda(fit, type = "distance") X sub.prior <- fit$sub.prior # # Now recompute weights # X for(j in seq(J)) { X TT <- dmat[g == j, Assign[[j]], drop = F] X TT <- exp(-0.5 * (TT - rowmin(TT))) X TT <- TT * outer(rep(1, nrow(TT)), sub.prior[[j]]) X weights[[j]][] <- TT/drop(TT %*% rep(1, ncol(TT))) X } # #Compute posterior probabilities # X pclass <- matrix(1, n, J) X dmat <- exp(-0.5 * (dmat - rowmin(dmat))) X for(j in seq(J)) { X priorj <- sub.prior[[j]] X ass <- Assign[[j]] X TT <- dmat[, ass, drop = F] * outer(rep(1, n), priorj) X TTot <- drop(TT %*% rep(1, length(ass))) X pclass[, j] <- prior[j] * TTot X } X pclass <- pclass/drop(pclass %*% rep(1, J)) X if(trace) X cat("Iteration", counter, "\tDeviance(multinomial)", X format(round(ll <- llmult(pclass, g), 5)), "\n" X ) X } X if(!trace) ll <- llmult(pclass, g) X # get rid of the fitted values; these take up too much space X if(!keep.fitted) fit$fit$fitted.values <- NULL # X dimnames(pclass) <- list(NULL, names(Assign)) X conf <- confusion(softmax(pclass), fg) X fit <- c(fit, list(weights = weights, prior = prior, assign.theta = X Assign, deviance = ll, confusion = conf, terms = Terms)) X fit$call <- this.call X fit$sub.df <- sub.df X fit$tot.df <- tot.df X class(fit) <- c("mda", "fda") X fit } "mda.fit"<- function(x, g, weights, theta, assign.theta, Rj, sub.df, tot.df, dimension = R - X 1, eps = .Machine$single.eps, method = polyreg, ...) { X this.call <- match.call() X n <- nrow(x) X cnames <- names(weights) X J <- length(cnames) # # now extract the lengths and names of the appropriate matrices in # weights X R <- sum(Rj) # # Get the total weight for each subclass (to be used in weighting the scores) # Keep it in list form # X wtots <- lapply(weights, function(x) X apply(x, 2, sum)) # X sub.prior <- lapply(wtots, function(x) X x/sum(x)) # #dp is the unconditional probability of being in sublcass r of class j # X dp <- unlist(wtots) X subclass.names <- names(dp) X dp <- dp/sum(dp) # # Initialization uses orthogonal contrasts for theta # Even if contrasts are supplied, they need to be orthogonalized X if(missing(theta)) X theta <- contr.helmert(R) X theta <- contr.fda(dp, theta) # # Here is the fix for shrinking # # X if(!(is.null(sub.df) & is.null(tot.df))) { X Q <- diag(dp) + transform.penalty(prior = dp, cl = rep(seq(J), X Rj), df = sub.df, tot.df = tot.df) X theta <- fix.theta(theta, Q) X } # # Now compute the weighted score means for each observation # X Theta <- matrix(0, n, R - 1) X obs.weights <- double(n) X for(j in seq(J)) { X Theta[g == j, ] <- weights[[j]] %*% theta[assign.theta[[j]], , X drop = F] X obs.weights[g == j] <- weights[[j]] %*% rep(1, Rj[j]) X } #Theta is now an n x K matrix of responses for the regression, normalized wrt # the weights (i.e. theta normalized wrt dp) # where K is min(R-1, ncol of starting theta) # # X fit <- method(x, Theta, obs.weights, ...) X Theta <- Theta * obs.weights X ssm <- t(Theta) %*% fitted(fit)/n # # now we need the svd of ssm (unweighted) # X ed <- svd(ssm, nu = 0) X thetan <- ed$v X lambda <- ed$d X # Note: the discriminant eigenvalues are a transformation # of the optimal scaling values X lambda[lambda > 1 - eps] <- 1 - eps # # If lambda is one we get errors, so we illiminate this possibility X discr.eigen <- lambda/(1 - lambda) X pe <- (100 * cumsum(discr.eigen))/sum(discr.eigen) X dimension <- min(dimension, sum(lambda > eps)) X if(dimension == 0) { X warning("degenerate problem; no discrimination") X return(structure(list(dimension = 0, fit = fit, call = X this.call), class = "fda")) X } X thetan <- thetan[, seq(dimension), drop = F] # X pe <- pe[seq(dimension)] #Now produce projected centroids # X alpha <- sqrt(lambda[seq(dimension)]) X sqima <- sqrt(1 - lambda[seq(dimension)]) # # package up results # X vnames <- paste("v", seq(dimension), sep = "") X means <- scale(theta %*% thetan, F, sqima/alpha) X dimnames(means) <- list(subclass.names, vnames) X names(lambda) <- c(vnames, rep("", length(lambda) - dimension)) # X names(pe) <- vnames X structure(list(percent.explained = pe, values = lambda, means = means, X theta.mod = thetan, dimension = dimension, sub.prior = X sub.prior, fit = fit, call = this.call)) } "mda.means"<- function(object, x, y) { X weights <- means <- object$weights X nn <- names(object$weights) X for(i in nn) { X xx <- x[y == i, ] X ww <- weights[[i]] X nc <- ncol(ww) X xm <- matrix(0, ncol(x), nc) X for(j in seq(nc)) { X www <- ww[, j] X www <- www/sum(www) X xm[, j] <- apply(xx * www, 2, sum) X } X means[[i]] <- xm X } X means } "mda.start"<- function(x, g, subclasses = 3, trace.mda.start = F, start.method = c("kmeans", X "lvq"), tries = 5, criterion = c("misclassification", "deviance"), ... X ) { X if(!length(find("lvqtest"))) X stop("mda requires functions in the classif collection donated by Brian Ripley to statlib/S" X ) X start.method <- match.arg(start.method) X if((start.method == "kmeans") && !length(find("kmeans"))) X stop("mda with start.method=kmeans requires the kmeans() function, currently only available in Splus" X ) X criterion <- match.arg(criterion) X name.criterion <- switch(criterion, X misclassification = "Misclassification Error", X deviance = "Deviance(multinomial)") X starter <- get(paste(start.method, "start", sep = "."), mode = X "function") X fg <- factor(g) X cnames <- levels(fg) X prior <- table(fg) X J <- length(cnames) X n <- length(g) X g <- as.numeric(fg) # Now loop over tries # X best.ll <- 1/.Machine$single.eps X for(try in seq(tries)) { X start <- starter(x, fg, subclasses) X weights <- start$weights X if(criterion == "misclassification") { X pg <- lvqtest(start, x) X ll <- attr(confusion(pg, fg), "error") X } X else { X subclasses <- sapply(weights, ncol) # X Assign <- split(seq(sum(subclasses)), rep(seq(J), X subclasses)) X names(Assign) <- cnames # X fit <- mda.fit(x, g, weights, assign.theta = Assign, Rj X = subclasses, eps = .Machine$single.eps, X method = polyreg, ...) # # predict.fda works on a mda.fit object # X dmat <- exp(-0.5 * predict.fda(fit, type = "distance")) X sub.prior <- fit$sub.prior # # #Now compute the weights and the posteriors # X pclass <- matrix(1, n, J) X for(j in seq(J)) { X priorj <- sub.prior[[j]] X ass <- Assign[[j]] X TT <- dmat[, ass, drop = F] X TT <- TT * outer(rep(1, n), priorj) X TTot <- drop(TT %*% rep(1, length(ass))) X wmj <- TT[g == j, , drop = F]/TTot[g == j] X pclass[, j] <- prior[j] * TTot X dimnames(wmj) <- list(NULL, paste("s", seq( X along = ass), sep = "")) X weights[[j]] <- wmj X } X pclass <- pclass/drop(pclass %*% rep(1, J)) X ll <- llmult(pclass, g) X } X if(trace.mda.start) X cat(start.method, "start \t", name.criterion, format( X round(ll, 5)), "\n") X if(ll < best.ll) { X keep.weights <- weights X best.ll <- ll X } X } X structure(keep.weights, criterion = best.ll, name.criterion = X name.criterion) } "mean.penalty"<- function(prior, cl) { X n <- length(prior) X Q <- matrix(0, n, n) X ll <- unique(cl) X for(lll in ll) { X which <- cl == lll X pp <- prior[which] X npp <- length(pp) X pp <- diag(npp) - outer(rep(1, npp), pp/sum(pp)) X Q[which, which] <- t(pp) %*% pp X } X attr(Q, "cl") <- cl X attr(Q, "prior") <- prior X Q } "model.matrix.mars"<- function(object, x, which = object$selected.terms, full = F, ...) { # # make mars design matrix from output of mars # x can be any matrix of predictors (without a column of 1s) # # if x is missing, the x component of the mars object is returned # which terms to compute; if not given, the component selected.terms # from the mars object is used. which is a vector of indexes, ranging # from 1 thru nrow(object$factor) # # if full =T, a full matrix (including 0 columns # or unused terms) is returned X if(missing(x)) return(object$x) X x <- as.matrix(x) X dd <- dim(x) X n <- dd[1] X p <- dd[2] X nterms <- length(which) X dir <- object$factor X cut <- object$cuts X if(full) { X bx <- matrix(0, nrow = n, ncol = object$lenb) X bx[, 1] <- 1 X } X else bx <- matrix(1, nrow = n, ncol = nterms) X which <- which[-1] X for(i in seq(along = which)) { X j <- which[i] X if(all(dir[j, ] == 0)) { X stop("error in factor or which") X } X temp1 <- 1 X for(k in 1:p) { X if(dir[j, k] != 0) { X temp2 <- dir[j, k] * (x[, k] - cut[j, k]) X temp1 <- temp1 * temp2 * (temp2 > 0) X } X } X if(full) X bx[, j] <- temp1 X else bx[, i + 1] <- temp1 X } X bx } "plot.fda"<- function(object, data, g, coords = c(1, 2), group = c("true", "predicted"), X colors, pch, mcolors = max(colors) + 1, mpch, pcex = 0.5, mcex = 2.5, X ...) { X group <- match.arg(group) X if(missing(data)) { X vars <- predict(object, type = "var") X g <- predict(object) X group <- "predict" X } X else { X ff <- terms(object) X attr(ff, "intercept") <- 0 X m <- model.frame(ff, data) X x <- model.matrix(ff, m) X vars <- predict(object, x, type = "var") X if(group == "true") X g <- model.extract(m, response) X else g <- predict(object, x) X } X means <- object$means X g <- as.factor(g) X cc <- as.numeric(g) X np <- seq(levels(g)) X if(missing(colors)) X colors <- np X else colors <- rep(colors, length = length(np)) X if(missing(pch)) X pch <- paste(np) X else pch <- rep(paste(pch), length = length(np)) X mcolors <- rep(mcolors, length = length(np)) X if(missing(mpch)) X mpch <- pch X else mpch <- rep(paste(mpch), length = length(np)) X assign <- object$assign X if(is.null(assign)) X assign <- split(seq(np), seq(np)) X if(!is.matrix(coords)) { X coords <- matrix(coords, length(coords), length(coords)) X tt <- lower.tri(coords) X coords <- cbind(t(coords)[tt], coords[tt]) X } X for(ii in seq(nrow(coords))) { X coord.pair <- coords[ii, ] X plot(rbind(vars[, coord.pair], means[, coord.pair]), ..., type X = "n", xlab = paste("Discriminant Var", coord.pair[1]), X ylab = paste("Discriminant Var", coord.pair[2]), main X = paste("Discriminant Plot for", group, "classes")) X for(i in np) { X which <- cc == i X if(any(which)) X points(vars[which, coord.pair, drop = F], col X = colors[i], pch = pch[i], cex = pcex) X points(means[assign[[i]], coord.pair, drop = F], col = X mcolors[i], pch = 1, cex = mcex) X points(means[assign[[i]], coord.pair, drop = F], col = X mcolors[i], pch = mpch[i], cex = mcex/2) X } X } X invisible() } "polybasis"<- function(x, degree = 1, monomial = F) { X if(degree >= 3) X warning("This is not a smart polynomial routine. You may get numerical problems with a degree of 3 or more" X ) X x <- as.matrix(x) X dn <- dimnames(x) X dd <- dim(x) X np <- dd[2] X if(np == 1) X monomial <- T X if(degree > 1) { X if(monomial) { X px <- x X cc <- sapply(split(paste(diag(np)), rep(seq(np), rep(np, X np))), paste, collapse = "") X tx <- x X for(i in 2:degree) { X px <- px * tx X x <- cbind(x, px) X cc <- c(cc, sapply(split(paste(diag(np) * i), X rep(seq(np), rep(np, np))), paste, collapse X = "")) X } X } X else { X matarray <- array(x, c(dd, degree)) X for(i in 2:degree) X matarray[, , i] <- x^i X matarray <- aperm(matarray, c(1, 3, 2)) X x <- matarray[, , np] X ad0 <- seq(degree) X ad <- ad0 X ncol.mat0 <- degree X ncol.x <- degree X d0 <- paste(ad0) X cc <- d0 X for(ii in seq(np - 1, 1)) { X index0 <- rep(seq(ncol.mat0), ncol.x) X index <- rep(seq(ncol.x), rep(ncol.mat0, ncol.x X )) X newad <- ad0[index0] + ad[index] X retain <- newad <= degree X mat0 <- matarray[, , ii] X if(any(retain)) X newmat <- mat0[, index0[retain], drop = F] * X x[, index[retain], drop = F] X else newmat <- NULL X ddn <- paste(d0[index0[retain]], cc[index[ X retain]], sep = "") X zeros <- paste(rep(0, nchar(cc[1])), collapse X = "") X cc <- paste(0, cc, sep = "") X d00 <- paste(d0, zeros, sep = "") X x <- cbind(mat0, x, newmat) X cc <- c(d00, cc, ddn) X ad <- c(ad0, ad, newad[retain]) X ncol.x <- length(ad) X } X } X if(!is.null(dn)) X dn[[2]] <- cc X else dn <- list(NULL, cc) X dimnames(x) <- dn X } X cbind(Intercept = 1, x) } "polyreg"<- function(x, y, w, degree = 1, monomial = F, ...) { X x <- polybasis(x, degree, monomial) X if(iswt <- !missing(w)) { X if(any(w <= 0)) X stop("only positive weights") X w <- sqrt(w) X y <- y * w X x <- x * w X } X qrx <- qr(x) X coef <- as.matrix(qr.coef(qrx, y)) X fitted <- qr.fitted(qrx, y) X if((df <- qrx$rank) < ncol(x)) X coef[qrx$pivot, ] <- coef X if(iswt) X fitted <- fitted/w X structure(list(fitted.values = fitted, coefficients = coef, degree = X degree, monomial = monomial, df = df), class = "polyreg") } "pplot"<- function(x, g, colors, pch, add = F, type = "p", ...) { X g <- as.factor(g) X cc <- codes(g) X np <- seq(levels(g)) X if(missing(colors)) X colors <- np + 1 X else colors <- rep(colors, length = length(np)) X if(missing(pch)) X pch <- paste(np) X else pch <- rep(pch, length = length(np)) X if(!add) X plot(x, type = "n", ...) X for(i in unique(cc)) X points(x[cc == i, , drop = F], col = colors[i], pch = pch[i], X type = type) } "pplot4"<- function(x, ...) { X par(mfrow = c(3, 2)) X for(i in 1:3) { X for(j in (i + 1):4) X pplot(x[, c(i, j)], xlab = paste("var", i), ylab = X paste("var", j), ...) X } } "predict.bruto"<- function(object, x, type = c("fitted", "terms")) { X if(missing(x)) { X z <- fitted(object) X if(is.null(z)) X stop("need to supply x") X else return(z) X } X d <- as.integer(dim(x)) X type <- match.arg(type) X nq <- d[2] X n <- d[1] X if(nq != length(object$df)) X stop("x should have the same number of columns as the df component of object" X ) X ybar <- object$ybar X np <- as.integer(length(ybar)) X eta <- matrix(double(n * np), n, np) X Type <- codes(object$type) X storage.mode(Type) <- "integer" X storage.mode(x) <- "double" X if(!is.loaded(Fortran.symbol("bruto"))) X dyn.load(.bruto.object) X if(type == "fitted") { X .Fortran("pbruto", X x, X n, X nq, X ybar, X np, X object$knot, X object$nkmax, X object$nk, X object$coef, X Type, X object$xrange, X eta = eta, X eta)$eta X } X else { X ob <- as.list(seq(nq)) X names(ob) <- dimnames(x)[[2]] X knot <- object$knot X nk <- object$nk X xrange <- object$xrange X coef <- object$coef X fitm <- matrix(double(n * np), n, np) X dimnames(fitm) <- list(dimnames(x)[[1]], names(ybar)) X for(i in seq(nq)) { X if(Type[i] > 1) X fit <- .Fortran("psspl2", X x[, i], X n, X np, X knot[, i], X nk[i], X xrange[, i], X coef[, i], X coef[, i], X fit = fitm, X as.integer(0), X Type[i])$fit X else fit <- fitm X ob[[i]] <- list(x = x[, i], y = fit) X } X ob X } } "predict.fda"<- function(object, x, type = c("class", "variates", "posterior", "hierarchical", X "distances"), prior, dimension = J - 1) { X dist <- function(x, mean, m = ncol(mean)) X (scale(x, mean, F)^2) %*% rep(1, m) X type <- match.arg(type) X means <- object$means X Jk <- dim(means) X J <- Jk[1] X k <- Jk[2] # # Note for type=="hierarchical" dimension can be a vector # X if(type == "hierarchical") { X if(missing(dimension)) X dimension.set <- seq(k) X else { X dimension.set <- dimension[dimension <= k] X if(!length(dimension.set)) X dimension.set <- k X dimension <- max(dimension.set) X } X } X else dimension <- min(max(dimension), k) X if(missing(x)) X y <- predict(object$fit) X else { X if(inherits(x, "data.frame") || is.list(x)) { X Terms <- delete.response(terms(object)) X attr(Terms, "intercept") <- 0 X x <- model.matrix(Terms, x) X } X y <- predict(object$fit, x) X } X y <- y %*% object$theta[, seq(dimension), drop = F] X lambda <- object$values X alpha <- sqrt(lambda[seq(dimension)]) X sqima <- sqrt(1 - lambda[seq(dimension)]) X x <- scale(y, F, sqima * alpha) X if(missing(prior)) X prior <- object$prior X else { X if(any(prior < 0) | round(sum(prior), 5) != 1) X stop("innappropriate prior") X } X means <- means[, seq(dimension), drop = F] X switch(type, X variates = return(x), X class = { X n <- nrow(x) X prior <- 2 * log(prior) X mindist <- dist(x, means[1, ], dimension) - prior[1] X pclass <- rep(1, n) X for(i in seq(2, J)) { X ndist <- dist(x, means[i, ], dimension) - X prior[i] X l <- ndist < mindist X pclass[l] <- i X mindist[l] <- ndist[l] X } X levels(pclass) <- dimnames(means)[[1]] X return(factor(pclass)) X } X , X posterior = { X pclass <- matrix(0, nrow(x), J) X for(i in seq(J)) X pclass[, i] <- exp(-0.5 * dist(x, means[i, ], X dimension)) * prior[i] X dimnames(pclass) <- list(dimnames(x)[[1]], dimnames( X means)[[1]]) X return(pclass/drop(pclass %*% rep(1, J))) X } X , X hierarchical = { X prior <- 2 * log(prior) X Pclass <- vector("list", length(dimension.set)) X names(Pclass) <- paste("D", dimension.set, sep = "") X for(ad in seq(along = dimension.set)) { X d <- dimension.set[ad] X dd <- seq(d) X mindist <- dist(x[, dd, drop = F], means[1, dd, X drop = F], d) - prior[1] X pclass <- rep(1, nrow(x)) X for(i in seq(2, J)) { X ndist <- dist(x[, dd, drop = F], means[i, dd, X drop = F], d) - prior[i] X l <- ndist < mindist X pclass[l] <- i X mindist[l] <- ndist[l] X } X levels(pclass) <- dimnames(means)[[1]] X Pclass[[ad]] <- pclass X } X rownames <- dimnames(x)[[1]] X if(is.null(rownames)) X rownames <- paste(seq(nrow(x))) X return(structure(Pclass, class = "data.frame", X row.names = rownames, dimensions = X dimension.set)) X } X , X distances = { X dclass <- matrix(0, nrow(x), J) X for(i in seq(J)) X dclass[, i] <- dist(x, means[i, ], dimension) X dimnames(dclass) <- list(dimnames(x)[[1]], dimnames( X means)[[1]]) X return(dclass) X } X ) } "predict.gen.ridge"<- function(object, x, ...) { X if(missing(x)) X fitted(object) X else scale(x, object$xmeans, F) %*% object$coef } "predict.mars"<- function(object, x) { # # computes fitted values for design points x, based on mars fit # in object # X if(missing(x)) { X z <- fitted(object) X if(is.null(z)) X stop("need to supply x") X else return(z) X } X model.matrix.mars(object, x) %*% object$coefficients } "predict.mda"<- function(object, x, type = c("class", "variates", "posterior", "hierarchical", X "weights"), prior = NULL, dimension = R - 1, g, ...) { X type <- match.arg(type) # # Note for type=="hierarchical" dimension can be a vector # X Rk <- dim(object$means) X R <- Rk[1] X k <- Rk[2] X if(type == "hierarchical") { X if(missing(dimension)) X dimension.set <- seq(k) X else { X dimension.set <- dimension[dimension <= k] X if(!length(dimension.set)) X dimension.set <- k X dimension <- max(dimension.set) X } X Pclass <- vector("list", length(dimension.set)) X names(Pclass) <- paste("D", dimension.set, sep = "") X for(ad in seq(along = dimension.set)) { X d <- dimension.set[ad] X Pclass[[ad]] <- if(missing(x)) Recall(object, prior = X prior, dimension = d, ...) else Recall(object, X x, prior = prior, dimension = d, ...) X } X rownames <- names(Pclass[[1]]) X if(is.null(rownames)) X rownames <- paste(seq(along = Pclass[[1]])) X return(structure(Pclass, class = "data.frame", row.names = X rownames, dimensions = dimension.set)) X } X else dimension <- min(max(dimension), k) X if(is.null(prior)) X prior <- object$prior X else { X if(any(prior < 0) | round(sum(prior), 5) != 1) X stop("innappropriate prior") X } X if(type == "variates") return(NextMethod("predict")) X # Define a function needed later X rowmin <- function(mat) X { X ncc <- ncol(mat) X if(ncc == 1) X return(drop(mat)) X rowm <- pmin(mat[, 1], mat[, 2]) X if(ncc == 2) X return(rowm) X else { X for(j in seq(from = 3, to = ncc)) X rowm <- pmin(rowm, mat[, j]) X } X rowm X } X dmat <- if(missing(x)) predict.fda(object, type = "distances", X dimension = dimension, ...) else predict.fda(object, x, X type = "distances", dimension = dimension, ...) X Assign <- object$assign X sub.prior <- object$sub.prior X J <- length(Assign) X if(type == "weights") { X if(missing(x)) X return(object$weights) X g <- as.numeric(g) X weights <- Assign X for(j in seq(J)) { X TT <- dmat[g == j, Assign[[j]], drop = F] X TT <- exp(-0.5 * (TT - rowmin(TT))) X TT <- TT * outer(rep(1, nrow(TT)), sub.prior[[j]]) X weights[[j]] <- TT/drop(TT %*% rep(1, ncol(TT))) X } X return(weights) X } X pclass <- matrix(1, nrow(dmat), J) X dmat <- exp(-0.5 * (dmat - rowmin(dmat))) X for(j in seq(J)) { X TT <- dmat[, Assign[[j]], drop = F] X TT <- TT * outer(rep(1, nrow(TT)), sub.prior[[j]]) X pclass[, j] <- prior[j] * drop(TT %*% rep(1, ncol(TT))) X } # dimnames(pclass) <- list(dimnames(x)[[1]], names(Assign)) X dimnames(pclass) <- list(NULL, names(Assign)) X switch(type, X class = softmax(pclass), X posterior = pclass/drop(pclass %*% rep(1, J))) } "predict.polyreg"<- function(object, x, ...) { X if(missing(x)) { X z <- fitted(object) X if(is.null(z)) X stop("need to supply x") X else return(z) X } X degree <- object$degree X monomial <- object$monomial X polybasis(x, degree, monomial) %*% object$coef } "print.fda"<- function(x, ...) { X if(!is.null(cl <- x$call)) { X cat("Call:\n") X dput(cl) X } X cat("\nDimension:", format(x$dimension), "\n") X cat("\nPercent Between-Group Variance Explained:\n") X print(round(x$percent, 2)) X error <- x$confusion X df <- x$fit X if(!is.null(df)) X df <- df$df X if(!is.null(df)) { X cat("\nDegrees of Freedom (per dimension):", format(sum(df)), X "\n") X } X if(!is.null(error)) { X n <- as.integer(sum(error)) X error <- format(round(attr(error, "error"), 5)) X cat("\nTraining Misclassification Error:", error, "( N =", n, X ")\n") X } X invisible(x) } "print.mda"<- function(x, ...) { X NextMethod("print") X if(!is.null(x$deviance)) X cat("\nDeviance:", format(round(x$deviance, 3)), "\n") X invisible(x) } "shrink"<- function(object, ...) UseMethod("shrink") "shrink.mda"<- function(object, sub.df = NULL, tot.df = NULL, ...) { # # This function takes an MDA object and shrinks # it, creating a new mda object # # First check for a null condition # X if(is.null(sub.df) & is.null(tot.df)) { X warning("No shrinkage parameters (sub.df or tot.df) given") X return(object) X } # # First recover theta X dimension <- object$dimension X lambda <- object$values[seq(dimension)] X theta.mod <- object$theta.mod X theta <- object$means X alpha <- sqrt(lambda) X sqima <- sqrt(1 - lambda) # X theta <- scale(theta, F, alpha/sqima) # # Construct penalty # X sub.prior <- unlist(object$sub.prior) X prior <- object$prior X Rj <- sapply(object$assign.theta, length) X dp <- sub.prior * rep(prior, Rj) X cl <- rep(seq(Rj), Rj) X P <- diag(dp) + transform.penalty(prior = dp, cl = cl, df = sub.df, X tot.df = tot.df) X K <- t(theta) %*% P %*% theta X TT <- chol((K + t(K))/2) X Tinv <- backsolve(TT, diag(length(lambda))) # # lambda*Tinv = diag(lambda)%*% Tinv # X M <- t(Tinv) %*% (lambda * Tinv) X ed <- svd(M) X thetan <- ed$v X lambda <- ed$d X # Note: the discriminant eigenvalues are a transformation # of the optimal scaling values X discr.eigen <- lambda/(1 - lambda) X pe <- (100 * cumsum(discr.eigen))/sum(discr.eigen) X dimension <- min(dimension, sum(lambda > .Machine$single.eps)) X if(dimension == 0) { X warning("degenerate problem; no discrimination") X return(structure(list(dimension = 0, fit = fit, call = X this.call), class = "fda")) X } X thetan <- thetan[, seq(dimension), drop = F] # X pe <- pe[seq(dimension)] #Now produce projected centroids # X alpha <- sqrt(lambda[seq(dimension)]) X sqima <- sqrt(1 - lambda[seq(dimension)]) # # package up results # X dm <- dimnames(theta) X vnames <- dm[[2]][seq(dimension)] X means <- scale(theta %*% Tinv %*% thetan, F, sqima/alpha) X dimnames(means) <- list(dm[[1]], vnames) X names(lambda) <- c(vnames, rep("", length(lambda) - dimension)) # X names(pe) <- vnames X theta.mod <- theta.mod %*% Tinv %*% thetan X object$confusion <- object$deviance <- NULL X incl.names <- c("percent.explained", "values", "means", "theta.mod", X "dimension") X rl <- list(pe, lambda, means, theta.mod, dimension) X names(rl) <- incl.names X object$sub.df <- sub.df X object$tot.df <- tot.df X object[incl.names] <- rl X object$weights <- NULL X update(object, sub.df = sub.df, tot.df = tot.df, weights = object, ...) } "softmax"<- function(x, gap = F) { X d <- dim(x) X maxdist <- x[, 1] X pclass <- rep(1, d[1]) X for(i in seq(2, d[2])) { X l <- x[, i] > maxdist X pclass[l] <- i X maxdist[l] <- x[l, i] X } X dd <- dimnames(x)[[2]] X if(gap) { X x <- abs(maxdist - x) X x[cbind(seq(d[1]), pclass)] <- drop(x %*% rep(1, d[2])) X gaps <- do.call("pmin", data.frame(x)) X } X pclass <- if(is.null(dd) || !length(dd)) pclass else factor(pclass, X levels = seq(d[2]), labels = dd) X if(gap) X list(class = pclass, gaps = gaps) X else pclass } "transform.penalty"<- function(Q, prior, cl, df = NULL, tot.df = NULL) { X if(missing(Q)) X Q <- mean.penalty(prior, cl) X if(missing(prior)) X prior <- attr(Q, "prior") X if(missing(cl)) X cl <- attr(Q, "cl") X transform.pen <- function(Q, prior, df) X { X df.inv <- function(d, df, lambda = NULL, iterations = 10) X { # #this function solves for lambda such that sum(1/(1 + d*lambda)) = df X if(is.null(lambda)) { X lambda <- 0.10000000000000001 X while(sum(1/(1 + d * lambda)) >= df) lambda <- X lambda * 2 X } X df.diriv <- function(d, lambda) X - sum((d * lambda)/(1 + d * lambda)^2) X current.df <- sum(1/(1 + d * lambda)) X if(abs((df - current.df)/df) < 0.0001 | iterations == 1 X ) X return(lambda, df = current.df) X else { #cat(df, lambda, current.df,"\n") X lambda <- exp(log(lambda) - (current.df - df)/ X df.diriv(d, lambda)) X Recall(d, df, lambda, iterations - 1) X } X } X pQp <- Q/outer(sqrt(prior), sqrt(prior)) X d <- svd(pQp)$d X lambda <- df.inv(d, df)$lambda X lambda * Q X } X if(!is.null(tot.df)) { X if(tot.df >= length(prior)) X return(Q * 0) X else return(transform.pen(Q, prior, tot.df)) X } X else { X ncl <- unique(cl) X df <- rep(df, length = length(ncl)) X for(i in seq(along = ncl)) { X which <- cl == ncl[i] X Q[which, which] <- Recall(Q[which, which, drop = F], X prior[which], tot.df = df[i]) X } X return(Q) X } } SHAR_EOF $shar_touch -am 0425070698 'dumpdata.mda' && chmod 0644 'dumpdata.mda' || $echo 'restore of' 'dumpdata.mda' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'dumpdata.mda:' 'MD5 check failed' 210a5cc2f4b4b4aac752897c3a31af1b dumpdata.mda SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'dumpdata.mda'`" test 48056 -eq "$shar_count" || $echo 'dumpdata.mda:' 'original size' '48056,' 'current size' "$shar_count!" fi fi # ============= helpfile.shar ============== if test -f 'helpfile.shar' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'helpfile.shar' '(file already exists)' else $echo 'x -' extracting 'helpfile.shar' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'helpfile.shar' && #!/bin/sh # This is a shell archive (produced by GNU sharutils 4.2). # To extract the files from this archive, save it to some FILE, remove # everything before the `!/bin/sh' line above, then type `sh FILE'. # # Made on 1998-04-25 07:06 PDT by . # Source directory was `/f7/trevor/MDA'. # # Existing files will *not* be overwritten unless `-c' is specified. # # This shar contains: # length mode name # ------ ---------- ------------------------------------------ # 2695 -rw-r--r-- .Data/.Help/bruto # 643 -rw-r--r-- .Data/.Help/confusion # 4562 -rw-r--r-- .Data/.Help/fda # 814 -rw-r--r-- .Data/.Help/gen.ridge # 4438 -rw-r--r-- .Data/.Help/junk # 4436 -rw-r--r-- .Data/.Help/junk~ # 2983 -rw-r--r-- .Data/.Help/mars # 2916 -rw-r--r-- .Data/.Help/mars~ # 6491 -rw-r--r-- .Data/.Help/mda # 681 -rw-r--r-- .Data/.Help/mda.start # 593 -rw-r--r-- .Data/.Help/model.matrix.mars # 451 -rw-r--r-- .Data/.Help/polyreg # 610 -rw-r--r-- .Data/.Help/predict.bruto # 1525 -rw-r--r-- .Data/.Help/predict.fda # 237 -rw-r--r-- .Data/.Help/predict.mars # 1975 -rw-r--r-- .Data/.Help/predict.mda # 590 -rw-r--r-- .Data/.Help/softmax # save_IFS="${IFS}" IFS="${IFS}:" gettext_dir=FAILED locale_dir=FAILED first_param="$1" for dir in $PATH do X if test "$gettext_dir" = FAILED && test -f $dir/gettext \ X && ($dir/gettext --version >/dev/null 2>&1) X then X set `$dir/gettext --version 2>&1` X if test "$3" = GNU X then X gettext_dir=$dir X fi X fi X if test "$locale_dir" = FAILED && test -f $dir/shar \ X && ($dir/shar --print-text-domain-dir >/dev/null 2>&1) X then X locale_dir=`$dir/shar --print-text-domain-dir` X fi done IFS="$save_IFS" if test "$locale_dir" = FAILED || test "$gettext_dir" = FAILED then X echo=echo else X TEXTDOMAINDIR=$locale_dir X export TEXTDOMAINDIR X TEXTDOMAIN=sharutils X export TEXTDOMAIN X echo="$gettext_dir/gettext -s" fi touch -am 1231235999 $$.touch >/dev/null 2>&1 if test ! -f 1231235999 && test -f $$.touch; then X shar_touch=touch else X shar_touch=: X echo X $echo 'WARNING: not restoring timestamps. Consider getting and' X $echo "installing GNU \`touch', distributed in GNU File Utilities..." X echo fi rm -f 1231235999 $$.touch # if mkdir _sh20746; then X $echo 'x -' 'creating lock directory' else X $echo 'failed to create lock directory' X exit 1 fi # ============= .Data/.Help/bruto ============== if test ! -d '.Data'; then X $echo 'x -' 'creating directory' '.Data' X mkdir '.Data' fi if test ! -d '.Data/.Help'; then X $echo 'x -' 'creating directory' '.Data/.Help' X mkdir '.Data/.Help' fi if test -f '.Data/.Help/bruto' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/bruto' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/bruto' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/bruto' && XX.BG XX.FN bruto XX.TL fit an additive spline model by adaptive backfitting XX.CS bruto(x, y, w, wp, dfmax, cost, maxit.select, maxit.backfit, XX thresh=0.0001, trace=T, start.linear=T, fit.object, ...) XX.PP XX.AG x a matrix of numeric predictors (does not include the column of 1s) XX.AG y a vector or matrix of responses XX.AG w optional observation weight vector XX.AG wp optional weight vector for each column of y; the RSS and GCV criteria use a weighted sum of squared residuals. XX.AG dfmax a vector of maximum df (degrees of freedom) for each term XX.AG cost cost per degree of freedom; default is 2. XX.AG maxit.select maximum number of iterations during the selection stage XX.AG maxit.backfit maximum number of iterations for the final backfit stage (with fixed lambda) XX.AG thresh convergence threshold (default is 0.0001); iterations cease when the relative change in GCV is below this threshold XX.AG trace logical flag. If TRUE (default) a progress report is printed during the fitting. XX.AG start.linear logical flag. If TRUE, the model starts with the linear fit. XX.AG fit.object This the object returned by `bruto()'; if supplied, the same model is fit to the presumeably new y. XX.RT a multiresponse additive model fit object of class `bruto' is returned. The model is fit by adaptive backfitting using smoothing splines. If there are `np' columns in `y', then `np' additive models are fit, but the same amount of smoothing (df) is used for the jth term of each. The procedure chooses between `df = 0' (term omitted), `df = 1' (term linear) or `df > 0' (term fitted by smoothing spline. The model selection is based on an approximation to the GCV criterion, which is used at each step of the backfitting procedure. Once the selection process stops, the model is backfit using the chosen amount of smoothing. XX.PP XX A bruto object has the following components of interest: XX.RC lambda a vector of chosen smoothing parameters, one for each column of x XX.RC df the df chosen for each column of x XX.RC type a factor with levels `excluded', `linear' or `smooth', indicating the status of each column of x. XX.RC gcv.select XX.RC gcv.backfit XX.RC df.select The sequence of gcv values and df selected during the execution of the function. XX.RC nit The number of iterations used XX.RC fitted.values a matrix of fitted values XX.RC residuals a matrix of residuals XX.RC call the call that produced this object XX.SH REFERENCES Trevor Hastie and Rob Tibshirani, XX.ul Generalized Additive Models, Chapman and Hall, 1990 (page 262). Trevor Hastie, Rob Tibshirani and Andreas Buja "Flexible Discriminant Analysis by Optimal Scoring" AT\&T Bell Laboratories Technical Memorandum, February 1993. XX.SA `predict.bruto' XX.EX bruto(x,y) XX.WR XX EEEEE X $shar_touch -am 0410235395 '.Data/.Help/bruto' && X chmod 0644 '.Data/.Help/bruto' || X $echo 'restore of' '.Data/.Help/bruto' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/bruto:' 'MD5 check failed' c3d8690a562f1e1dfac6755a6f5a128f .Data/.Help/bruto XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/bruto'`" X test 2695 -eq "$shar_count" || X $echo '.Data/.Help/bruto:' 'original size' '2695,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/confusion ============== if test -f '.Data/.Help/confusion' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/confusion' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/confusion' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/confusion' && XX.BG XX.FN confusion XX.TL compute the confusion matrix between two factors, or for an fda or mda object XX.CS confusion(predict, true) XX.AG predict the predicted factor, or an fda or mda model XX.AG true the true factor, or else (in the case where predict is a model object) a data frame (list) containing the test data. XX.RT essential table(predict, true), but with some useful attribute(s) XX.SA `fda', `predict.fda' XX.EX > confusion(predict(irisfit,iris$x), iris$g) XX Setosa Versicolor Virginica XX Setosa 50 0 0 Versicolor 0 48 1 XX Virginica 0 2 49 attr(, "error"): [1] 0.02 XX.WR EEEEE X $shar_touch -am 0410235395 '.Data/.Help/confusion' && X chmod 0644 '.Data/.Help/confusion' || X $echo 'restore of' '.Data/.Help/confusion' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/confusion:' 'MD5 check failed' afe4abac5dce7a262d03c90da4ed3676 .Data/.Help/confusion XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/confusion'`" X test 643 -eq "$shar_count" || X $echo '.Data/.Help/confusion:' 'original size' '643,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/fda ============== if test -f '.Data/.Help/fda' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/fda' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/fda' '(binary)' X sed 's/^X//' << 'EEEEE' | uudecode && begin 600 .Data/.Help/fda M+D)'"BY&3B!F9&$*+E1,"F9D83H@9FQE>&EB;&4@9&ES8W)I;6EN86YT(&%N M86QY7YL;V"DK>B<@971C("AT>7!E"F`_ M9F]R;75L82<@9F]R(&UO'!L86EN960@8GD@96%C:"!D:6UE M;G-I;VX@*')E;&%T:79E('1O('1H92!T;W1A;"!E>'!L86EN960N*2`*+E)# M('9A;'5E7-I"!W:&5N(&-L87-S:69Y:6YG('1H M92!T"<@86YD(&!Y)PIW:&5R M92!B;W1H(&-A;B!B92!M871R:6-E2!C86X@=&%K92!A9&1I=&EO;F%L"F%R9W5M96YT7)E9R2!(87-T:64L"E1I8G-H:7)A;FD@86YD($)U:F$L M(#$Y.30L($I!4T$L(#$R-34M,3(W,"X*(E!E;F%L:7IE9"!$:7-C&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/fda:' 'MD5 check failed' a279315a875681dbac72f4f6e9f552ce .Data/.Help/fda XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/fda'`" X test 4562 -eq "$shar_count" || X $echo '.Data/.Help/fda:' 'original size' '4562,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/gen.ridge ============== if test -f '.Data/.Help/gen.ridge' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/gen.ridge' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/gen.ridge' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/gen.ridge' && XX.BG XX.FN gen.ridge XX.TL Perform a penalized regression, as used by pda() XX.CS gen.ridge(x, y, weights, lambda=1, omega, df, ...) XX.AG x XX.AG y XX.AG weights the x and y matrix and possibly a weight vector XX.AG lambda the shrinkage penalty coefficient XX.AG omega a penalty object; omega is the eigendecomposition of the penalty matrix, and need not have full rank. By default, standard ridge is used. XX.AG df an alternative way to prescribe lambda, using the notion of equivalent degrees of freedom XX.RT A generalized ridge regression, where the coefficients are penalized according to omega. See the function definition for further details. No functions are provided for producing one dimensional penalty objects (omega). laplacian() creates a two dimensional penalty object, suitable for (small) images. XX.SA `laplacian' XX.WR EEEEE X $shar_touch -am 0410235395 '.Data/.Help/gen.ridge' && X chmod 0644 '.Data/.Help/gen.ridge' || X $echo 'restore of' '.Data/.Help/gen.ridge' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/gen.ridge:' 'MD5 check failed' 2e9c6325160f037a7725c575d86e5dc9 .Data/.Help/gen.ridge XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/gen.ridge'`" X test 814 -eq "$shar_count" || X $echo '.Data/.Help/gen.ridge:' 'original size' '814,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/junk ============== if test -f '.Data/.Help/junk' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/junk' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/junk' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/junk' && XX.BG XX.FN junk XX.TL XX~function to do ??? XX.DN XX~brief description XX.CS junk(formula=formula(data), data=sys.parent(), weights, theta, dimension=J - 1, eps=.Machine$single.eps, method=polyreg, keep.fitted=(n * dimension < 1000), ...) XX.RA XX.OA XX~move the above line to just above the first optional argument XX.AG formula XX~Describe formula here XX.AG data XX~Describe data here XX.AG weights XX~Describe weights here XX.AG theta XX~Describe theta here XX.AG dimension XX~Describe dimension here XX.AG eps XX~Describe eps here XX.AG method XX~Describe method here XX.AG keep.fitted XX~Describe keep.fitted here XX.AG ... XX~Describe ... here XX.RT XX~Describe the value returned XX.SE XX~describe any side effects if they exist XX.DT XX~explain details here. XX.SH REFERENCES XX~put references here, make other sections like NOTE and WARNING with .SH XX.SA XX~put functions to SEE ALSO here XX.EX # The function is currently defined as function(formula = formula(data), data = sys.parent(), XX weights, theta, dimension = J - 1, eps = XX .Machine$single.eps, method = polyreg, XX keep.fitted = (n * dimension < 1000), ...) { # # Flexible Discriminant Analysis # Function for fitting models described in # Hastie, Tibshirani and Buja, 1994, JASA # "Flexible Discriminant Analysis by Optimal Scoring" # Hastie, Buja and Tibshirani, 1995, Annals of Statistics # "Penalized Discriminant Analysis" # Modified 2/15/95 by T Hastie # XX this.call <- match.call() XX # This extracts the x and g from the formula or data frame # # -------< not for human consumption >-------- XX m <- match.call(expand = F) XX m[[1]] <- as.name("model.frame") XX m <- m[match(names(m), c("", "formula", "data", XX "weights"), 0)] XX m <- eval(m, sys.parent()) XX Terms <- attr(m, "terms") XX g <- model.extract(m, response) XX attr(Terms, "intercept") <- 0 XX x <- model.matrix(Terms, m) XX dd <- dim(x) XX n <- dd[1] XX m <- dd[2] # XX weights <- model.extract(m, weights) XX if(!length(weights)) XX weights <- rep(1, n) XX else if(any(weights < 0)) XX stop("negative weights not allowed") XX if(length(g) != n) XX stop("g should have length nrow(x)") XX fg <- factor(g) # # if some levels are missing, this gets rid of them # XX prior <- table(fg) XX prior <- prior/sum(prior) XX cnames <- levels(fg) XX g <- as.numeric(fg) XX J <- length(cnames) # # # construct indicator matrix for response # # Initialization uses orthogonal contrasts for theta # Even if contrasts are supplied, they need to be orthogonalized XX iswt <- F XX if(missing(weights)) XX dp <- table(g)/n XX else { XX weights <- (n * weights)/sum(weights) XX dp <- tapply(weights, g, sum)/n XX iswt <- T XX } XX if(missing(theta)) XX theta <- contr.helmert(J) XX theta <- contr.fda(dp, theta) XX Theta <- theta[g, , drop = F] # #Theta is now an n x K matrix of responses for the (nonparametric) # regression, normalized wrt the data (i.e. theta normalized wrt dp) # where K is min(J-1, ncol of starting theta) # XX fit <- method(x, Theta, weights, ...) XX if(iswt) XX Theta <- Theta * weights XX ssm <- t(Theta) %*% fitted(fit)/n # # This n could be (n-1) to make it "unbiassed" # # now we need the svd of ssm (unweighted) # XX ed <- svd(ssm, nu = 0) XX thetan <- ed$v XX lambda <- ed$d XX # Note: the discriminant eigenvalues are a transformation # of the optimal scaling values XX lambda[lambda > 1 - eps] <- 1 - eps # # If lambda is one we get errors, so we illiminate this possibility XX discr.eigen <- lambda/(1 - lambda) XX pe <- (100 * cumsum(discr.eigen))/sum( XX discr.eigen) XX dimension <- min(dimension, sum(lambda > eps)) XX if(dimension == 0) { XX warning("degenerate problem; no discrimination" XX ) XX return(structure(list(dimension = 0, XX fit = fit, call = this.call), XX class = "fda")) XX } XX thetan <- thetan[, seq(dimension), drop = F] # XX pe <- pe[seq(dimension)] XX #Now produce projected centroids # XX alpha <- sqrt(lambda[seq(dimension)]) XX sqima <- sqrt(1 - lambda[seq(dimension)]) # # package up results # XX vnames <- paste("v", seq(dimension), sep = "") XX means <- scale(theta %*% thetan, F, sqima/ XX alpha) XX dimnames(means) <- list(cnames, vnames) XX names(lambda) <- c(vnames, rep("", length( XX lambda) - dimension)) # XX names(pe) <- vnames XX obj <- structure(list(percent.explained = pe, XX values = lambda, means = means, XX theta.mod = thetan, dimension = XX dimension, prior = prior, fit = fit, XX call = this.call, terms = Terms), XX class = "fda") XX obj$confusion <- confusion(predict(obj), fg) XX # get rid of the fitted values; these take up too much space XX if(!keep.fitted) obj$fit$fitted.values <- NULL XX # XX obj } XX.KW ~keyword XX.WR EEEEE X $shar_touch -am 0229123996 '.Data/.Help/junk' && X chmod 0644 '.Data/.Help/junk' || X $echo 'restore of' '.Data/.Help/junk' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/junk:' 'MD5 check failed' b48cf3e9c4df192579147f3d6ac49cdf .Data/.Help/junk XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/junk'`" X test 4438 -eq "$shar_count" || X $echo '.Data/.Help/junk:' 'original size' '4438,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/junk~ ============== if test -f '.Data/.Help/junk~' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/junk~' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/junk~' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/junk~' && XX.BG XX.FN fda XX.TL XX~function to do ??? XX.DN XX~brief description XX.CS fda(formula=formula(data), data=sys.parent(), weights, theta, dimension=J - 1, eps=.Machine$single.eps, method=polyreg, keep.fitted=(n * dimension < 1000), ...) XX.RA XX.OA XX~move the above line to just above the first optional argument XX.AG formula XX~Describe formula here XX.AG data XX~Describe data here XX.AG weights XX~Describe weights here XX.AG theta XX~Describe theta here XX.AG dimension XX~Describe dimension here XX.AG eps XX~Describe eps here XX.AG method XX~Describe method here XX.AG keep.fitted XX~Describe keep.fitted here XX.AG ... XX~Describe ... here XX.RT XX~Describe the value returned XX.SE XX~describe any side effects if they exist XX.DT XX~explain details here. XX.SH REFERENCES XX~put references here, make other sections like NOTE and WARNING with .SH XX.SA XX~put functions to SEE ALSO here XX.EX # The function is currently defined as function(formula = formula(data), data = sys.parent(), XX weights, theta, dimension = J - 1, eps = XX .Machine$single.eps, method = polyreg, XX keep.fitted = (n * dimension < 1000), ...) { # # Flexible Discriminant Analysis # Function for fitting models described in # Hastie, Tibshirani and Buja, 1994, JASA # "Flexible Discriminant Analysis by Optimal Scoring" # Hastie, Buja and Tibshirani, 1995, Annals of Statistics # "Penalized Discriminant Analysis" # Modified 2/15/95 by T Hastie # XX this.call <- match.call() XX # This extracts the x and g from the formula or data frame # # -------< not for human consumption >-------- XX m <- match.call(expand = F) XX m[[1]] <- as.name("model.frame") XX m <- m[match(names(m), c("", "formula", "data", XX "weights"), 0)] XX m <- eval(m, sys.parent()) XX Terms <- attr(m, "terms") XX g <- model.extract(m, response) XX attr(Terms, "intercept") <- 0 XX x <- model.matrix(Terms, m) XX dd <- dim(x) XX n <- dd[1] XX m <- dd[2] # XX weights <- model.extract(m, weights) XX if(!length(weights)) XX weights <- rep(1, n) XX else if(any(weights < 0)) XX stop("negative weights not allowed") XX if(length(g) != n) XX stop("g should have length nrow(x)") XX fg <- factor(g) # # if some levels are missing, this gets rid of them # XX prior <- table(fg) XX prior <- prior/sum(prior) XX cnames <- levels(fg) XX g <- as.numeric(fg) XX J <- length(cnames) # # # construct indicator matrix for response # # Initialization uses orthogonal contrasts for theta # Even if contrasts are supplied, they need to be orthogonalized XX iswt <- F XX if(missing(weights)) XX dp <- table(g)/n XX else { XX weights <- (n * weights)/sum(weights) XX dp <- tapply(weights, g, sum)/n XX iswt <- T XX } XX if(missing(theta)) XX theta <- contr.helmert(J) XX theta <- contr.fda(dp, theta) XX Theta <- theta[g, , drop = F] # #Theta is now an n x K matrix of responses for the (nonparametric) # regression, normalized wrt the data (i.e. theta normalized wrt dp) # where K is min(J-1, ncol of starting theta) # XX fit <- method(x, Theta, weights, ...) XX if(iswt) XX Theta <- Theta * weights XX ssm <- t(Theta) %*% fitted(fit)/n # # This n could be (n-1) to make it "unbiassed" # # now we need the svd of ssm (unweighted) # XX ed <- svd(ssm, nu = 0) XX thetan <- ed$v XX lambda <- ed$d XX # Note: the discriminant eigenvalues are a transformation # of the optimal scaling values XX lambda[lambda > 1 - eps] <- 1 - eps # # If lambda is one we get errors, so we illiminate this possibility XX discr.eigen <- lambda/(1 - lambda) XX pe <- (100 * cumsum(discr.eigen))/sum( XX discr.eigen) XX dimension <- min(dimension, sum(lambda > eps)) XX if(dimension == 0) { XX warning("degenerate problem; no discrimination" XX ) XX return(structure(list(dimension = 0, XX fit = fit, call = this.call), XX class = "fda")) XX } XX thetan <- thetan[, seq(dimension), drop = F] # XX pe <- pe[seq(dimension)] XX #Now produce projected centroids # XX alpha <- sqrt(lambda[seq(dimension)]) XX sqima <- sqrt(1 - lambda[seq(dimension)]) # # package up results # XX vnames <- paste("v", seq(dimension), sep = "") XX means <- scale(theta %*% thetan, F, sqima/ XX alpha) XX dimnames(means) <- list(cnames, vnames) XX names(lambda) <- c(vnames, rep("", length( XX lambda) - dimension)) # XX names(pe) <- vnames XX obj <- structure(list(percent.explained = pe, XX values = lambda, means = means, XX theta.mod = thetan, dimension = XX dimension, prior = prior, fit = fit, XX call = this.call, terms = Terms), XX class = "fda") XX obj$confusion <- confusion(predict(obj), fg) XX # get rid of the fitted values; these take up too much space XX if(!keep.fitted) obj$fit$fitted.values <- NULL XX # XX obj } XX.KW ~keyword XX.WR EEEEE X $shar_touch -am 0229123896 '.Data/.Help/junk~' && X chmod 0644 '.Data/.Help/junk~' || X $echo 'restore of' '.Data/.Help/junk~' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/junk~:' 'MD5 check failed' 943c2a201c992b8942d3dc9aed660b69 .Data/.Help/junk~ XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/junk~'`" X test 4436 -eq "$shar_count" || X $echo '.Data/.Help/junk~:' 'original size' '4436,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/mars ============== if test -f '.Data/.Help/mars' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/mars' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/mars' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/mars' && XX.BG XX.FN mars XX.TL mars: Multivariate Additive Regression Splines XX.CS mars(x, y, w, wp, degree, nk, penalty, thresh, prune, trace, forward.step, prevfit) XX.PP XX.AG x a matrix containing the independent variables. These should not include the constant term( a column of all ones). XX.AG y a vector containing the response variable, or in the case of multiple responses, a matrix whose columns are the response values for each variable XX.AG w an optional vector of observation weights. XX.AG wp an optional vector of response weights. XX.AG degree an optional integer specifying maximum interaction degree (default is 1) XX.AG nk an optional integer specifying the maximum number of model terms XX.AG penalty an optional value specifying the cost per degree of freedom charge (default is 2) XX.AG thresh an optional value specifying forward stepwise stopping threshold (default is 0.001---see dmarss.r file for details). XX.AG prune an optional logical value specifying whether the model should be pruned in a backward stepwise fashion (default is TRUE). XX.AG trace an optional logical value specifying whether info should be printed along the way (default is FALSE). XX.AG forward.step an optional logical value specifying whether forward stepwise process should be carried out (default is TRUE). XX.AG prevfit optional data structure from previous fit. To see the effect of changing the penalty paramater, one can use prevfit with `forward.step=F' XX XX.RT structure with the following components: XX.RC call call used to `mars()' XX.RC all.terms term numbers in full model. 1=constant term. remaining terms are in pairs- 2 3, 4 5 etc. all.terms indicates nonsingular set of terms XX.RC selected.terms term numbers in selected model. XX.RC penalty the input penalty value XX.RC degree the input degree value XX.RC thresh the input threshold value XX.RC gcv gcv of chosen model XX.RC factor XX matrix with ijth element =1 if term i has a factor of the form xj > c, =-1 XX if term i has a factor of the form xj <= c, and 0 if xj is not in term i XX.RC cuts XX matrix with ijth element = cut point c for var j in term i XX.RC residuals XX residuals from fit XX.RC fitted values fitted values from fit XX.RC lenb length of full model XX.RC coefficients least squares coefficients for final model XX.RC x the input x matrix. XX.PP This function was coded from scratch, and did not use any of Friedman's mars code. It gives quite similar results to Friedman's program in our tests, but not exactly the same results. We have not implemented Friedman's anova decomposition nor are categorical predictors handled properly yet. Our version does handle multiple response variables, however. As it is not well-tested, we would like to hear of any bugs. XX.PP XX Trevor Hastie and Robert Tibshirani XX XX.SH REFERENCE J. Friedman "Multivariate Additive Regression Splines". Annals of Statistic 1991 XX.SA `predict.mars', `model.matrix.mars' XX.EX a <- mars(x,y) aa <- predict.mars(a,x) # get fitted values for predictor matrix x XX.KW non-parametric regression XX.WR EEEEE X $shar_touch -am 1127211697 '.Data/.Help/mars' && X chmod 0644 '.Data/.Help/mars' || X $echo 'restore of' '.Data/.Help/mars' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/mars:' 'MD5 check failed' d014abc4778662fe68428d62bfc03df6 .Data/.Help/mars XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/mars'`" X test 2983 -eq "$shar_count" || X $echo '.Data/.Help/mars:' 'original size' '2983,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/mars~ ============== if test -f '.Data/.Help/mars~' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/mars~' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/mars~' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/mars~' && XX.BG XX.FN mars XX.TL mars: Multivariate Additive Regression Splines XX.CS mars(x, y, w, wp, degree, nk, penalty, thresh, prune, trace, forward.step, prevfit) XX.PP XX.AG x a matrix containing the independent variables. XX.AG y a vector containing the response variable, or in the case of multiple responses, a matrix whose columns are the response values for each variable XX.AG w an optional vector of observation weights. XX.AG wp an optional vector of response weights. XX.AG degree an optional integer specifying maximum interaction degree (default is 1) XX.AG nk an optional integer specifying the maximum number of model terms XX.AG penalty an optional value specifying the cost per degree of freedom charge (default is 2) XX.AG thresh an optional value specifying forward stepwise stopping threshold (default is 0.001---see dmarss.r file for details). XX.AG prune an optional logical value specifying whether the model should be pruned in a backward stepwise fashion (default is TRUE). XX.AG trace an optional logical value specifying whether info should be printed along the way (default is FALSE). XX.AG forward.step an optional logical value specifying whether forward stepwise process should be carried out (default is TRUE). XX.AG prevfit optional data structure from previous fit. To see the effect of changing the penalty paramater, one can use prevfit with `forward.step=F' XX XX.RT structure with the following components: XX.RC call call used to `mars()' XX.RC all.terms term numbers in full model. 1=constant term. remaining terms are in pairs- 2 3, 4 5 etc. all.terms indicates nonsingular set of terms XX.RC selected.terms term numbers in selected model. XX.RC penalty the input penalty value XX.RC degree the input degree value XX.RC thresh the input threshold value XX.RC gcv gcv of chosen model XX.RC factor XX matrix with ijth element =1 if term i has a factor of the form xj > c, =-1 XX if term i has a factor of the form xj <= c, and 0 if xj is not in term i XX.RC cuts XX matrix with ijth element = cut point c for var j in term i XX.RC residuals XX residuals from fit XX.RC fitted values fitted values from fit XX.RC lenb length of full model XX.RC coefficients least squares coefficients for final model XX.RC x the input x matrix. XX.PP This function was coded from scratch, and did not use any of Friedman's mars code. It gives quite similar results to Friedman's program in our tests, but not exactly the same results. We have not implemented Friedman's anova decomposition nor are categorical predictors handled properly yet. Our version does handle multiple response variables, however. As it is not well-tested, we would like to hear of any bugs. XX.PP XX Trevor Hastie and Robert Tibshirani XX XX.SH REFERENCE J. Friedman "Multivariate Additive Regression Splines". Annals of Statistic 1991 XX.SA `predict.mars', `model.matrix.mars' XX.EX a <- mars(x,y) aa <- predict.mars(a,x) # get fitted values for predictor matrix x XX.KW non-parametric regression XX.WR EEEEE X $shar_touch -am 0410235395 '.Data/.Help/mars~' && X chmod 0644 '.Data/.Help/mars~' || X $echo 'restore of' '.Data/.Help/mars~' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/mars~:' 'MD5 check failed' 47ada327ae254c9e103d1c85612cef86 .Data/.Help/mars~ XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/mars~'`" X test 2916 -eq "$shar_count" || X $echo '.Data/.Help/mars~:' 'original size' '2916,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/mda ============== if test -f '.Data/.Help/mda' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/mda' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/mda' '(binary)' X sed 's/^X//' << 'EEEEE' | uudecode && begin 600 .Data/.Help/mda M+D)'"BY&3B!M9&$*+E1,"DUI>'1U7YX)R!I="!D97-C0IV96-T;W(@=&AA="!C86X@8F4@8V]E2!T2!A(&-A;&P@=&\@8&UD82YS=&%R="AX+"!G+"!S=6)C M;&%S2<@=V5I9VAT M2!F=6YC=&EO;G,@=7-E9G5L M(&9O2X*+E)#(&1I;65N6EN9R!T:&4@=')A:6YI;F<@9&%T80HN M4D,@=V5I9VAT2!O9B!T:&4@;65T:&]D*"D@=7-E9"X*5&AE(&1E=FEA;F-E(&-A;B!I;F-R M96%S92!W:71H('1H92!I=&5R871I;VYS+"!B=70@9V5N97)A;&QY(&1O97,@ M;F]T+@HN4%`*5&AE(&!M971H;V0G(&9U;F-T:6]N2<*=VAE22!A7)E9R@I)R!H87,@82!@9&5G6YO;6EA;"!R96=R97-S:6]N(&]F('1H92!R97%U:7)E9"!T;W1A M;"!D96=R964N(`I3964@=&AE(&1O8W5M96YT871I;VX@9F]R(&!P2!B=6=S+@HN4%`*(%1R M979O2!(87-T:64L"E1I8G-H:7)A M;FD@86YD($)U:F$L(#$Y.30L($I!4T$L(#$R-34M,3(W,"X*(E!E;F%L:7IE M9"!$:7-C"P@9&%T M82`](&=L87-S+G1R86EN*0H^(&=L87-S+FUD80I#86QL.@IM9&$H9F]R;75L M82`](&<@?B!X+"!D871A(#T@9VQA7!E/2)P;W-T(BD@(R!A8F)R979I871I;VYS(&%R M92!A;&QO=V5D"CX@8V]N9G5S:6]N*&=L87-S+FUD82QG;&%S&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/mda:' 'MD5 check failed' 99c64d175092dd946cecd12f662a6ae6 .Data/.Help/mda XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/mda'`" X test 6491 -eq "$shar_count" || X $echo '.Data/.Help/mda:' 'original size' '6491,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/mda.start ============== if test -f '.Data/.Help/mda.start' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/mda.start' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/mda.start' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/mda.start' && XX.BG XX.FN mda.start XX.TL Provide starting weights for the mda function XX.CS mda.start(x, g, subclasses=3, trace.mda.start=F, start.method=c("kmeans", "lvq"), tries=5, criterion=c("misclassification", "deviance"), ...) XX.AG x The x data, or an mda object XX.AG g The response vector g XX.AG subclasses number of subclasses per class, as in mda. XX.AG trace.mda.start Show results of each iteration XX.AG start.method Either "kmeans" or "lvq". The latter requires the statlib contribution "classif". XX.AG tries Number of random starts XX.AG criterion By default, classification errors on the training data. Posterior deviance is also an option. XX.RT A list of weight matrices, one for each class XX.WR EEEEE X $shar_touch -am 0410235395 '.Data/.Help/mda.start' && X chmod 0644 '.Data/.Help/mda.start' || X $echo 'restore of' '.Data/.Help/mda.start' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/mda.start:' 'MD5 check failed' e06388bf9b835feee488baeb1b1bf6f6 .Data/.Help/mda.start XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/mda.start'`" X test 681 -eq "$shar_count" || X $echo '.Data/.Help/mda.start:' 'original size' '681,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/model.matrix.mars ============== if test -f '.Data/.Help/model.matrix.mars' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/model.matrix.mars' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/model.matrix.mars' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/model.matrix.mars' && XX.BG XX.FN model.matrix.mars XX.TL produce a model matrix from a `mars' object XX.CS model.matrix.mars(object, x, which, full=F, ...) XX.AG object a `mars' object XX.AG x optional argument; if supplied, the mars basis functions are evaluated at these new observations. XX.AG which which columns should be used. The default is to use the columns described by the component `selected.terms' on `object' XX.AG full if `TRUE' the entire set of columns are selected, even redundant ones. This is used for updating a mars fit XX.RT a model matrix corresponding to the selected columns XX.SA `mars', `predict.mars' XX.WR EEEEE X $shar_touch -am 0410235395 '.Data/.Help/model.matrix.mars' && X chmod 0644 '.Data/.Help/model.matrix.mars' || X $echo 'restore of' '.Data/.Help/model.matrix.mars' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/model.matrix.mars:' 'MD5 check failed' b661cee6aca65240a8ac0f0dae3382ec .Data/.Help/model.matrix.mars XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/model.matrix.mars'`" X test 593 -eq "$shar_count" || X $echo '.Data/.Help/model.matrix.mars:' 'original size' '593,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/polyreg ============== if test -f '.Data/.Help/polyreg' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/polyreg' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/polyreg' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/polyreg' && XX.BG XX.FN polyreg XX.TL simple minded polynomial regression XX.CS polyreg(x, y, w, degree=1, monomial=F, ...) XX.PP XX.AG x predictor matrix XX.AG y response matrix XX.AG w optional (positive) weights XX.AG degree total degree of polynomial basis (default is 1) XX.AG monomial if `TRUE' a monomial basis is used (no cross terms). default is `FALSE' XX.RT a polynomial regression fit, containing the essential ingredients for its predict method. XX.SA `predict.polyreg' XX.WR EEEEE X $shar_touch -am 0410235395 '.Data/.Help/polyreg' && X chmod 0644 '.Data/.Help/polyreg' || X $echo 'restore of' '.Data/.Help/polyreg' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/polyreg:' 'MD5 check failed' 6d884ec6cb0fb09d1fa36f3848f7dc4f .Data/.Help/polyreg XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/polyreg'`" X test 451 -eq "$shar_count" || X $echo '.Data/.Help/polyreg:' 'original size' '451,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/predict.bruto ============== if test -f '.Data/.Help/predict.bruto' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/predict.bruto' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/predict.bruto' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/predict.bruto' && XX.BG XX.FN predict.bruto XX.TL a method for making predictions from a `bruto' object. XX.CS predict.bruto(object, x, type=c("fitted", "terms")) XX.AG object a fitted `bruto' object XX.AG x values at which predictions are to be made. XX.AG type if type is `fitted', the fitted values are returned. If type is `terms', a list of fitted terms is returned, each with an `x' and `y' component. These can be used to show the fitted functions. XX.RT either a fit matrix or a list of fitted terms. XX.SA `bruto', `predict' XX.EX fit <- bruto(x, y) fitted.terms <- predict(fit, newx, type = "terms") for( tt in fitted.terms) plot(tt) XX.WR EEEEE X $shar_touch -am 0410235395 '.Data/.Help/predict.bruto' && X chmod 0644 '.Data/.Help/predict.bruto' || X $echo 'restore of' '.Data/.Help/predict.bruto' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/predict.bruto:' 'MD5 check failed' f3b8933148de26452c71b6ad83b61167 .Data/.Help/predict.bruto XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/predict.bruto'`" X test 610 -eq "$shar_count" || X $echo '.Data/.Help/predict.bruto:' 'original size' '610,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/predict.fda ============== if test -f '.Data/.Help/predict.fda' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/predict.fda' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/predict.fda' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/predict.fda' && XX.BG XX.FN predict.fda XX.TL make predictions from an fda object XX.CS predict.fda(object, x, type, prior, dimension) XX.AG object An object of class `fda'. XX.AG x new data at which to make predictions. If missing, the training data is used. XX.AG type kind of predictions: `type = class' (default) produces a fitted factor, `type = variates' produces a matrix of discriminant variables, `type = posterior' produces a matrix of posterior probabilities (based on a gaussian assumption), and `type=hierarchical' produces the predicted class in sequence for models of all dimensions. XX.AG prior the prior probabability vector for each class; the default is the training sample proportions. XX.AG dimension the dimension of the space to be used, no larger than the dimension component XX of `object'. XX.RT An appropriate object depending on `type'. `object' has a component `"fit'" which is regression fit produced by the `method' argument to `fda'. There should be a `predict' method for this object which is invoked. This method should itself take as input `object' and optionally `x'. XX.SA `fda', `mars', `bruto', `polyreg', `softmax', `confusion' XX.EX > irisfit <- fda(iris$x, iris$g) > irisfit Call: fda(x = iris$x, g = iris$g) XX Dimension: 2 XX Percent Between-Group Variance Explained: XX v1 v2 XX 99.12 100 > confusion(predict(irisfit,iris$x), iris$g) XX Setosa Versicolor Virginica XX Setosa 50 0 0 Versicolor 0 48 1 XX Virginica 0 2 49 attr(, "error"): [1] 0.02 XX.WR EEEEE X $shar_touch -am 0410235395 '.Data/.Help/predict.fda' && X chmod 0644 '.Data/.Help/predict.fda' || X $echo 'restore of' '.Data/.Help/predict.fda' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/predict.fda:' 'MD5 check failed' 7f9ff1bac3d2423df214118b139ac1f0 .Data/.Help/predict.fda XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/predict.fda'`" X test 1525 -eq "$shar_count" || X $echo '.Data/.Help/predict.fda:' 'original size' '1525,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/predict.mars ============== if test -f '.Data/.Help/predict.mars' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/predict.mars' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/predict.mars' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/predict.mars' && XX.BG XX.FN predict.mars XX.TL make predictions from a `mars' object XX.CS predict.mars(object, x) XX.AG object a `mars' object XX.AG x values at which predictions are to be made. XX.RT the fitted values XX.SA `mars', `predict', `model.matrix.mars' XX.WR EEEEE X $shar_touch -am 0410235395 '.Data/.Help/predict.mars' && X chmod 0644 '.Data/.Help/predict.mars' || X $echo 'restore of' '.Data/.Help/predict.mars' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/predict.mars:' 'MD5 check failed' d946255f7104f8dd805b306f17842db1 .Data/.Help/predict.mars XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/predict.mars'`" X test 237 -eq "$shar_count" || X $echo '.Data/.Help/predict.mars:' 'original size' '237,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/predict.mda ============== if test -f '.Data/.Help/predict.mda' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/predict.mda' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/predict.mda' '(text)' X sed 's/^X//' << 'EEEEE' > '.Data/.Help/predict.mda' && XX.BG XX.FN predict.mda XX.TL make predictions from an mda object XX.CS predict.mda(object, x, type, prior, dimension, ...) XX.AG object a fitted mda object XX.AG x new data at which to make predictions. If missing, the training data is used. XX.AG type kind of predictions: `type = class' (default) produces a fitted factor, `type = variates' produces a matrix of discriminant variables (note that the maximal dimension is determined by the number of subclasses), `type = posterior' produces a matrix of posterior probabilities (based on a gaussian assumption), `type=hierarchical' produces the predicted class in sequence for models of dimensions specified by `dimension' argument. XX.AG prior the prior probabability vector for each class; the default is the training sample proportions. XX.AG dimension the dimension of the space to be used, no larger than the dimension component XX of `object', and in general < R, the number of subclasses. Dimension XX can be a vector for use with `type="hierarchical"' XX.RT An appropriate object depending on `type'. `object' has a component `"fit'" which is regression fit produced by the `method' argument to `mda'. There should be a `predict' method for this object which is invoked. This method should itself take as input `object' and optionally `x'. XX.SA `mda',`fda', `mars', `bruto', `polyreg', `softmax', `confusion' XX.EX > glass.mda <- mda(g ~ x, data = glass.train) > glass.mda Call: mda(formula = g ~ x, data = glass.train) XX Dimension: 9 XX Percent Between-Group Variance Explained: XX v1 v2 v3 v4 v5 v6 v7 v8 v9 XX 68.13 86.45 91.5 95.85 97.92 99.13 99.58 99.96 100 XX Degrees of Freedom (per dimension): 10 XX Training Misclassification Error: 0.16854 ( N = 89 ) XX Deviance: 75.363 > predict(glass.mda, type="post") # abbreviations are allowed > confusion(glass.mda,glass.test) XX O VW WF WNF XX O 9 0 0 4 XX VW 0 3 4 3 XX WF 0 4 27 5 WNF 2 1 9 25 attr(, "error"): [1] 0.3333333 XX XX.KW classification XX.WR XX XX XX XX XX XX XX XX XX XX EEEEE X $shar_touch -am 0410235395 '.Data/.Help/predict.mda' && X chmod 0644 '.Data/.Help/predict.mda' || X $echo 'restore of' '.Data/.Help/predict.mda' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/predict.mda:' 'MD5 check failed' de6c336edbabd847a077c2c0b9224099 .Data/.Help/predict.mda XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/predict.mda'`" X test 1975 -eq "$shar_count" || X $echo '.Data/.Help/predict.mda:' 'original size' '1975,' 'current size' "$shar_count!" X fi fi # ============= .Data/.Help/softmax ============== if test -f '.Data/.Help/softmax' && test "$first_param" != -c; then X $echo 'x -' SKIPPING '.Data/.Help/softmax' '(file already exists)' else X $echo 'x -' extracting '.Data/.Help/softmax' '(binary)' X sed 's/^X//' << 'EEEEE' | uudecode && begin 600 .Data/.Help/softmax M+D)'"BY&3B!S;V9T;6%X"BY43`II9&5N=&EF>2!T:&4@;6%X:6UU;2!I;B!E M86-H(')O=R!O9B!A(&UA=')I>`HN0U,*"AX+"!G87`]1BD*+D%' M('@*82!M871R:7@@;V8@;6]D92!N=6UE7!E(#T@(G!O M&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo '.Data/.Help/softmax:' 'MD5 check failed' b9433aad8a02c407d50d5a63d3445934 .Data/.Help/softmax XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < '.Data/.Help/softmax'`" X test 590 -eq "$shar_count" || X $echo '.Data/.Help/softmax:' 'original size' '590,' 'current size' "$shar_count!" X fi fi rm -fr _sh20746 Xexit 0 SHAR_EOF $shar_touch -am 0425070698 'helpfile.shar' && chmod 0644 'helpfile.shar' || $echo 'restore of' 'helpfile.shar' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'helpfile.shar:' 'MD5 check failed' f233f91970cca8aab116aab590d844de helpfile.shar SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'helpfile.shar'`" test 60886 -eq "$shar_count" || $echo 'helpfile.shar:' 'original size' '60886,' 'current size' "$shar_count!" fi fi # ============= ratfor.shar ============== if test -f 'ratfor.shar' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'ratfor.shar' '(file already exists)' else $echo 'x -' extracting 'ratfor.shar' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'ratfor.shar' && #!/bin/sh # This is a shell archive (produced by GNU sharutils 4.2). # To extract the files from this archive, save it to some FILE, remove # everything before the `!/bin/sh' line above, then type `sh FILE'. # # Made on 1998-04-25 07:06 PDT by . # Source directory was `/f7/trevor/MDA'. # # Existing files will *not* be overwritten unless `-c' is specified. # # This shar contains: # length mode name # ------ ---------- ------------------------------------------ # 3441 -rw-r--r-- bruto.r # 778 -rw-r--r-- dcalcvar.r # 20070 -rw-r--r-- dmarss.r # 527 -rw-r--r-- dorthreg.r # 807 -rw-r--r-- dqrreg.r # 22241 -rw-r--r-- mspline.r # 24539 -rw-r--r-- sortdi.f # save_IFS="${IFS}" IFS="${IFS}:" gettext_dir=FAILED locale_dir=FAILED first_param="$1" for dir in $PATH do X if test "$gettext_dir" = FAILED && test -f $dir/gettext \ X && ($dir/gettext --version >/dev/null 2>&1) X then X set `$dir/gettext --version 2>&1` X if test "$3" = GNU X then X gettext_dir=$dir X fi X fi X if test "$locale_dir" = FAILED && test -f $dir/shar \ X && ($dir/shar --print-text-domain-dir >/dev/null 2>&1) X then X locale_dir=`$dir/shar --print-text-domain-dir` X fi done IFS="$save_IFS" if test "$locale_dir" = FAILED || test "$gettext_dir" = FAILED then X echo=echo else X TEXTDOMAINDIR=$locale_dir X export TEXTDOMAINDIR X TEXTDOMAIN=sharutils X export TEXTDOMAIN X echo="$gettext_dir/gettext -s" fi touch -am 1231235999 $$.touch >/dev/null 2>&1 if test ! -f 1231235999 && test -f $$.touch; then X shar_touch=touch else X shar_touch=: X echo X $echo 'WARNING: not restoring timestamps. Consider getting and' X $echo "installing GNU \`touch', distributed in GNU File Utilities..." X echo fi rm -f 1231235999 $$.touch # if mkdir _sh20785; then X $echo 'x -' 'creating lock directory' else X $echo 'failed to create lock directory' X exit 1 fi # ============= bruto.r ============== if test -f 'bruto.r' && test "$first_param" != -c; then X $echo 'x -' SKIPPING 'bruto.r' '(file already exists)' else X $echo 'x -' extracting 'bruto.r' '(text)' X sed 's/^X//' << 'EEEEE' > 'bruto.r' && subroutine bruto(x,n,q,y,p,w,knot,nkmax,nk,wp,match,nef, XX dfmax,cost,lambda,df,coef,type,xrange, XX gcvsel,gcvbak,dfit,maxit,nit, XX eta,resid,thresh,work,trace) implicit double precision(a-h,o-z) double precision x(n,q),y(n,p),w(n),knot(nkmax+4,q),wp(p), XX dfmax(q),cost,lambda(q),df(q),coef(nkmax*p,q),xrange(2,q), XX gcvsel(q,maxit(1)),gcvbak(q,maxit(2)),dfit(q,maxit(1)), XX eta(n,p),resid(n,p),thresh,work(1) integer n,q,p,nkmax,nk(q),match(n,q),nef(q),type(q),maxit(2), XX nit(2),ier,ntype logical trace, select XX #Compute residuals do j=1,p{ XX do i=1,n{resid(i,j)=y(i,j)-eta(i,j)} } XX #bruto backfitting with selection select=.true. call bakssp(select,x,n,q,y,p,w,knot,nkmax,nk,wp,match,nef, XX dfmax,cost,lambda,df,coef,type,xrange, XX gcvsel,dfit,maxit(1),nit(1),eta,resid,thresh*10d0,work,trace) #regular backfitting select=.false. call bakssp(select,x,n,q,y,p,w,knot,nkmax,nk,wp,match,nef, XX dfmax,cost,lambda,df,coef,type,xrange, XX gcvbak,dfit,maxit(2),nit(2),eta,resid,thresh,work,trace) do j=1,p{ XX do i=1,n{eta(i,j)=y(i,j)-resid(i,j)} } return end subroutine bakssp(select,x,n,q,y,p,w,knot,nkmax,nk,wp,match,nef, XX dfmax,cost,lambda,df,coef,type,xrange, XX gcv,dfit,maxit,nit,s,resid,thresh,work,trace) implicit double precision(a-h,o-z) double precision x(n,q),y(n,p),w(n),knot(nkmax+4,q),wp(p), XX dfmax(q),cost,lambda(q),df(q),coef(nkmax*p,q), XX xrange(2,q),gcv(q,maxit),dfit(q,maxit), XX s(n,p),resid(n,p),thresh,work(1) integer n,q,p,nkmax,nk(q),match(n,q),nef(q),type(q),maxit,nit,ier,ntype double precision dfoff, gcv0, ndf,gcv1,gcvrat,ndfoff,wmean,sbar,rss,tol logical center, select, trace center=.true. tol=1d-3 # this is the convergence criterion for gcv #remove the mean from the residuals, and compute rss rss=0d0 do j=1,p{ XX sbar=wmean(n,resid(1,j),w) XX do i=1,n{ XX resid(i,j)=resid(i,j) - sbar XX work(i)=resid(i,j)**2 XX } rss=rss+wp(j)*wmean(n,work,w) } dfoff=0 do k=1,q {dfoff=dfoff+df(k)} gcv1=rss/((1-(1+dfoff*cost)/n)**2) gcvrat=1d0 nit=0 while(nitthresh ){ XX gcv0=gcv1 XX nit=nit+1 XX do k=1,q{ XX gcv(k,nit)=gcv1 XX if(!select&type(k)==1)next # form partial residuals if necessary XX if(type(k)>1){ #get the fitted values XX call psspl2(x(1,k),n,p,knot(1,k),nk(k),xrange(1,k), XX coef(1,k),coef(1,k),s,0,type(k)) XX do j=1,p{ XX sbar=wmean(n,s(1,j),w) XX do i=1,n{resid(i,j)=resid(i,j) + s(i,j)-sbar} XX } XX } XX XX ndfoff=dfoff-df(k) XX if(select) {ntype=0} XX else {ntype=type(k)} XX call sspl2(x(1,k),resid,w,n,p,knot(1,k),nk(k),wp,match(1,k), XX nef(k),ndfoff,dfmax(k),cost,lambda(k),ndf, XX gcv1,coef(1,k),s,ntype,center, XX xrange(1,k),work,tol,ier) XX gcv(k,nit)=gcv1 XX if(select){ XX dfit(k,nit)=ndf XX df(k)=ndf XX dfoff=ndfoff+ndf XX type(k)=ntype XX } XX if(type(k)>1){ XX do j=1,p{ XX do i=1,n{ XX resid(i,j)=resid(i,j) - s(i,j) XX } XX } XX } XX } gcvrat=dabs(gcv1-gcv0)/gcv0 if(trace){ XX call intpr("outer iteration",15,nit,1) XX call dblepr("gcv ",5,gcv1,1) XX call dblepr("ratio",5,gcvrat,1) } } return end subroutine pbruto(x,n,q,ybar,p,knot,nkmax,nk,coef,type,xrange, XX eta,work) implicit double precision(a-h,o-z) double precision x(n,q),ybar(p),knot(nkmax+4,q),coef(nkmax*p,q), XX xrange(2,q),eta(n,p),work(n,p) integer n,q,p,nkmax,nk(q),type(q) #initialization do j=1,p{ XX do i =1,n {eta(i,j)=ybar(j)} } do k=1,q{ XX if(type(k)==1)next XX call psspl2(x(1,k),n,p,knot(1,k),nk(k),xrange(1,k), XX coef(1,k),coef(1,k),work,0,type(k)) XX do j=1,p{ XX do i=1,n{eta(i,j)=eta(i,j) + work(i,j)} XX } } return end EEEEE X $shar_touch -am 0405184995 'bruto.r' && X chmod 0644 'bruto.r' || X $echo 'restore of' 'bruto.r' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo 'bruto.r:' 'MD5 check failed' 66ee33d048c4fbb03e00fe45cd5b8654 bruto.r XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'bruto.r'`" X test 3441 -eq "$shar_count" || X $echo 'bruto.r:' 'original size' '3441,' 'current size' "$shar_count!" X fi fi # ============= dcalcvar.r ============== if test -f 'dcalcvar.r' && test "$first_param" != -c; then X $echo 'x -' SKIPPING 'dcalcvar.r' '(file already exists)' else X $echo 'x -' extracting 'dcalcvar.r' '(text)' X sed 's/^X//' << 'EEEEE' > 'dcalcvar.r' && subroutine calcvar(nx,n,px,qr,qrank,qpivot,cov,tmpcov,work) XX implicit double precision (a-h,o-z) integer n,px,qrank,qpivot(px) double precision qr(nx,px),cov(px,px), tmpcov(px,px),work(1) XX double precision dn, wmean, dsum,dwrss,xbar,xcov integer i,j,km # compute the unscaled coviance matrix for the linear coefficients do i=1,qrank{ XX do j=1,qrank{ XX tmpcov(i,j)=0d0 XX cov(i,j)=qr(i,j) XX } XX tmpcov(i,i)=1e0 XX } info=0 XX # # notice I have put px in as the row dim of cov # call dbksl(cov,px,qrank,tmpcov,px,info) XX do i=1,qrank{ XX do j=i,qrank{ XX dsum=0e0 XX km=max(i,j) XX for(k=km;k<=qrank;k=k+1){ XX dsum=dsum+tmpcov(i,k)*tmpcov(j,k) XX } XX tmpcov(i,j)=dsum XX tmpcov(j,i)=dsum XX } XX } # no do i=1,qrank{ XX do j=1,qrank{ XX cov(i,j)=tmpcov(i,j) XX } XX } return end XX EEEEE X $shar_touch -am 0405184995 'dcalcvar.r' && X chmod 0644 'dcalcvar.r' || X $echo 'restore of' 'dcalcvar.r' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo 'dcalcvar.r:' 'MD5 check failed' 3a6afb85142745ff3b56524089efb49f dcalcvar.r XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'dcalcvar.r'`" X test 778 -eq "$shar_count" || X $echo 'dcalcvar.r:' 'original size' '778,' 'current size' "$shar_count!" X fi fi # ============= dmarss.r ============== if test -f 'dmarss.r' && test "$first_param" != -c; then X $echo 'x -' SKIPPING 'dmarss.r' '(file already exists)' else X $echo 'x -' extracting 'dmarss.r' '(text)' X sed 's/^X//' << 'EEEEE' > 'dmarss.r' && subroutine marss(nx,n,p,nclass,y,x,w,tagx,maxorder,mmax,penalty,thresh,forwstep,interms,prune,bx,fullin,lenb, bestgcv, bestin, flag, cut,dir, res,alpha,beta, scrat,iscrat,trace) # # multivariate mars # # input # # integer nx- number of rows of x, bx, and tagx # integer n- number of observations # integer p- number of variables # integer nclass - number of classes (resposne variables) # double precision y(n,nclass)- matrix of response values # double precision x(nx,p) -matrix of predictors # double precision w(n)- weights (ge 0) # integer tagx(nx,p) - tag array for sorting x # integer maxorder- maxumum order of interaction (mi in Friedman's notation) # # integer mmax- maximum number of basis functions (nk in Friedman's notation) # if it is desired to add say 3 terms to an existing model, one can set mmax=interms+3 # and set interms appropriately (see below) # # double precision penalty- cost per knot selection (df in Friedman's notation) # # double precision thresh- forward stepwise threshold - forw step stops if (rss/rssnull).le.thresh) # to prevent numerical problems with overfit solutions # # logical forwstep- if true, forward stepwise is carried out. if false, procedure # starts with basis matrix bx and carried out backward stepwise. This saves # a great deal of time in choosing the penalty parameter # # integer interms- number of terms in existing model. for a fresh start interms=1. # # logical prune- true if pruning (backward stepwise deletion) is desired # # # # # # double precision bx(nx,mmax)- basis matrix- used on input only if interms>1 or forwstep=false # # integer fullin(mmax)- basis column indicator. columns marked with a 1 indicate active # columns- used on input only if interms>1 or forwstep=false # # integer lenb- number of active columns of bx - # used on input only if interms>1 or forwstep=false # # double precision res(n,nclass) - residuals - used as input only if interms>1, in which case they # should be the residuals from the current fitted model # # double precision scrat(1)- scratch array of length at least # 1+n+ 2*n*mmax +4*mmax*mmax +3*mmax + 3*mmax*nclass +3*nclass # +28n+51 (for linpack) # # integer iscrat(1)- scratch array of length at least 4*mmax # #logical trace - should there be a printing of progress # # # output # # double precision bx(nx,mmax)- basis matrix. only first lenb columns are active. first column is 1 # integer fullin(mmax)- basis column indicator. columns marked with a 1 indicate # the linearly indep set # integer lenb- number of active columns of bx # double precision bestgcv- gcv for best model # integer bestin(mmax) basis column indicator for best model. 1 indicates term is in model # integer flag(mmax,p)- ijth element is 1 if predictor j is in model term i # (zero otherwise). note for a zero column of bx (fullin=0), the row of flag # will be all 1, for a technical reason, it should be ignored of course. # # double precision cut(mmax,p) - ijth element cutpoint for variable j in model term i # double precision dir(mmax,p) - ijth element is 1 if factor involving predictor j in term i is # of the form (x-t)+ ; it is -1 if it has the form (t-x)+ # # double precision res(n,nclass) - residuals from final fit # double precision alpha(nclass)- intercepts from final fit # double precision beta(mmax,nclass)- slopes from final fit XX XX XX implicit double precision (a-h,o-z) XX integer nx, n,p,nclass,tagx(nx,p),maxorder,mmax,bestin(mmax),flag(mmax,p), fullin(mmax) double precision y(n,nclass),x(nx,p),w(n),bx(nx,mmax),bestgcv,cut(mmax,p),dir(mmax,p),res(nx,nclass), alpha(nclass),beta(mmax,nclass) XX double precision scrat(1), tem integer iscrat(1) logical forwstep, prune, trace, tracec common tracec tracec=trace XX XX XX # iscrat must be of len at least 3*mmax # scrat must be of length at least # 1+n+ 2*n*mmax +4*mmax*mmax +3*mmax + 3*mmax*nclass +3*nclass XX #1. bxorth #2. bxorthm #3. cov #4. covsy #5. ybar #6. scr1 #7. scr5 #8. scr6 #9. temp XX #10. bxsc #11. r #12. betasc, #13 #14 #15 work XX len1=n*mmax len2=mmax len3=mmax*mmax len4=mmax*nclass len5=nclass len6=mmax len7=mmax len8=nclass len9=n XX XX len10=n*mmax len11=mmax*mmax len12=mmax*nclass len13=mmax*mmax len14=mmax*mmax XX n1=1 n2=n1+len1 n3=n2+len2 n4=n3+len3 n5=n4+len4 n6=n5+len5 n7=n6+len6 n8=n7+len7 n9=n8+len8 n10=n9+len9 n11=n10+len10 n12=n11+len11 n13=n12+len12 n14=n13+len13 n15=n14+len14 XX XX XX XX call marsnew1(nx,n,p,nclass,y,x,w,tagx,maxorder,mmax,bx,bestgcv, bestin,fullin,lenb,flag, cut,dir,res, alpha,beta,penalty,thresh,forwstep,interms,prune, scrat,scrat(n2),scrat(n3),scrat(n4),scrat(n5), scrat(n6),scrat(n7),scrat(n8),scrat(n9), scrat(n10),scrat(n11),scrat(n12),scrat(n13),scrat(n14),scrat(n15), iscrat,iscrat(1+mmax),iscrat(1+2*mmax),iscrat(1+3*mmax)) XX return end XX subroutine marsnew1(nx,n,p,nclass,y,x,w,tagx,maxorder,mmax,bx,bestgcv, bestin,fullin,lenb,flag, cut,dir,res,alpha,beta,penalty,thresh,forwstep,interms,prune, bxorth,bxorthm,cov,covsy,ybar, scr1,scr5,scr6, temp, XX bxsc, r, betasc,varsc,var,work, termlen,in, tempin, qpivot) # input n,p,nclass,x,tagx,y,w,mmax,maxorder XX # nclass is # of response variables XX # output- bestin, bestgcv, flag, cut,dir, bx,res XX implicit double precision (a-h,o-z) XX integer n,nterms2,p,m,mmax,flag(mmax,p),v,tagx(nx,p),termlen(mmax), nclass,fullin(mmax) XX double precision cov(mmax,mmax),covsy(mmax,nclass), critmax,x(nx,p),bx(nx,mmax),bxorth(n,mmax),bxorthm(mmax),y(n,nclass),ybar(nclass), scr1(mmax),scr5(mmax),scr6(nclass) XX double precision temp(n),w(n), cut(mmax,p),dir(mmax,p),alpha(nclass),beta(mmax,nclass), bxsc(n,mmax), r(mmax,mmax), dofit, XX res(nx,nclass),betasc(mmax,nclass), varsc(mmax,mmax), var(mmax,mmax), stopfac,work(1) XX integer tempin(mmax), bestin(mmax),qrank, qpivot(mmax) XX logical forwstep,go, prune, newform, cvar, trace common trace XX # # forw stepwise stops if gcv/gcvnull > stopfac, to avoid adding lots of # noisy terms # double precision rtemp(4) integer itemp(4) XX tolbx=.01 stopfac=10.0 #stopfac=10e9 prevcrit=10e9 #write(6,*) "in marsnew mmax=",mmax, " forwstep=",forwstep, "interms=",interms XX #call intpr("M1",2,n,1) #call dblepr("penalty",7,penalty,1) #call intpr("maxo",4,maxorder,1) XX if(interms.eq.1) {dofit=0} XX else{dofit=0;do j=2,lenb {dofit=dofit+fullin(j)} XX nterms=interms } XX XX if(forwstep){ fullin(1)=1 do i=2,mmax{ XX fullin(i)=0 } XX do i=1,n { XX w(i)=1 } XX do i=1, mmax{ XX termlen(i)=0 XX do j=1, p{ XX flag(i,j)=0 XX cut(i,j)=0 }} XX XX XX nterms=1; nterms2=2; do i=1,n { XX bx(i,1)=1 XX bxorth(i,1)=1.0/dsqrt(dfloat(n)) } XX bxorthm(1)=1/dsqrt(dfloat(n)) XX do i=1,n { XX do j=1, mmax{ XX bx(i,j)=0.0 }} XX do i=1,n { XX bx(i,1)=1 } XX do k=1, nclass{ XX ybar(k)=0.0 XX do i=1,n { XX ybar(k)=ybar(k)+y(i,k)/n }} if(interms.eq.1){ rssnull=0.0 XX do k=1, nclass{ XX do i=1,n { XX rssnull=rssnull+(y(i,k)-ybar(k))**2 XX }} } else{ XX rssnull=0.0 XX do k=1, nclass{ XX do i=1,n { XX rssnull=rssnull+res(i,k)**2 XX }} XX } XX XX XX rss=rssnull XX cmm= (1+dofit) + penalty*(.5*dofit) gcvnull=(rssnull/n)/(1.0-cmm/n)**2 XX #write(6,*) "initial rss=", rssnull, " initial gcv=", gcvnull if(trace)call dblepr("initial rss=",11,rssnull,1) if(trace)call dblepr("initial gcv=",11,gcvnull,1) XX lenb=1 ii=interms-1 go=.true. while( (ii.lt.(mmax-1)).and.((rss/rssnull).gt.thresh).and.go){ XX ii=ii+2 XX # write(6,*) "bef addtrm",ii, nterms XX XX do i1=1, nterms{ XX do i2=1, nterms{ XX cov(i1,i2)=0 XX }} XX XX do j=1, nterms{ XX cov(j,j)=0.0 XX do i=1,n { XX cov(j,j)=cov(j,j)+(bxorth(i,j)-bxorthm(j))*(bxorth(i,j)-bxorthm(j)) XX } XX } XX XX XX XX XX do k=1,nclass{ XX XX do j=1, nterms{ XX covsy(j,k)=0.0 XX do i=1,n { XX covsy(j,k)=covsy(j,k)+(y(i,k)-ybar(k))*bxorth(i,j) XX } XX } XX } XX XX XX XX XX do ik=1,mmax { tempin(ik)=fullin(ik)} XX XX call addtrm(nx,bx,tempin,bxorth,bxorthm,p,n,nclass,rss,prevcrit, XX cov,covsy,y,ybar,x,tagx,w,termlen,mmax,tolbx, XX nterms,flag,maxorder,scr1,scr5,scr6,imax,jmax,kmax,critmax, newform,bxsc, r, betasc, temp) XX # check whether to accept term doftemp=dofit XX doftemp=doftemp+1 if((imax>1).and.(newform)) {doftemp=doftemp+1} XX XX XX XX temprss=rss-critmax cmm= (1+doftemp) + penalty*(.5*doftemp) gcv=(temprss/n)/(1.0-cmm/n)**2 XX go=.false. if(((critmax/rss).gt.thresh).and.((gcv/gcvnull).lt.stopfac)){ XX go=.true. XX dofit=doftemp XX XX rss=rss-critmax XX kk=tagx(imax,jmax) #write(6,256) jmax, imax ,kmax, critmax,x(kk,jmax), rss,gcv, dofit 256 format(" ","adding term"," jmax=",i3, " imax=",i3 ," kmax=",i3, " critmax= ",f8.2, " cutp=", f9.5," rss=",f8.2, " gcv=",f8.2, " dofit=",f9.3) # write(6,*) itemp(1)=jmax;itemp(2)=imax;itemp(3)=kmax; rtemp(1)=critmax;rtemp(2)=x(kk,jmax);rtemp(3)=rss;rtemp(4)=gcv XX if(trace) call intpr("adding term ",12,ii,1) if(trace)call intpr("var, sp index, parent",21,itemp,3) if(trace)call dblepr("critmax cut rss gcv",19,rtemp,4) XX prevcrit=critmax XX do j=1,p { XX flag(ii,j)=flag(kmax,j) XX flag(ii+1,j)=flag(kmax,j) XX cut(ii,j)=cut(kmax,j) XX cut(ii+1,j)=cut(kmax,j) XX dir(ii,j)=dir(kmax,j) XX dir(ii+1,j)=dir(kmax,j) XX XX XX } XX termlen(ii)=termlen(kmax)+1 XX termlen(ii+1)=termlen(kmax)+1 XX XX XX XX XX XX do i=1,n { temp(i)=x(tagx(i,jmax),jmax)} XX temp1=temp(imax) XX XX fullin(ii)=1; XX XX if((imax.gt.1).and.(newform)) {fullin(ii+1)=1 } XX XX flag(ii,jmax)=1; XX flag(ii+1,jmax)=1; XX cut(ii,jmax)=temp1 XX cut(ii+1,jmax)=temp1 XX dir(ii,jmax)=1 XX dir(ii+1,jmax)=-1 # #this is to prevent trying to split a 0 column later XX XX if(fullin(ii+1).eq.0) XX { #do j=1,p {flag(ii+1,j)=1} # XX termlen(ii+1)=maxorder+1 XX } XX XX XX XX do i=1,n { if( (x(i,jmax)-temp1).gt.0) bx(i,ii)=bx(i,kmax)*(x(i,jmax)-temp1) XX if((temp1-x(i,jmax)).ge.0) bx(i,ii+1)=bx(i,kmax)*(temp1-x(i,jmax)) XX } XX XX XX XX XX if(nterms.eq.1){ XX temp1=0.0; do i=1,n { temp1=temp1+bx(i,2)/n; } XX do i=1,n {bxorth(i,2)=bx(i,2)-temp1;} XX } XX else{ ## write(6,*) "bef call to orthreg" XX XX XX call orthreg(n,n,nterms,bxorth,fullin, bx(1,ii),bxorth(1,nterms2)) XX } XX XX if(fullin(ii+1).eq.1){ XX XX XX call orthreg(n,n,nterms+1,bxorth,fullin, bx(1,ii+1),bxorth(1,nterms2+1)) # XX } XX else {do i=1,n {bxorth(i,nterms2+1)=0}} XX XX bxorthm(nterms2)=0.0 ; XX bxorthm(nterms2+1)=0.0 XX XX do i=1,n { XX bxorthm(nterms2)=bxorthm(nterms2)+bxorth(i,nterms2)/n XX bxorthm(nterms2+1)=bxorthm(nterms2+1)+bxorth(i,nterms2+1)/n XX } XX XX XX XX temp1=0.0;temp2=0.0;do i=1,n {temp1=temp1+bxorth(i,nterms2)**2 XX temp2=temp2+bxorth(i,nterms2+1)**2 } ##write(6,*) "norm of new columns",temp1,temp2 XX if (temp1.gt.0.0) {do i=1,n { bxorth(i,nterms2) =bxorth(i,nterms2)/dsqrt(temp1) }} XX if (temp2.gt.0.0) {do i=1,n { bxorth(i,nterms2+1)=bxorth(i,nterms2+1)/dsqrt(temp2) }} XX lenb=lenb+2 XX nterms=nterms+2; nterms2=nterms2+2; XX }} XX #write(6,*) "stopping forw stepwise" #if((rss/rssnull).le.thresh) write(6,*) "rss ratio=",rss/rssnull #if((critmax/rss).le.thresh) write(6,*) "crit ratio=",critmax/rss #if((gcv/gcvnull).gt.stopfac) write(6,*) "gcv ratio=",gcv/gcvnull XX rtemp(1)=rss/rssnull;rtemp(2)=critmax/rss;rtemp(3)=gcv/gcvnull if(trace)call dblepr("stopping forw step; rss crit and gcv ratios",43,rtemp,3) if(trace){ XX if((rss/rssnull).le.thresh) call dblepr("rss ratio=",10,rss/rssnull,1) XX if((critmax/rss).le.thresh) call dblepr ("crit ratio=",11,critmax/rss,1) call dblepr("critmax",7,critmax,1) call dblepr("rss",3,rss,1) XX if((gcv/gcvnull).gt.stopfac) call dblepr("gcv ratio=",10,gcv/gcvnull,1) } XX # write(6,*) "after forward step ","go=",go," rss=",rss, "nterms=", nterms XX XX } XX XX #call intpr("M2",2,fullin,10) #call intpr("nterms",6,nterms,1) dofit= -1 do i=1,nterms{ XX bestin(i)=fullin(i) XX dofit=dofit+fullin(i) } #call intpr("M21",3,n,1) XX XX if(trace)call intpr("aft forw step",13,nterms,1) XX XX call qrreg(nx,n,mmax,lenb,nclass,bx,bxsc,bestin,y,qpivot,qrank,beta,res,rss,cvar,var,varsc,scr1, work) XX XX #call dblepr("rss",3,rss,1) #call dblepr("beta",4,beta,4) #call intpr("M3",2,n,1) XX nt=dofit+1 if(qrank< nt) { #write(6,*) "singular matrix" XX do i=qrank+1,nt{ XX bestin(qpivot(i))=0 XX fullin(qpivot(i))=0 #write(6,*) "removing col ",qpivot(i) XX dofit=dofit-1; } } XX XX ##write(6,*) "bef mlin", nterms,dofit,(bestin(j),j=1,nterms),y(1,1),n,mmax,nclass,w(1) XX cvar=.true. #call intpr("M23",3,fullin,11) #call intpr("bestin",6,bestin,11) # #call intpr("lenb",4,lenb,1) #call intpr("mmax",4,mmax,1) #call dblepr("res",3,res,10) #call dblepr("bx",2,bx(1,1),5) #call dblepr("bx",2,bx(1,2),5) #call dblepr("bx",2,bx(n,11),1) #call dblepr("y",1,y,5) #call intpr("nx",2,nx,1) #call intpr("n",1,n,1) XX XX XX XX #write(6,*) "rank=",qrank, "rss=",rss XX XX XX XX XX XX XX XX rssfull=rss XX # dofit is the number of nonconstant, indep basis fns in the model XX XX XX XX cmm= (1+dofit) + penalty*(.5*dofit) bestgcv=(rss/n)/(1.0-cmm/n)**2 XX #write(6,*) "full model,gcv=",bestgcv, " rss=",rssfull, "dofit=",dofit rtemp(1)=bestgcv; rtemp(2)=rssfull;rtemp(3)=dofit if(trace)call dblepr("full model: gcv rss dofit",25,rtemp,3) if(trace)call intpr("terms",5,fullin,lenb) XX # after forward stepwise, nterms and lenb = # of columns of bx # dofit is number of nonzero columns XX XX XX if(prune){ XX do i=1,mmax{ tempin(i)=bestin(i)} #write(6,*) "initial tempin", (tempin(i),i=1,nterms) XX #call intpr("M5",2,n,1) while(dofit>0 ) { XX XX XX jo=1;rsstemp=10e99;minterm=0 XX do ii=2, lenb { XX XX if(tempin(ii).eq.1) { XX jo=jo+1 XX temp7=0.0 XX do kc=1,nclass{ XX temp7=temp7+beta(jo,kc)**2/var(jo,jo) XX } XX XX if(temp7 < rsstemp) {minterm=ii; rsstemp=temp7} XX XX XX XX }} XX rss=rss+rsstemp XX dofit=dofit-1; XX cmm= (1.0+dofit) + penalty*(.5*dofit) XX gcv=(rss/n)/(1.0-cmm/n)**2 XX XX tempin(minterm)=0 XX XX # write(6,100) minterm,gcv, rss,dofit,(tempin(ik),ik=1,lenb) 100 format(" ","pruning, minterm= ",i4, " gcv=",f9.3,2x, " rss=",f9.3,2x," dof=",f9.3, " model= ",60(i1,1x)) XX XX XX if(gcv< bestgcv) {bestgcv=gcv;do i=1,mmax {bestin(i)=tempin(i);}} XX XX XX if(dofit > 0) { cvar=.true. call qrreg(nx,n,mmax,lenb,nclass,bx,bxsc,tempin,y,qpivot,qrank,beta,res,rss,cvar,var,varsc,scr1,work) XX XX XX } } XX XX XX call qrreg(nx,n,mmax,lenb,nclass,bx,bxsc,bestin,y,qpivot,qrank,beta,res,rss,cvar,var,varsc, scr1, work) #write(6,101) bestgcv,rss,(bestin(ik),ik=1,lenb) 101 format(" ","best model gcv=",f9.3," rss=",f9.3,2x,"model= ",60(i1,1x)) if(trace)call intpr("best model",10,bestin,lenb) if(trace)call dblepr(" gcv=",4,bestgcv,1) XX XX } XX XX XX XX XX return end subroutine addtrm(nx,bx,tempin,bxorth,bxorthm,p,n,nclass, rss,prevcrit,cov,covsy,y,ybar,x,tagx,w,termlen,mmax,tolbx, nterms,flag,maxorder,scr1,scr5,scr6,imax,jmax,kmax,critmax, newform,bxsc,r,betasc, scrat) XX XX implicit double precision (a-h,o-z) XX integer n,nterms,nterms2,p,mmax,flag(mmax,p),v,tagx(nx,p),termlen(mmax), nclass, tempin(mmax), flagtot, minspan, iendspan XX double precision cov(mmax,mmax),covsy(mmax,nclass), critmax,x(nx,p),bx(nx,mmax),bxorth(n,mmax),bxorthm(mmax),y(n,nclass),ybar(nclass), scr1(mmax),scr5(mmax),scr6(nclass), bxsc(n,mmax), r(mmax,mmax), XX betasc(mmax,nclass), scrat(n),w(n) XX double precision temp1, temp2, scr2,sumb, sumbx, su, st, tem, temm logical newform, tnewform, trace common trace critmax=0;jmax=0;imax=0;kmax=0; XX XX # finds best term to add to current mars model XX XX XX XX XX XX XX ##write(6,*) "in addtrm n,p, nclass, nterms,mmax, y(1,1), x(1,1), tagx(1,1),flag(1,1),bx(1,1#),bxorth(1,1)", n,p, nclass, nterms,mmax, y(1,1), x(1,1), tagx(1,1),flag(1,1),bx(1,1),bxor#th(1,1) XX XX XX XX XX do m=1,nterms { XX XX nm=0 XX do jjj=1,n { XX if(bx(jjj,m).gt.0) {nm=nm+1} XX } XX XX tem=-(1d0/(n*nm))*dlog(1d0 - 5d-2) XX XX minspan= -1d0*(dlog(tem)/dlog(2d0))/2.5 XX XX tem=(5d-2)/n XX iendspan=3d0-dlog(tem)/dlog(2d0) XX # # write(6,*) "minspan,iendspan",minspan,iendspan XX if(termlen(m)< maxorder){ XX XX do v=1,p { XX XX if(flag(m,v).eq.0){ # check if model term type is already in model, to avoid a linear dependence XX XX XX tnewform=.true. XX mm=1 XX while((mm.le.nterms).and.tnewform) { XX mm=mm+1 XX if(tempin(mm).eq.1){ XX tnewform=.false. XX if(flag(mm,v).ne.1){tnewform=.true.;go to 9911} XX do j=1,p {if(j.ne.v) { if(flag(mm,j).ne.flag(m,j)) {tnewform=.true.;go to 9911 }}} XX } 9911 continue } XX # if new form of term, fit term bm*x # should it be bxorth below? XX if(tnewform) { XX nterms2=nterms+1 XX XX do i=1,n { XX scrat(i)=x(i,v)*bx(i,m) XX } # # XX if(nterms>1){ XX XX call orthreg(n,n,nterms,bxorth,tempin, scrat,bxorth(1,nterms2)) XX } else{ XX tem=0;do i=1,n {tem=tem+scrat(i)/n;} XX do i=1,n{ bxorth(i,2)=scrat(i)-tem} XX } XX XX XX bxorthm(nterms2)=0.0 ; XX do i=1,n { XX bxorthm(nterms2)=bxorthm(nterms2)+bxorth(i,nterms2)/n XX } XX temp1=0.0;do i=1,n {temp1=temp1+bxorth(i,nterms2)**2} XX if (temp1.gt.tolbx) {do i=1,n { bxorth(i,nterms2)=bxorth(i,nterms2)/dsqrt(temp1) }} XX else {do i=1,n {bxorth(i,nterms2)=0};tnewform=.false.;} XX XX do i1=1, nterms2{ XX cov(i1,nterms2)=0.0 XX cov(nterms2, i1)=0.0 XX } XX XX XX cov(nterms2,nterms2)=1 XX XX XX do kc=1,nclass{ XX covsy(nterms2,kc)=0.0 XX do i=1,n { XX covsy(nterms2,kc)=covsy(nterms2,kc)+(y(i,kc)-ybar(kc))*bxorth(i,nterms2) XX XX } XX } XX XX XX critnew=0.0 XX do kc=1,nclass { XX temp1=0 XX do i=1,n { temp1=temp1+y(i,kc)*bxorth(i,nterms2)} XX critnew=critnew+temp1**2 XX } ## write(6,*) "geoff",nterms,v,critnew XX if(critnew.gt.critmax) { XX jmax=v XX critmax=critnew XX imax=1 XX kmax=m XX } XX } XX XX if(tnewform) {nterms2=nterms+1; nterms21=nterms+2} else{nterms2=nterms; nterms21=nterms+1;critnew=0.0} XX #if((nterms.eq.7).and.(v.eq.2)) { # call dblepr("bxorth",6,bxorth(1,nterms2),10) # call dblepr("temp1",5,temp1,1) # call dblepr("critmax",7,critmax,1) XX #} XX XX XX XX XX # try other knot locations XX XX XX XX XX XX do kc=1, nclass{ XX covsy(nterms21,kc)=0 XX } XX do ii=1,nterms21 { XX cov(ii,nterms21)=0 XX cov(nterms21,ii)=0 XX } XX do kc=1,nclass { XX scr6(kc)=0 XX } XX do ii=1,nterms21{ XX scr1(ii)=0; XX } scr2=0;su=0;st=0;sumbx2=0;sumb=0.0;sumbx=0.0; XX for (k=n-1;k>0;k=k-1) { XX XX do i=1,nterms2 { XX kk=tagx(k,v) XX kk1=tagx(k+1,v) XX scr1(i)=scr1(i)+(bxorth(kk1,i)-bxorthm(i))*bx(kk1,m) XX cov(i,nterms21)=cov(i,nterms21)+ (x(kk1,v)-x(kk,v))*scr1(i) XX cov(nterms21,i)=cov(i,nterms21) XX } XX scr2=scr2+(bx(kk1,m)**2)*x(kk1,v) XX XX sumbx2=sumbx2+bx(kk1,m)**2 XX sumb=sumb+bx(kk1,m) XX sumbx=sumbx+bx(kk1,m)*x(kk1,v) XX su=st XX st=sumbx-sumb*x(kk,v) XX cov(nterms21,nterms21)= cov(nterms21,nterms21)+ XX (x(kk1,v)-x(kk,v))*(2*scr2-sumbx2*(x(kk,v)+x(kk1,v)))+ XX ( (su*su)-(st*st) )/n XX XX XX crittemp=critnew XX do kc=1, nclass{ XX XX scr6(kc)=scr6(kc)+(y(kk1,kc)-ybar(kc))*bx(kk1,m) XX covsy(nterms21,kc)=covsy(nterms21,kc )+(x(kk1,v)-x(kk,v))*scr6(kc) XX XX temp1=covsy(nterms21,kc) XX temp2=cov(nterms21,nterms21) XX XX do jk=1,nterms2 { XX temp1=temp1-covsy(jk,kc)*cov(jk,nterms21) XX temp2=temp2-cov(jk,nterms21)*cov(jk,nterms21) XX } # # this has to be fixed!! # XX XX XX XX if(cov(nterms21,nterms21)>0){ XX if((temp2/cov(nterms21,nterms21)) > tolbx) {critadd=(temp1*temp1)/temp2} else {critadd=0.0}} XX else{critadd=0} XX XX crittemp=crittemp+critadd XX #if((nterms.eq.7).and.(v.eq.2).and.(critadd.gt.0)) { XX #call intpr("k",1,k,1) XX #call dblepr("temp2",5,temp2,1) #call dblepr("cov",3,cov(nterms21,nterms21),1) # call dblepr("critadd",7,critadd,1) # call dblepr("crittemp",8,crittemp,1) XX #} XX XX XX XX XX if(crittemp.gt.(1.01*rss)) {crittemp=0.0} XX if(crittemp.gt.(2*prevcrit)) {crittemp=0.0} XX XX } XX XX XX XX if(k.gt.1){ XX k0=tagx(k-1,v) } # # decide whether to accept term # XX #if((k.eq.70).and.(v.eq.10)){ ##write(6,*) "HI",k,v,critmax,minspan,mod(k,minspan),iendspan,n-iendspan,x(kk,v),x(kk0,v) #} #call intpr("k",1,k,1) #call intpr("v",1,v,1) #call dblepr("crittemp",8,crittemp,1) #call intpr("minspan",7,minspan,1) #call intpr("iendspan",8,iendspan,1) XX if((crittemp.gt.critmax).and.(mod(k,minspan).eq.0).and. (k.ge.iendspan).and.(k.le.(n-iendspan)).and.(bx(kk1,m).gt.0).and. (.not.((k.gt.1).and.(x(kk,v).eq.x(k0,v))) )) { XX jmax=v XX critmax=crittemp XX imax=k XX kmax=m XX newform=tnewform XX XX XX } XX XX XX XX } } 9999 continue XX }} XX XX } XX XX XX return end XX XX XX XX EEEEE X $shar_touch -am 0405184995 'dmarss.r' && X chmod 0644 'dmarss.r' || X $echo 'restore of' 'dmarss.r' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo 'dmarss.r:' 'MD5 check failed' c1a19ba270625cbef0206481a102893e dmarss.r XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'dmarss.r'`" X test 20070 -eq "$shar_count" || X $echo 'dmarss.r:' 'original size' '20070,' 'current size' "$shar_count!" X fi fi # ============= dorthreg.r ============== if test -f 'dorthreg.r' && test "$first_param" != -c; then X $echo 'x -' SKIPPING 'dorthreg.r' '(file already exists)' else X $echo 'x -' extracting 'dorthreg.r' '(text)' X sed 's/^X//' << 'EEEEE' > 'dorthreg.r' && subroutine orthreg(nx,n,p,x,in, y,res) # # does lin reg of y on x. assumes that x is orthogonal with cols > 1 having mean 0 # # # "in" is a vector of column indicators (0 means term is to be deleted) XX implicit double precision (a-h,o-z) XX XX integer n,nx,p, in(p) double precision x(nx,p),y(n),res(n) do i=1,n { XX res(i)=y(i) } do j=1,p { XX if(in(j).eq.1){ XX temp1=0 XX temp2=0 XX do i=1,n { XX temp1=temp1+res(i)*x(i,j) XX temp2=temp2+x(i,j)*x(i,j) XX } XX beta=temp1/temp2 XX do i=1,n { XX res(i)=res(i)-beta*x(i,j) XX } }} return end XX XX EEEEE X $shar_touch -am 0405184995 'dorthreg.r' && X chmod 0644 'dorthreg.r' || X $echo 'restore of' 'dorthreg.r' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo 'dorthreg.r:' 'MD5 check failed' 5468af4a41b4e3653bb8732e00eeeb1a dorthreg.r XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'dorthreg.r'`" X test 527 -eq "$shar_count" || X $echo 'dorthreg.r:' 'original size' '527,' 'current size' "$shar_count!" X fi fi # ============= dqrreg.r ============== if test -f 'dqrreg.r' && test "$first_param" != -c; then X $echo 'x -' SKIPPING 'dqrreg.r' '(file already exists)' else X $echo 'x -' extracting 'dqrreg.r' '(text)' X sed 's/^X//' << 'EEEEE' > 'dqrreg.r' && subroutine qrreg(nx,n,px,p,nclass,x,xsc,in,y,qpivot,qrank,beta,res,rss,cvar,var,varsc,scr1,work) XX implicit double precision (a-h,o-z) XX integer nx,n,p,px, qpivot(p),qrank,nclass,in(p) double precision x(nx,p), xsc(n,p), y(n,nclass),res(nx,nclass),beta(px,nclass),work(1),scr1(p),var(px,p),varsc(px,p) logical cvar XX ii=0 do j=1,p { XX if(in(j).eq.1){ XX ii=ii+1 XX do i=1,n { XX xsc(i,ii)=x(i,j) XX } }} nt=ii ijob=101 info=1 temp3=1d-2 do i=1,p {qpivot(i)=i} call dqrdca(xsc,n,n,nt,scr1,qpivot,work,qrank,temp3) # computes both fits and beta rss=0.0 XX do k=1,nclass{ XX call dqrsl(xsc,n,n,qrank,scr1,y(1,k),work(1),work(1),beta(1,k), XX work(1),res(1,k),ijob,info) XX do i=1,n { res(i,k)=y(i,k)-res(i,k); rss=rss+res(i,k)*res(i,k)} XX } XX if(cvar) {call calcvar(nx,n,px,xsc,qrank,qpivot,var,varsc,work)} return end EEEEE X $shar_touch -am 0405184995 'dqrreg.r' && X chmod 0644 'dqrreg.r' || X $echo 'restore of' 'dqrreg.r' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo 'dqrreg.r:' 'MD5 check failed' 87e00acb8f53e360d0e7b19403cc0360 dqrreg.r XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'dqrreg.r'`" X test 807 -eq "$shar_count" || X $echo 'dqrreg.r:' 'original size' '807,' 'current size' "$shar_count!" X fi fi # ============= mspline.r ============== if test -f 'mspline.r' && test "$first_param" != -c; then X $echo 'x -' SKIPPING 'mspline.r' '(file already exists)' else X $echo 'x -' extracting 'mspline.r' '(text)' X sed 's/^X//' << 'EEEEE' > 'mspline.r' && subroutine sspl(x,y,w,n,ldy,p,knot,nk,method,tol,wp,ssy, XX dfoff,dfmax,cost,lambda,df,cv,gcv,coef,s,lev, XX xwy,hs,sg,abd,p1ip,ier) implicit double precision(a-h,o-z) XX # A Cubic B-spline Smoothing routine. # # The algorithm minimises: # # (1/n) * sum w(i)* (y(i)-s(i))**2 + lambda* int ( s"(x) )**2 dx # # for each of p response variables in y #Input XX # x(n) vector containing the ordinates of the observations # y(ldy,p) matrix (n x p) of responses (ldy can be greater than n) # w(n) vector containing the weights given to each data point # n number of data points # ldy leading dimension of y # p number of columns in y # knot(nk+4) vector of knot points defining the cubic b-spline basis. # nk number of b-spline coefficients to be estimated # nk <= n+2 # method method for selecting amount of smoothing, lambda # 1 = fixed lambda # 2 = fixed df # 3 = gcv # 4 = cv # tol used in Golden Search routine # wp(p) weights, length p, used to combine cv or gcv in 3 or 4 above # ssy(p) offsets for weighted sum of squares for y; can be all zero, # else should be the variability lost due to collapsing # onto unique values # dfoff offset df used in gcv calculations (0 is good default) # dfmax maximum value for df allowed when gcv or cv are used # routine simply returns the value at dfmax if it was exceeded # cost cost per df (1 is good default) #Input/Output # lambda penalised likelihood smoothing parameter # df trace(S) #Output # cv omnibus cv criterion # gcv omnibus gcv criterion (including penalty and offset) # coef(nk,p) vector of spline coefficients # s(ldy,p) matrix of smoothed y-values # lev(n) vector of leverages # Working arrays/matrix # xwy(nk,p) X'Wy # hs(nk,4) the diagonals of the X'WX matrix # sg(nk,4) the diagonals of the Gram matrix # abd(4,nk) [ X'WX+lambda*SIGMA] in diagonal form # p1ip(4,nk) inner products between columns of L inverse # ier error indicator # ier = 0 ___ everything fine # ier = 1 ___ spar too small or too big # problem in cholesky decomposition double precision x(n),y(ldy,p),w(n),knot(nk+4),tol,wp(p),ssy(p), XX dfoff,dfmax,cost,lambda,df,cv,gcv,coef(nk,p),s(ldy,p),lev(n), XX xwy(nk,p),hs(nk,4),sg(nk,4),abd(4,nk),p1ip(4,nk) integer n,p,ldy,nk,method,ier XX # Compute SIGMA, X'WX, X'WY, trace, ratio, s0, s1. XX # SIGMA-> sg[] # X'WX -> hs[] # X'WY -> xwy[] call sgram(sg(1,1),sg(1,2),sg(1,3),sg(1,4),knot,nk) call stxwx2(x,y,w,n,ldy,p,knot,nk,xwy,hs(1,1),hs(1,2),hs(1,3),hs(1,4)) XX # Compute estimate XX if(method==1) {# Value of lambda supplied XX XX call sslvr2(x,y,w,n,ldy,p,knot,nk,method,tol,wp,ssy, XX dfoff,cost,lambda,df,cv,gcv,coef,s,lev, XX xwy,hs(1,1),hs(1,2),hs(1,3),hs(1,4), XX sg(1,1),sg(1,2),sg(1,3),sg(1,4),abd,p1ip,ier) XX } XX else {# Use Forsythe, Malcom and Moler routine to minimise criterion XX call fmm(x,y,w,n,ldy,p,knot,nk,method,tol,wp,ssy, XX dfoff,cost,lambda,df,cv,gcv,coef,s,lev, XX xwy,hs,sg,abd,p1ip,ier) XX if(method>2&df>dfmax){ XX df=dfmax XX call fmm(x,y,w,n,ldy,p,knot,nk,2,tol,wp,ssy, XX dfoff,cost,lambda,df,cv,gcv,coef,s,lev, XX xwy,hs,sg,abd,p1ip,ier) XX } XX } return end XX subroutine fmm(xs,ys,ws,n,ldy,nvar,knot,nk,method,tol,wp,ssy, XX dfoff,cost,lambda,df,cv,gcv,coef,s,lev, XX xwy,hs,sg,abd,p1ip,ier) double precision xs(n),ys(ldy,nvar),ws(n),knot(nk+4),tol,wp(nvar),ssy(nvar), XX dfoff,cost,lambda,df,cv,gcv,coef(nk,nvar),s(ldy,nvar),lev(n), XX xwy(nk,nvar),hs(nk,4),sg(nk,4),abd(4,nk),p1ip(4,nk) integer n,ldy,nvar,nk,method,ier XX # Local variables XX double precision t1,t2,ratio, XX a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w, XX fu,fv,fw,fx,x,targdf, XX ax,bx XX integer i ax=1d-10 #used to be lspar bx=1.5 #used to be uspar t1=0. ; t2=0. targdf=df do i=3,nk-3 { t1 = t1 + hs(i,1) } do i=3,nk-3 { t2 = t2 + sg(i,1) } ratio = t1/t2 # # an approximation x to the point where f attains a minimum on # the interval (ax,bx) is determined. # # # input.. # # ax left endpoint of initial interval # bx right endpoint of initial interval # f function subprogram which evaluates f(x) for any x # in the interval (ax,bx) # tol desired length of the interval of uncertainty of the final # result ( .ge. 0.0) # # # output.. # # fmin abcissa approximating the point where f attains a minimum # # # the method used is a combination of golden section search and # successive parabolic interpolation. convergence is never much slower # than that for a fibonacci search. if f has a continuous second # derivative which is positive at the minimum (which is not at ax or # bx), then convergence is superlinear, and usually of the order of # about 1.324.... # the function f is never evaluated at two points closer together # than eps*dabs(fmin) + (tol/3), where eps is approximately the square # root of the relative machine precision. if f is a unimodal # function and the computed values of f are always unimodal when # separated by at least eps*dabs(x) + (tol/3), then fmin approximates # the abcissa of the global minimum of f on the interval ax,bx with # an error less than 3*eps*dabs(fmin) + tol. if f is not unimodal, # then fmin may approximate a local, but perhaps non-global, minimum to # the same accuracy. # this function subprogram is a slightly modified version of the # algol 60 procedure localmin given in richard brent, algorithms for # minimization without derivatives, prentice - hall, inc. (1973). # # # double precision a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w # double precision fu,fv,fw,fx,x # # c is the squared inverse of the golden ratio # XX c = 0.5*(3. - dsqrt(5d0)) # # eps is approximately the square root of the relative machine # precision. # XX eps = 1d0 XX 10 eps = eps/2d0 XX tol1 = 1d0 + eps XX if (tol1 .gt. 1d0) go to 10 XX eps = dsqrt(eps) # # initialization # XX a = ax XX b = bx XX v = a + c*(b - a) XX w = v XX x = v XX e = 0.0 XX lambda = ratio*16.**(-2. + x*(6.)) XX call sslvr2(xs,ys,ws,n,ldy,nvar,knot,nk,method,tol,wp,ssy, XX dfoff,cost,lambda,df,cv,gcv,coef,s,lev, XX xwy,hs(1,1),hs(1,2),hs(1,3),hs(1,4), XX sg(1,1),sg(1,2),sg(1,3),sg(1,4),abd,p1ip,ier) switch(method){ XX case 2: XX fx=3d0+(targdf-df)**2 XX case 3: XX fx=gcv XX case 4: XX fx=cv XX } XX fv = fx XX fw = fx # # main loop starts here # XX 20 xm = 0.5*(a + b) XX tol1 = eps*dabs(x) + tol/3d0 XX tol2 = 2d0*tol1 # # check stopping criterion # XX if (dabs(x - xm) .le. (tol2 - 0.5*(b - a))) go to 90 # # is golden-section necessary # XX if (dabs(e) .le. tol1) go to 40 # # fit parabola # XX r = (x - w)*(fx - fv) XX q = (x - v)*(fx - fw) XX p = (x - v)*q - (x - w)*r XX q = 2.00*(q - r) XX if (q .gt. 0.0) p = -p XX q = dabs(q) XX r = e XX e = d # # is parabola acceptable # XX 30 if (dabs(p) .ge. dabs(0.5*q*r)) go to 40 XX if (p .le. q*(a - x)) go to 40 XX if (p .ge. q*(b - x)) go to 40 # # a parabolic interpolation step # XX d = p/q XX u = x + d # # f must not be evaluated too close to ax or bx # XX if ((u - a) .lt. tol2) d = dsign(tol1, xm - x) XX if ((b - u) .lt. tol2) d = dsign(tol1, xm - x) XX go to 50 # # a golden-section step # XX 40 if (x .ge. xm) e = a - x XX if (x .lt. xm) e = b - x XX d = c*e # # f must not be evaluated too close to x # XX 50 if (dabs(d) .ge. tol1) u = x + d XX if (dabs(d) .lt. tol1) u = x + dsign(tol1, d) XX XX lambda = ratio*16.**(-2. + u*(6.)) XX call sslvr2(xs,ys,ws,n,ldy,nvar,knot,nk,method,tol,wp,ssy, XX dfoff,cost,lambda,df,cv,gcv,coef,s,lev, XX xwy,hs(1,1),hs(1,2),hs(1,3),hs(1,4), XX sg(1,1),sg(1,2),sg(1,3),sg(1,4),abd,p1ip,ier) switch(method){ XX case 2: XX fu=3d0+(targdf-df)**2 XX case 3: XX fu=gcv XX case 4: XX fu=cv XX } # # update a, b, v, w, and x # XX if (fu .gt. fx) go to 60 XX if (u .ge. x) a = x XX if (u .lt. x) b = x XX v = w XX fv = fw XX w = x XX fw = fx XX x = u XX fx = fu XX go to 20 XX 60 if (u .lt. x) a = u XX if (u .ge. x) b = u XX if (fu .le. fw) go to 70 XX if (w .eq. x) go to 70 XX if (fu .le. fv) go to 80 XX if (v .eq. x) go to 80 XX if (v .eq. w) go to 80 XX go to 20 XX 70 v = w XX fv = fw XX w = u XX fw = fu XX go to 20 XX 80 v = u XX fv = fu XX go to 20 # # end of main loop # XX 90 continue if(method==2){call sslvr2(xs,ys,ws,n,ldy,nvar,knot,nk,1,tol,wp,ssy, XX dfoff,cost,lambda,df,cv,gcv,coef,s,lev, XX xwy,hs(1,1),hs(1,2),hs(1,3),hs(1,4), XX sg(1,1),sg(1,2),sg(1,3),sg(1,4),abd,p1ip,ier)} XX return end subroutine stxwx2(x,z,w,k,ldy,pz,xknot,n,y,hs0,hs1,hs2,hs3) implicit double precision(a-h,o-z) double precision z(ldy,pz),w(k),x(k),xknot(n+4), XX y(n,pz),hs0(n),hs1(n),hs2(n),hs3(n), XX eps,vnikx(4,1),work(16) # local integer k,n,pz,ldy, j,i,pp,ileft,mflag # local XX # Initialise the output vectors XX do i=1,n { XX hs0(i)=0d0 XX hs1(i)=0d0 XX hs2(i)=0d0 XX hs3(i)=0d0 XX do j=1,pz { y(i,j)=0d0 } XX } XX # Compute X'WX -> hs0,hs1,hs2,hs3 and X'WZ -> y XX eps = .1d-9 do i=1,k { XX call interv(xknot(1),(n+1),x(i),ileft,mflag) XX if(mflag== 1) { XX if(x(i)<=(xknot(ileft)+eps)){ileft=ileft-1} XX else{return} XX } XX call bsplvd (xknot,4,x(i),ileft,work,vnikx,1) XX XX j= ileft-4+1 XX do pp=1,pz {y(j,pp) = y(j,pp)+w(i)*z(i,pp)*vnikx(1,1)} XX hs0(j)=hs0(j)+w(i)*vnikx(1,1)**2 XX hs1(j)=hs1(j)+w(i)*vnikx(1,1)*vnikx(2,1) XX hs2(j)=hs2(j)+w(i)*vnikx(1,1)*vnikx(3,1) XX hs3(j)=hs3(j)+w(i)*vnikx(1,1)*vnikx(4,1) XX XX j= ileft-4+2 XX do pp=1,pz {y(j,pp) = y(j,pp)+w(i)*z(i,pp)*vnikx(2,1)} XX hs0(j)=hs0(j)+w(i)*vnikx(2,1)**2 XX hs1(j)=hs1(j)+w(i)*vnikx(2,1)*vnikx(3,1) XX hs2(j)=hs2(j)+w(i)*vnikx(2,1)*vnikx(4,1) XX XX j= ileft-4+3 XX do pp=1,pz {y(j,pp) = y(j,pp)+w(i)*z(i,pp)*vnikx(3,1)} XX hs0(j)=hs0(j)+w(i)*vnikx(3,1)**2 XX hs1(j)=hs1(j)+w(i)*vnikx(3,1)*vnikx(4,1) XX XX j= ileft-4+4 XX do pp=1,pz {y(j,pp) = y(j,pp)+w(i)*z(i,pp)*vnikx(4,1)} XX hs0(j)=hs0(j)+w(i)*vnikx(4,1)**2 XX } XX return end subroutine sslvr2(x,y,w,n,ldy,p,knot,nk,method,tol,wp,ssy, XX dfoff,cost,lambda,df,cv,gcv,coef,sz,lev, XX xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, XX abd,p1ip,info) XX implicit double precision(a-h,o-z) double precision x(n),y(ldy,p),w(n),knot(nk+4),tol,wp(p),ssy(p), XX dfoff,cost,lambda,df,cv,gcv,coef(nk,p),sz(ldy,p),lev(n), XX xwy(nk,p), XX hs0(nk),hs1(nk),hs2(nk),hs3(nk), XX sg0(nk),sg1(nk),sg2(nk),sg3(nk), XX abd(4,nk),p1ip(4,nk) integer n,p,ldy,nk,method,info #local storage double precision b0,b1,b2,b3,eps,vnikx(4,1),work(16), XX xv,bvalue,rss,tssy integer ld4,i,icoef,ileft,ilo,j,mflag logical fittoo fittoo= (method!=2) XX XX ilo = 1 ; eps = .1d-10 ; ld4=4 XX # Purpose : Solves the smoothing problem and computes the # criterion functions (CV and GCV). # The coeficients of estimated smooth XX if(fittoo){ XX do i=1,nk { do j=1,p {coef(i,j) = xwy(i,j) } } XX } do i=1,nk { abd(4,i) = hs0(i)+lambda*sg0(i) } do i=1,(nk-1) { abd(3,i+1) = hs1(i)+lambda*sg1(i) } do i=1,(nk-2) { abd(2,i+2) = hs2(i)+lambda*sg2(i) } do i=1,(nk-3) { abd(1,i+3) = hs3(i)+lambda*sg3(i) } XX call dpbfa(abd,ld4,nk,3,info) if(info.ne.0) { return } if(fittoo){ XX do j=1,p{ call dpbsl(abd,ld4,nk,3,coef(1,j)) } XX # Value of smooths at the data points XX XX icoef = 1 XX do i=1,n { XX xv = x(i) XX do j=1,p{ sz(i,j) = bvalue(knot,coef(1,j),nk,4,xv,0) } XX } XX } # Compute the criterion functions XX # Get Leverages First XX call sinrp2(abd,ld4,nk,p1ip) XX do i=1,n { XX xv = x(i) XX call interv(knot(1),(nk+1),xv,ileft,mflag) XX if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } XX if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } XX j=ileft-3 XX call bsplvd(knot,4,xv,ileft,work,vnikx,1) XX b0=vnikx(1,1);b1=vnikx(2,1);b2=vnikx(3,1);b3=vnikx(4,1) XX lev(i) = (p1ip(4,j)*b0**2 + 2.*p1ip(3,j)*b0*b1 + XX 2.*p1ip(2,j)*b0*b2 + 2.*p1ip(1,j)*b0*b3 + XX p1ip(4,j+1)*b1**2 + 2.*p1ip(3,j+1)*b1*b2 + XX 2.*p1ip(2,j+1)*b1*b3 + XX p1ip(4,j+2)*b2**2 + 2.*p1ip(3,j+2)*b2*b3 + XX p1ip(4,j+3)*b3**2 )*w(i) XX } # Evaluate Criteria rss = 0d0 ; df = 0d0 ; sumw=0d0;gcv=0d0;cv=0d0; do i=1,n { XX df = df + lev(i) XX } if(fittoo){ XX do i=1,n { XX sumw = sumw + w(i) XX do j=1,p{ XX rss = rss + w(i)*wp(j)*(y(i,j)-sz(i,j))**2 XX cv = cv +w(i)*wp(j)*((y(i,j)-sz(i,j))/(1-lev(i)))**2 XX } XX } tssy=0d0 do j=1,p{tssy=tssy+wp(j)*ssy(j)} XX gcv=((rss+tssy)/sumw)/((1d0-((dfoff+df-1)*cost+1)/sumw)**2) XX #note: the weights should sum to n (the number of original observations) XX cv=(cv+tssy)/sumw XX } #lev includes the weights #Note that this version of cv omits ALL observations at #tied x values, since the data are already collapsed here return end subroutine sinrp2(abd,ld4,nk,p1ip) implicit double precision(a-h,o-z) double precision abd(ld4,nk),p1ip(ld4,nk), XX wjm3(3),wjm2(2),wjm1(1),c0,c1,c2,c3 XX integer ld4,nk,i,j,k XX XX # Purpose : Computes Inner Products between columns of L(-1) XX # where L = abd is a Banded Matrix with 3 subdiagonals XX XX # A refinement of Elden's trick is used. XX # Coded by Finbarr O'Sullivan wjm3(1)=0d0; wjm3(2)=0d0; wjm3(1)=0d0 wjm2(1)=0d0; wjm2(2)=0d0 wjm1(1)=0d0 XX do i=1,nk { XX j=nk-i+1 XX c0 = 1d0/abd(4,j) XX if(j<=nk-3) { XX c1 = abd(1,j+3)*c0 XX c2 = abd(2,j+2)*c0 XX c3 = abd(3,j+1)*c0 XX } XX else if(j==nk-2) { XX c1 = 0d0 XX c2 = abd(2,j+2)*c0 XX c3 = abd(3,j+1)*c0 XX } XX else if(j==nk-1) { XX c1 = 0d0 XX c2 = 0d0 XX c3 = abd(3,j+1)*c0 XX } XX else if(j==nk) { XX c1 = 0d0 XX c2 = 0d0 XX c3 = 0d0 XX } XX p1ip(1,j) = 0d0- (c1*wjm3(1)+c2*wjm3(2)+c3*wjm3(3)) XX p1ip(2,j) = 0d0- (c1*wjm3(2)+c2*wjm2(1)+c3*wjm2(2)) XX p1ip(3,j) = 0d0- (c1*wjm3(3)+c2*wjm2(2)+c3*wjm1(1)) XX XX p1ip(4,j) = c0**2 + XX c1**2*wjm3(1)+2.*c1*c2*wjm3(2)+2.*c1*c3*wjm3(3) + XX c2**2*wjm2(1)+2.*c2*c3*wjm2(2) + XX c3**2*wjm1(1) XX XX wjm3(1)=wjm2(1) ; wjm3(2)=wjm2(2) ; wjm3(3)=p1ip(2,j) XX wjm2(1)=wjm1(1) ; wjm2(2)=p1ip(3,j); wjm1(1)=p1ip(4,j) XX } return end subroutine suff2(n,p,ny,match,y,w,ybar,wbar,ssy,work) integer match(n),n,ny,p,i,j double precision y(n,ny),ybar(p+1,ny),w(n),wbar(p+1),ssy(ny),work(n) double precision tsum #ssy is the within response variability that is lost by collapsing call pack(n,p,match,w,wbar) do j=1,ny{ XX do i=1,n XX work(i)=y(i,j)*w(i) XX call pack(n,p,match,work,ybar(1,j)) XX do i=1,p{ XX if(wbar(i)>0d0) XX ybar(i,j)=ybar(i,j)/wbar(i) XX else ybar(i,j)=0d0 XX } XX call unpack(n,p,match,ybar(1,j),work) XX tsum=0d0 XX do i=1,n XX tsum=tsum+ w(i)*(y(i,j)-work(i))**2 XX ssy(j)=tsum XX } return end subroutine sspl0(x,y,w,n,p,knot,nk,method,tol,wp,match,nef,center, XX dfoff,dfmax,cost,lambda,df,cv,gcv,coef,s,lev, XX xrange,work,ier) double precision x(n),y(n,p),w(n),knot(nk+4),tol,wp(p), XX dfoff,dfmax,cost,lambda,df,cv,gcv,coef(1),s(n,p),lev(nef), XX xrange(2),work(1) #workspace must be (2*p+2)*nefp1 + (p+16)*nk + n +p integer n,p,nk,method,ier,nef, nefp1, n2,match(1) logical center double precision xmiss,sigtol,xdiff,wmean,temp if(nef==0){# match has not been initialized XX xmiss=1d20 XX sigtol=1d-5 XX call namat(x,match,n,nef,work,work(n+1),xmiss,sigtol) XX xrange(1)=work(1) #work is actually the sorted unique xs XX xrange(2)=work(nef) XX } else{ XX do i=1,n {work(match(i))=x(i)} XX } xdiff=xrange(2)-xrange(1) do i=1,nef {work(i)=(work(i)-xrange(1))/xdiff} if(nk==0){ XX call sknotl(work,nef,knot,nk) XX nk=nk-4 XX } if(dfmax > dble(nk))dfmax=dble(nk) if(cost>0){ XX temp=dble(n-center)/cost - dfoff XX if(dfmax>temp)dfmax=temp XX } nefp1=nef+1 n2=nefp1*(2*p+2)+1 call sspl1(x,y,w,n,p,knot,nk,method,tol,wp,match,nef,nefp1,center, XX dfoff,dfmax,cost,lambda,df,cv,gcv,coef,s,lev, XX work(1),work(nefp1+1), #xin,yin XX work(nefp1*(p+1)+1),work(nefp1*(p+2)+1), #win, sout XX work(n2),work(n2+p*nk),work(n2+(p+4)*nk), #xwy, hs,sg XX work(n2+(p+8)*nk),work(n2+(p+12)*nk),work(n2+(p+16)*nk), XX work(n2+(p+16)*nk+p),ier) XX return end XX #Memory management subroutine subroutine sspl1(x,y,w,n,p,knot,nk,method,tol,wp,match,nef,nefp1,center, XX dfoff,dfmax,cost,lambda,df,cv,gcv,coef,s,lev, XX xin,yin,win,sout, XX xwy,hs,sg,abd,p1ip,ssy,work,ier) double precision x(n),y(n,p),w(n),knot(nk+4),tol,wp(p), XX dfoff,dfmax,cost,lambda,df,cv,gcv,coef(nk,p),s(n,p),lev(nef), XX xin(nefp1),yin(nefp1,p),win(nefp1),sout(nefp1,p), XX xwy(nk,p),hs(nk,4),sg(nk,4),abd(4,nk),p1ip(4,nk), XX ssy(p),work(n) integer n,p,nefp1,nk,method,ier,match(n),nef logical center double precision xbar,sbar, wmean,tdfoff call suff2(n,nef,p,match,y,w,yin,win,ssy,work) if(center){ XX if(cost>0){tdfoff=dfoff-1/cost} } call sspl(xin,yin,win,nef,nefp1,p,knot,nk,method,tol,wp,ssy, XX tdfoff,dfmax,cost,lambda,df,cv,gcv,coef,sout,lev, XX xwy,hs,sg,abd,p1ip,ier) XX #now unpack the results do j=1,p{ XX call unpack(n,nef,match,sout(1,j),s(1,j)) XX if(center){ XX sbar=wmean(nef,sout(1,j),win) XX do i=1,n{ s(i,j)=s(i,j)-sbar} XX } } if(center)df=df-1 return end subroutine namat(x,match,n,nef,work,iwork,xmiss,tol) #returns match (order) and work(1:nef) is the sorted unique x implicit double precision(a-h,o-z) double precision x(1),xmiss,work(1),tol,xend,xstart integer match(1),n,nef,iwork(1),index do i=1,n { XX work(i)=x(i) XX iwork(i)=i XX } XX call sortdi(work,n,iwork,1,n) xstart=x(iwork(1)) index=n xend=x(iwork(n)) while(xend >= xmiss & index > 1){ XX index=index-1 XX xend=x(iwork(index)) XX } XX tol=tol*(xend-xstart) index=1 work(1)=xstart for(i=1;i<=n;i=i+1){ XX while((x(iwork(i))-xstart)n)goto 10 XX } XX xstart= x(iwork(i)) XX index=index+1 XX match(iwork(i))=index XX work(index)=xstart XX } 10 if(xstart >= xmiss) XX {nef=index-1} XX else {nef=index} return end XX subroutine simfit(x,y,w,n,p,dfoff,cost,wp,gcv,coef,s,type,center,work) # # computes constant and linear fits, and selects the best using gcv # implicit double precision (a-h,o-z) XX integer n,p,type double precision x(n),y(n,p),w(n),cost,dfoff,wp(p),gcv,coef(2,p), XX s(n,p),work(p) logical center double precision sx,sy,sumw, dcent,sxx,syy,sxy dcent=1-dble(center) #center is F for no centering, else T #Note: if type enters 1 or 2, no selection is made sumw=0d0;gcvc=0d0; do i=1,n { XX sumw=sumw+w(i) } if(type!=1){#either 0 or 2 in which case the linear is needed as well XX sx=0.0 ; sxx=0.0; gcvl=0d0; XX XX do i=1,n { XX sx=sx+w(i)*x(i) XX } XX xbar=sx/sumw XX do i=1,n { XX sxx=sxx+w(i)*(x(i)-xbar)*x(i) XX } } do j=1,p{ XX sy=0d0;syy=0d0; XX do i=1,n{ XX sy=sy+w(i)*y(i,j) XX } XX work(j)=sy/sumw XX do i=1,n{ XX syy=syy+w(i)*(y(i,j)-work(j))*y(i,j) XX } XX gcvc=gcvc+wp(j)*syy XX if(type!=1){ #once again, do for linear as well XX sxy=0.0; XX do i=1,n{ XX sxy=sxy+w(i)*(x(i)-xbar)*y(i,j) XX } XX coef(2,j)=sxy/sxx XX gcvl=gcvl+wp(j)*(syy -sxy*coef(2,j)) XX } } if(type==0){ XX gcvc =gcvc/ (sumw* (1 - (dfoff*cost + dcent)/sumw)**2 ) XX gcvl=gcvl/(sumw* (1 - (dcent + (dfoff +1)* cost)/sumw)**2) XX if(gcvc<=gcvl){ XX type=1 XX gcv=gcvc XX } XX else{ XX type=2 XX gcv=gcvl XX } } else { XX XX if(type==1) {gcv=gcvc/(sumw* (1 - (dfoff*cost + dcent)/sumw)**2 )} XX else {gcv=gcvl/(sumw* (1 - (dcent + (dfoff + 1)*cost)/sumw)**2)} } XX if(type==1){ XX do j=1,p{ XX coef(1,j)=work(j)*dcent XX coef(2,j)=0d0 XX do i=1,n {s(i,j)=coef(1,j)} XX } XX } else{ XX do j=1,p{ XX coef(1,j)=work(j)*dcent-xbar*coef(2,j) XX do i=1,n {s(i,j)=coef(1,j)+coef(2,j)*x(i)} XX } XX } return end XX XX subroutine sspl2(x,y,w,n,p,knot,nk,wp,match,nef, XX dfoff,dfmax,cost,lambda,df,gcv,coef,s,type,center, XX xrange,work,tol,ier) double precision x(n),y(n,p),w(n),knot(nk+4),wp(p), XX dfoff,dfmax,cost,lambda,df,gcv,coef(1),s(n,p), XX xrange(2),work(1),tol #this routine selects from the linear and constant model as well #see documentation for sspl #workspace must be (2*p+2)*nefp1 + (p+16)*nk + 2*n #if type>0 then no selection is performed; the fit is simply computed. integer n,p,nk,nef,type,match(n),ier,method double precision wmean, coef1,coef2,cv logical center #center is F for no centering, else T if(type==3){ XX method=1 call sspl0(x,y,w,n,p,knot,nk,method,tol,wp,match,nef,center, XX dfoff,dfmax,cost,lambda,df,cv,gcv,coef,s,work(1), XX xrange,work(n+1),ier) XX return } if(type>0){ XX call simfit(x,y,w,n,p,dfoff,cost,wp,gcv,coef,s, XX type,center,work) XX df=dble(type)- dble(center) return } #selection is being performed method=3 call sspl0(x,y,w,n,p,knot,nk,method,tol,wp,match,nef,center, XX dfoff,dfmax,cost,lambda,df,cv,gcv,coef,s,work(1), XX xrange,work(n+1),ier) gcv1=gcv call simfit(x,y,w,n,p,dfoff,cost,wp,gcv,work,work(2*p+1),type,center, XX work((n+2)*p+1)) if(gcv<=gcv1){ #the coef swapping is so as not to destroy the spline coefs if needed XX df=dble(type)- dble(center) XX do j=1,p{ XX coef1=work(1+(j-1)*2) XX coef2=work(2+(j-1)*2) XX if(type==1){ XX do i=1,n {s(i,j)=coef1} XX } XX else{ XX do i=1,n {s(i,j) =coef1+coef2*x(i)} XX } XX coef(1+(j-1)*2)=coef1 XX coef(2+(j-1)*2)=coef2 XX } XX } else{ XX type=3 XX gcv=gcv1 XX } return end subroutine psspl2(x,n,p,knot,nk,xrange,coef,coefl,s,order,type) implicit double precision(a-h,o-z) #make predictions from a fitted smoothing spline double precision x(n),knot(nk+4),xrange(2),coef(nk,p),coefl(2,1),s(n,p) integer n,p,nk,order, type double precision ytemp switch(type){ XX case 1:{ XX do j=1,p{ XX if(order>=1){ytemp=0d0} else {ytemp=coefl(1,j)} XX do i=1,n {s(i,j)=ytemp} XX } XX } XX case 2:{ XX if(order>=1){ XX do j=1,p{ XX if(order==1){ytemp=coefl(2,j)} XX else {ytemp=0d0} XX do i =1,n {s(i,j)=ytemp} XX } XX } XX else{ XX do j=1,p{ XX do i=1,n {s(i,j)=coefl(1,j)+coefl(2,j)*x(i)} XX } XX } XX } XX case 3: { XX call psspl(x,n,p,knot,nk,xrange,coef,s,order) XX } } return end XX XX XX subroutine psspl(x,n,p,knot,nk,xrange,coef,s,order) implicit double precision(a-h,o-z) #make predictions from a fitted smoothing spline, linear or constant double precision x(n),knot(nk+4),xrange(2),coef(nk,p),s(n,p) integer n,p,nk,order double precision xcs,xmin,xdif, endv(2),ends(2),xends(2),stemp double precision bvalue integer where if(order>2|order<0)return xdif=xrange(2)-xrange(1) xmin=xrange(1) xends(1)=0d0 xends(2)=1d0 do j=1,p{ XX if(order==0){ XX endv(1)=bvalue(knot,coef(1,j),nk,4,0d0,0) XX endv(2)=bvalue(knot,coef(1,j),nk,4,1d0,0) XX } XX if(order<=1){ XX ends(1)=bvalue(knot,coef(1,j),nk,4,0d0,1) XX ends(2)=bvalue(knot,coef(1,j),nk,4,1d0,1) XX } XX do i=1,n{ XX xcs=(x(i)-xmin)/xdif XX where=0 XX if(xcs<0d0){where=1} XX if(xcs>1d0){where=2} XX if(where>0){#beyond extreme knots XX switch(order){ XX case 0: XX stemp=endv(where)+ XX (xcs-xends(where))*ends(where) XX case 1: XX stemp=ends(where) XX case 2: XX stemp=0d0 XX } XX } XX else {stemp=bvalue(knot,coef(1,j),nk,4,xcs,order)} XX if(order>0){s(i,j)=stemp/(xdif**dble(order))} XX else s(i,j)=stemp XX } } return end EEEEE X $shar_touch -am 0405184995 'mspline.r' && X chmod 0644 'mspline.r' || X $echo 'restore of' 'mspline.r' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo 'mspline.r:' 'MD5 check failed' 2cb704627325910bac18e072a5bf87a1 mspline.r XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'mspline.r'`" X test 22241 -eq "$shar_count" || X $echo 'mspline.r:' 'original size' '22241,' 'current size' "$shar_count!" X fi fi # ============= sortdi.f ============== if test -f 'sortdi.f' && test "$first_param" != -c; then X $echo 'x -' SKIPPING 'sortdi.f' '(file already exists)' else X $echo 'x -' extracting 'sortdi.f' '(text)' X sed 's/^X//' << 'EEEEE' > 'sortdi.f' && XX subroutine sortdi (v,lenv,a,ii,jj) c c puts into a the permutation vector which sorts v into c increasing order. only elements from ii to jj are considered. c arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements c v is returned sorted c this is a modification of cacm algorithm #347 by r. c. singleton, c which is a modified hoare quicksort. c XX implicit double precision (a-h, o-z) XX integer lenv XX integer t,tt,ii,jj,iu(20),il(20) XX integer a(jj) XX double precision v(lenv), vt, vtt XX m=1 XX i=ii XX j=jj XX 10 if (i.ge.j) go to 80 XX 20 k=i XX ij=(j+i)/2 XX t=a(ij) XX vt=v(ij) XX if (v(i).le.vt) go to 30 XX a(ij)=a(i) XX a(i)=t XX t=a(ij) XX v(ij)=v(i) XX v(i)=vt XX vt=v(ij) XX 30 l=j XX if (v(j).ge.vt) go to 50 XX a(ij)=a(j) XX a(j)=t XX t=a(ij) XX v(ij)=v(j) XX v(j)=vt XX vt=v(ij) XX if (v(i).le.vt) go to 50 XX a(ij)=a(i) XX a(i)=t XX t=a(ij) XX v(ij)=v(i) XX v(i)=vt XX vt=v(ij) XX go to 50 XX 40 a(l)=a(k) XX a(k)=tt XX v(l)=v(k) XX v(k)=vtt XX 50 l=l-1 XX if (v(l).gt.vt) go to 50 XX tt=a(l) XX vtt=v(l) XX 60 k=k+1 XX if (v(k).lt.vt) go to 60 XX if (k.le.l) go to 40 XX if (l-i.le.j-k) go to 70 XX il(m)=i XX iu(m)=l XX i=k XX m=m+1 XX go to 90 XX 70 il(m)=k XX iu(m)=j XX j=l XX m=m+1 XX go to 90 XX 80 m=m-1 XX if (m.eq.0) return XX i=il(m) XX j=iu(m) XX 90 if (j-i.gt.10) go to 20 XX if (i.eq.ii) go to 10 XX i=i-1 XX 100 i=i+1 XX if (i.eq.j) go to 80 XX t=a(i+1) XX vt=v(i+1) XX if (v(i).le.vt) go to 100 XX k=i XX 110 a(k+1)=a(k) XX v(k+1)=v(k) XX k=k-1 XX if (vt.lt.v(k)) go to 110 XX a(k+1)=t XX v(k+1)=vt XX go to 100 XX end XX subroutine interv ( xt, lxt, x, left, mflag ) XX implicit double precision(a-h,o-z) computes left = max( i ; 1 .le. i .le. lxt .and. xt(i) .le. x ) . c c****** i n p u t ****** c xt.....a double precision sequence, of length lxt , assumed to be nondecreasing c lxt.....number of terms in the sequence xt . c x.....the point whose location with respect to the sequence xt is c to be determined. c c****** o u t p u t ****** c left, mflag.....both integers, whose value is c c 1 -1 if x .lt. xt(1) c i 0 if xt(i) .le. x .lt. xt(i+1) c lxt 1 if xt(lxt) .le. x c c in particular, mflag = 0 is the 'usual' case. mflag .ne. 0 c indicates that x lies outside the halfopen interval c xt(1) .le. y .lt. xt(lxt) . the asymmetric treatment of the c interval is due to the decision to make all pp functions cont- c inuous from the right. c c****** m e t h o d ****** c the program is designed to be efficient in the common situation that c it is called repeatedly, with x taken from an increasing or decrea- c sing sequence. this will happen, e.g., when a pp function is to be c graphed. the first guess for left is therefore taken to be the val- c ue returned at the previous call and stored in the l o c a l varia- c ble ilo . a first check ascertains that ilo .lt. lxt (this is nec- c essary since the present call may have nothing to do with the previ- c ous call). then, if xt(ilo) .le. x .lt. xt(ilo+1), we set left = c ilo and are done after just three comparisons. c otherwise, we repeatedly double the difference istep = ihi - ilo c while also moving ilo and ihi in the direction of x , until c xt(ilo) .le. x .lt. xt(ihi) , c after which we use bisection to get, in addition, ilo+1 = ihi . c left = ilo is then returned. c XX integer left,lxt,mflag, ihi,ilo,istep,middle XX double precision x,xt(lxt) XX data ilo /1/ c save ilo (a valid fortran statement in the new 1977 standard) XX ihi = ilo + 1 XX if (ihi .lt. lxt) go to 20 XX if (x .ge. xt(lxt)) go to 110 XX if (lxt .le. 1) go to 90 XX ilo = lxt - 1 XX ihi = lxt c XX 20 if (x .ge. xt(ihi)) go to 40 XX if (x .ge. xt(ilo)) go to 100 c c **** now x .lt. xt(ilo) . decrease ilo to capture x . XX 30 istep = 1 XX 31 ihi = ilo XX ilo = ihi - istep XX if (ilo .le. 1) go to 35 XX if (x .ge. xt(ilo)) go to 50 XX istep = istep*2 XX go to 31 XX 35 ilo = 1 XX if (x .lt. xt(1)) go to 90 XX go to 50 c **** now x .ge. xt(ihi) . increase ihi to capture x . XX 40 istep = 1 XX 41 ilo = ihi XX ihi = ilo + istep XX if (ihi .ge. lxt) go to 45 XX if (x .lt. xt(ihi)) go to 50 XX istep = istep*2 XX go to 41 XX 45 if (x .ge. xt(lxt)) go to 110 XX ihi = lxt c c **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval. XX 50 middle = (ilo + ihi)/2 XX if (middle .eq. ilo) go to 100 c note. it is assumed that middle = ilo in case ihi = ilo+1 . XX if (x .lt. xt(middle)) go to 53 XX ilo = middle XX go to 50 XX 53 ihi = middle XX go to 50 c**** set output and return. XX 90 mflag = -1 XX left = 1 XX return XX 100 mflag = 0 XX left = ilo XX return XX 110 mflag = 1 XX left = lxt XX return XX end XX subroutine bsplvd ( t, k, x, left, a, dbiatx, nderiv ) XX implicit double precision(a-h,o-z) calls bsplvb calculates value and deriv.s of all b-splines which do not vanish at x c c****** i n p u t ****** c t the knot array, of length left+k (at least) c k the order of the b-splines to be evaluated c x the point at which these values are sought c left an integer indicating the left endpoint of the interval of c interest. the k b-splines whose support contains the interval c (t(left), t(left+1)) c are to be considered. c a s s u m p t i o n - - - it is assumed that c t(left) .lt. t(left+1) c division by zero will result otherwise (in b s p l v b ). c also, the output is as advertised only if c t(left) .le. x .le. t(left+1) . c nderiv an integer indicating that values of b-splines and their c derivatives up to but not including the nderiv-th are asked c for. ( nderiv is replaced internally by the integer in (1,k) c closest to it.) c c****** w o r k a r e a ****** c a an array of order (k,k), to contain b-coeff.s of the derivat- c ives of a certain order of the k b-splines of interest. c c****** o u t p u t ****** c dbiatx an array of order (k,nderiv). its entry (i,m) contains c value of (m-1)st derivative of (left-k+i)-th b-spline of c order k for knot sequence t , i=m,...,k; m=1,...,nderiv. c c****** m e t h o d ****** c values at x of all the relevant b-splines of order k,k-1,..., c k+1-nderiv are generated via bsplvb and stored temporarily c in dbiatx . then, the b-coeffs of the required derivatives of the c b-splines of interest are generated by differencing, each from the c preceding one of lower order, and combined with the values of b- c splines of corresponding order in dbiatx to produce the desired c values. c XX integer k,left,nderiv, i,ideriv,il,j,jlow,jp1mid,kp1,kp1mm XX * ,ldummy,m,mhigh XX double precision a(k,k),dbiatx(k,nderiv),t(1),x XX double precision factor,fkp1mm,sum XX mhigh = max0(min0(nderiv,k),1) c mhigh is usually equal to nderiv. XX kp1 = k+1 XX call bsplvb(t,kp1-mhigh,1,x,left,dbiatx) XX if (mhigh .eq. 1) go to 99 c the first column of dbiatx always contains the b-spline values c for the current order. these are stored in column k+1-current c order before bsplvb is called to put values for the next c higher order on top of it. XX ideriv = mhigh XX do 15 m=2,mhigh XX jp1mid = 1 XX do 11 j=ideriv,k XX dbiatx(j,ideriv) = dbiatx(jp1mid,1) XX 11 jp1mid = jp1mid + 1 XX ideriv = ideriv - 1 XX call bsplvb(t,kp1-ideriv,2,x,left,dbiatx) XX 15 continue c c at this point, b(left-k+i, k+1-j)(x) is in dbiatx(i,j) for c i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the c first column of dbiatx is already in final form. to obtain cor- c responding derivatives of b-splines in subsequent columns, gene- c rate their b-repr. by differencing, then evaluate at x. c XX jlow = 1 XX do 20 i=1,k XX do 19 j=jlow,k XX 19 a(j,i) = 0d0 XX jlow = i XX 20 a(i,i) = 1d0 c at this point, a(.,j) contains the b-coeffs for the j-th of the c k b-splines of interest here. c XX do 40 m=2,mhigh XX kp1mm = kp1 - m XX fkp1mm = dble(kp1mm) XX il = left XX i = k c c for j=1,...,k, construct b-coeffs of (m-1)st derivative of c b-splines from those for preceding derivative by differencing c and store again in a(.,j) . the fact that a(i,j) = 0 for c i .lt. j is used.sed. XX do 25 ldummy=1,kp1mm XX factor = fkp1mm/(t(il+kp1mm) - t(il)) c the assumption that t(left).lt.t(left+1) makes denominator c in factor nonzero. XX do 24 j=1,i XX 24 a(i,j) = (a(i,j) - a(i-1,j))*factor XX il = il - 1 XX 25 i = i - 1 c c for i=1,...,k, combine b-coeffs a(.,i) with b-spline values c stored in dbiatx(.,m) to get value of (m-1)st derivative of c i-th b-spline (of interest here) at x , and store in c dbiatx(i,m). storage of this value over the value of a b-spline c of order m there is safe since the remaining b-spline derivat- c ive of the same order do not use this value due to the fact c that a(j,i) = 0 for j .lt. i . XX 30 do 40 i=1,k XX sum = 0. XX jlow = max0(i,m) XX do 35 j=jlow,k XX 35 sum = a(j,i)*dbiatx(j,m) + sum XX 40 dbiatx(i,m) = sum XX 99 return XX end XX double precision function bvalue ( t, bcoef, n, k, x, jderiv ) XX implicit double precision(a-h,o-z) calls interv c calculates value at x of jderiv-th derivative of spline from b-repr. c the spline is taken to be continuous from the right. c c****** i n p u t ****** c t, bcoef, n, k......forms the b-representation of the spline f to c be evaluated. specifically, c t.....knot sequence, of length n+k, assumed nondecreasing. c bcoef.....b-coefficient sequence, of length n . c n.....length of bcoef and dimension of s(k,t), c a s s u m e d positive . c k.....order of the spline . c c w a r n i n g . . . the restriction k .le. kmax (=20) is imposed c arbitrarily by the dimension statement for aj, dm, dm below, c but is n o w h e r e c h e c k e d for. c c x.....the point at which to evaluate . c jderiv.....integer giving the order of the derivative to be evaluated c a s s u m e d to be zero or positive. c c****** o u t p u t ****** c bvalue.....the value of the (jderiv)-th derivative of f at x . c c****** m e t h o d ****** c the nontrivial knot interval (t(i),t(i+1)) containing x is lo- c cated with the aid of interv . the k b-coeffs of f relevant for c this interval are then obtained from bcoef (or taken to be zero if c not explicitly available) and are then differenced jderiv times to c obtain the b-coeffs of (d**jderiv)f relevant for that interval. c precisely, with j = jderiv, we have from x.(12) of the text that c c (d**j)f = sum ( bcoef(.,j)*b(.,k-j,t) ) c c where c / bcoef(.), , j .eq. 0 c / c bcoef(.,j) = / bcoef(.,j-1) - bcoef(.-1,j-1) c / ----------------------------- , j .gt. 0 c / (t(.+k-j) - t(.))/(k-j) c c then, we use repeatedly the fact that c c sum ( a(.)*b(.,m,t)(x) ) = sum ( a(.,x)*b(.,m-1,t)(x) ) c with c (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1) c a(.,x) = --------------------------------------- c (x - t(.)) + (t(.+m-1) - x) c c to write (d**j)f(x) eventually as a linear combination of b-splines c of order 1 , and the coefficient for b(i,1,t)(x) must then c be the desired number (d**j)f(x). (see x.(17)-(19) of text). c XX parameter(kmax = 20) XX integer jderiv,k,n, i,ilo,imk,j,jc,jcmin,jcmax,jj,km1,mflag,nmi XX double precision bcoef(n),t(1),x XX double precision aj(kmax),dm(kmax),dp(kmax),fkmj c dimension t(n+k) current fortran standard makes it impossible to specify the length of c t precisely without the introduction of otherwise superfluous c additional arguments. XX bvalue = 0. XX if (jderiv .ge. k) go to 99 c c *** find i s.t. 1 .le. i .lt. n+k and t(i) .lt. t(i+1) and c t(i) .le. x .lt. t(i+1) . if no such i can be found, x lies c outside the support of the spline f and bvalue = 0. c (the asymmetry in this choice of i makes f rightcontinuous) XX if( (x.ne.t(n+1)) .or. (t(n+1).ne.t(n+k)) ) go to 700 XX i = n XX go to 701 XX 700 call interv ( t, n+k, x, i, mflag ) XX if (mflag .ne. 0) go to 99 XX 701 continue c *** if k = 1 (and jderiv = 0), bvalue = bcoef(i). XX km1 = k - 1 XX if (km1 .gt. 0) go to 1 XX bvalue = bcoef(i) XX go to 99 c c *** store the k b-spline coefficients relevant for the knot interval c (t(i),t(i+1)) in aj(1),...,aj(k) and compute dm(j) = x - t(i+1-j), c dp(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable c from input to zero. set any t.s not obtainable equal to t(1) or c to t(n+k) appropriately. XX 1 jcmin = 1 XX imk = i - k XX if (imk .ge. 0) go to 8 XX jcmin = 1 - imk XX do 5 j=1,i XX 5 dm(j) = x - t(i+1-j) XX do 6 j=i,km1 XX aj(k-j) = 0. XX 6 dm(j) = dm(i) XX go to 10 XX 8 do 9 j=1,km1 XX 9 dm(j) = x - t(i+1-j) c XX 10 jcmax = k XX nmi = n - i XX if (nmi .ge. 0) go to 18 XX jcmax = k + nmi XX do 15 j=1,jcmax XX 15 dp(j) = t(i+j) - x XX do 16 j=jcmax,km1 XX aj(j+1) = 0. XX 16 dp(j) = dp(jcmax) XX go to 20 XX 18 do 19 j=1,km1 XX 19 dp(j) = t(i+j) - x c XX 20 do 21 jc=jcmin,jcmax XX 21 aj(jc) = bcoef(imk + jc) c c *** difference the coefficients jderiv times. XX if (jderiv .eq. 0) go to 30 XX do 23 j=1,jderiv XX kmj = k-j XX fkmj = dble(kmj) XX ilo = kmj XX do 23 jj=1,kmj XX aj(jj) = ((aj(jj+1) - aj(jj))/(dm(ilo) + dp(jj)))*fkmj XX 23 ilo = ilo - 1 c c *** compute value at x in (t(i),t(i+1)) of jderiv-th derivative, c given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv). XX 30 if (jderiv .eq. km1) go to 39 XX jdrvp1 = jderiv + 1 XX do 33 j=jdrvp1,km1 XX kmj = k-j XX ilo = kmj XX do 33 jj=1,kmj XX aj(jj) = (aj(jj+1)*dm(ilo) + aj(jj)*dp(jj))/(dm(ilo)+dp(jj)) XX 33 ilo = ilo - 1 XX 39 bvalue = aj(1) c XX 99 return XX end XX subroutine bsplvb ( t, jhigh, index, x, left, biatx ) XX implicit double precision(a-h,o-z) calculates the value of all possibly nonzero b-splines at x of order c c jout = dmax( jhigh , (j+1)*(index-1) ) c c with knot sequence t . c c****** i n p u t ****** c t.....knot sequence, of length left + jout , assumed to be nonde- c creasing. a s s u m p t i o n . . . . c t(left) .lt. t(left + 1) . c d i v i s i o n b y z e r o will result if t(left) = t(left+1) c jhigh, c index.....integers which determine the order jout = max(jhigh, c (j+1)*(index-1)) of the b-splines whose values at x are to c be returned. index is used to avoid recalculations when seve- c ral columns of the triangular array of b-spline values are nee- c ded (e.g., in bvalue or in bsplvd ). precisely, c if index = 1 , c the calculation starts from scratch and the entire triangular c array of b-spline values of orders 1,2,...,jhigh is generated c order by order , i.e., column by column . c if index = 2 , c only the b-spline values of order j+1, j+2, ..., jout are ge- c nerated, the assumption being that biatx , j , deltal , deltar c are, on entry, as they were on exit at the previous call. c in particular, if jhigh = 0, then jout = j+1, i.e., just c the next column of b-spline values is generated. c c w a r n i n g . . . the restriction jout .le. jmax (= 20) is im- c posed arbitrarily by the dimension statement for deltal and c deltar below, but is n o w h e r e c h e c k e d for . c c x.....the point at which the b-splines are to be evaluated. c left.....an integer chosen (usually) so that c t(left) .le. x .le. t(left+1) . c c****** o u t p u t ****** c biatx.....array of length jout , with biatx(i) containing the val- c ue at x of the polynomial of order jout which agrees with c the b-spline b(left-jout+i,jout,t) on the interval (t(left), c t(left+1)) . c c****** m e t h o d ****** c the recurrence relation c c x - t(i) t(i+j+1) - x c b(i,j+1)(x) = -----------b(i,j)(x) + ---------------b(i+1,j)(x) c t(i+j)-t(i) t(i+j+1)-t(i+1) c c is used (repeatedly) to generate the (j+1)-vector b(left-j,j+1)(x), c ...,b(left,j+1)(x) from the j-vector b(left-j+1,j)(x),..., c b(left,j)(x), storing the new values in biatx over the old. the c facts that c b(i,1) = 1 if t(i) .le. x .lt. t(i+1) c and that c b(i,j)(x) = 0 unless t(i) .le. x .lt. t(i+j) c are used. the particular organization of the calculations follows al- c gorithm (8) in chapter x of the text. c XX parameter(jmax = 20) XX integer index,jhigh,left, i,j,jp1 XX double precision biatx(jhigh),t(1),x, deltal(jmax) XX double precision deltar(jmax),saved,term c dimension biatx(jout), t(left+jout) current fortran standard makes it impossible to specify the length of c t and of biatx precisely without the introduction of otherwise c superfluous additional arguments. XX data j/1/ c save j,deltal,deltar (valid in fortran 77) c XX go to (10,20), index XX 10 j = 1 XX biatx(1) = 1d0 XX if (j .ge. jhigh) go to 99 c XX 20 jp1 = j + 1 XX deltar(j) = t(left+j) - x XX deltal(j) = x - t(left+1-j) XX saved = 0d0 XX do 26 i=1,j XX term = biatx(i)/(deltar(i) + deltal(jp1-i)) XX biatx(i) = saved + deltar(i)*term XX 26 saved = deltal(jp1-i)*term XX biatx(jp1) = saved XX j = jp1 XX if (j .lt. jhigh) go to 20 c XX 99 return XX end XX subroutine dpbfa(abd,lda,n,m,info) XX integer lda,n,m,info XX double precision abd(lda,1) c c dpbfa factors a double precision symmetric positive definite c matrix stored in band form. c c dpbfa is usually called by dpbco, but it can be called c directly with a saving in time if rcond is not needed. c c on entry c c abd double precision(lda, n) c the matrix to be factored. the columns of the upper c triangle are stored in the columns of abd and the c diagonals of the upper triangle are stored in the c rows of abd . see the comments below for details. c c lda integer c the leading dimension of the array abd . c lda must be .ge. m + 1 . c c n integer c the order of the matrix a . c c m integer c the number of diagonals above the main diagonal. c 0 .le. m .lt. n . c c on return c c abd an upper triangular matrix r , stored in band c form, so that a = trans(r)*r . c c info integer c = 0 for normal return. c = k if the leading minor of order k is not c positive definite. c c band storage c c if a is a symmetric positive definite band matrix, c the following program segment will set up the input. c c m = (band width above diagonal) c do 20 j = 1, n c i1 = max0(1, j-m) c do 10 i = i1, j c k = i-j+m+1 c abd(k,j) = a(i,j) c 10 continue c 20 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas ddot c fortran max0,dsqrt c c internal variables c XX double precision ddot,t XX double precision s XX integer ik,j,jk,k,mu c begin block with ...exits to 40 c c XX do 30 j = 1, n XX info = j XX s = 0.0d0 XX ik = m + 1 XX jk = max0(j-m,1) XX mu = max0(m+2-j,1) XX if (m .lt. mu) go to 20 XX do 10 k = mu, m XX t = abd(k,j) - ddot(k-mu,abd(ik,jk),1,abd(mu,j),1) XX t = t/abd(m+1,jk) XX abd(k,j) = t XX s = s + t*t XX ik = ik - 1 XX jk = jk + 1 XX 10 continue XX 20 continue XX s = abd(m+1,j) - s c ......exit XX if (s .le. 0.0d0) go to 40 XX abd(m+1,j) = dsqrt(s) XX 30 continue XX info = 0 XX 40 continue XX return XX end XX subroutine dpbsl(abd,lda,n,m,b) XX integer lda,n,m XX double precision abd(lda,1),b(1) c c dpbsl solves the double precision symmetric positive definite c band system a*x = b c using the factors computed by dpbco or dpbfa. c c on entry c c abd double precision(lda, n) c the output from dpbco or dpbfa. c c lda integer c the leading dimension of the array abd . c c n integer c the order of the matrix a . c c m integer c the number of diagonals above the main diagonal. c c b double precision(n) c the right hand side vector. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains c a zero on the diagonal. technically this indicates c singularity but it is usually caused by improper subroutine c arguments. it will not occur if the subroutines are called c correctly and info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call dpbco(abd,lda,n,rcond,z,info) c if (rcond is too small .or. info .ne. 0) go to ... c do 10 j = 1, p c call dpbsl(abd,lda,n,c(1,j)) c 10 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas daxpy,ddot c fortran min0 c c internal variables c XX double precision ddot,t XX integer k,kb,la,lb,lm c c solve trans(r)*y = b c XX do 10 k = 1, n XX lm = min0(k-1,m) XX la = m + 1 - lm XX lb = k - lm XX t = ddot(lm,abd(la,k),1,b(lb),1) XX b(k) = (b(k) - t)/abd(m+1,k) XX 10 continue c c solve r*x = y c XX do 20 kb = 1, n XX k = n + 1 - kb XX lm = min0(k-1,m) XX la = m + 1 - lm XX lb = k - lm XX b(k) = b(k)/abd(m+1,k) XX t = -b(k) XX call daxpy(lm,t,abd(la,k),1,b(lb),1) XX 20 continue XX return XX end EEEEE X $shar_touch -am 0405184995 'sortdi.f' && X chmod 0644 'sortdi.f' || X $echo 'restore of' 'sortdi.f' 'failed' X if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ X && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then X md5sum -c << SHAR_EOF >/dev/null 2>&1 \ X || $echo 'sortdi.f:' 'MD5 check failed' f5694b73134c4f3789e64054344f1e34 sortdi.f XSHAR_EOF X else X shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'sortdi.f'`" X test 24539 -eq "$shar_count" || X $echo 'sortdi.f:' 'original size' '24539,' 'current size' "$shar_count!" X fi fi rm -fr _sh20785 Xexit 0 SHAR_EOF $shar_touch -am 0425070698 'ratfor.shar' && chmod 0644 'ratfor.shar' || $echo 'restore of' 'ratfor.shar' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'ratfor.shar:' 'MD5 check failed' 36e5701eb312bb117549475a278af6a7 ratfor.shar SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'ratfor.shar'`" test 81669 -eq "$shar_count" || $echo 'ratfor.shar:' 'original size' '81669,' 'current size' "$shar_count!" fi fi # ============= makefile ============== if test -f 'makefile' && test "$first_param" != -c; then $echo 'x -' SKIPPING 'makefile' '(file already exists)' else $echo 'x -' extracting 'makefile' '(text)' sed 's/^X//' << 'SHAR_EOF' > 'makefile' && X FFLAGS = -G0 all: mda.so mda.so: X Splus SHLIB -o mda.so bruto.r dcalcvar.r dmarss.r dorthreg.r dqrreg.r mspline.r sortdi.f MARS.o: dcalcvar.o dmarss.o dorthreg.o dqrreg.o X ld -r -d -o MARS.o dcalcvar.o dmarss.o dorthreg.o dqrreg.o #sometimes the compiler complains about the -G0 flag, # and suggests that it be used with ld instead # On SGI machines we need to use the -shared option to make dyn.load work, along # with a "fixed" version of dyn.load2 from utstat.toronto.edu BRUTO.o: mspline.o bruto.o sortdi.o X ld -r -d -o BRUTO.o mspline.o bruto.o sortdi.o the.S.dump: X echo "make.dumpdata.mda()" | Splus the.helpfile.shar: X shar -d EEEEE .Data/.Help/* > helpfile.shar the.ratfor.shar: bruto.r dcalcvar.r dmarss.r dorthreg.r dqrreg.r mspline.r sortdi.f X shar -d EEEEE bruto.r dcalcvar.r dmarss.r dorthreg.r dqrreg.r mspline.r sortdi.f> ratfor.shar mda.shar: the.S.dump the.helpfile.shar the.ratfor.shar makefile README X shar README dumpdata.mda helpfile.shar ratfor.shar makefile > mda.shar SHAR_EOF $shar_touch -am 0423103198 'makefile' && chmod 0644 'makefile' || $echo 'restore of' 'makefile' 'failed' if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \ && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then md5sum -c << SHAR_EOF >/dev/null 2>&1 \ || $echo 'makefile:' 'MD5 check failed' 18e75d1330f402f3c0713a4762a07525 makefile SHAR_EOF else shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'makefile'`" test 1000 -eq "$shar_count" || $echo 'makefile:' 'original size' '1000,' 'current size' "$shar_count!" fi fi rm -fr _sh20801 exit 0