#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh 'charmatch.c' <<'END_OF_FILE'
X/* SCCS: @(#)charmatch.c 1.1 1/26/90 */
X/*
X** The working part of the charmatch function in S
X*/
X#include
X
charmatch(input, n_input, target, n_target, result)
char *input[], *target[];
long *n_input, *n_target;
long result[];
X{
X register int i,j, k;
X register int match;
X int temp, perfect;
X
X for (i=0; i<*n_input; i++) {
X temp = strlen(input[i]);
X match= -1; perfect=0;
X for (j=0; j< *n_target; j++) {
X k = strncmp(input[i], target[j], temp);
X if (k==0) {
X if (strlen(target[j]) == temp) {
X if (perfect==1) match=0;
X else {
X perfect=1;
X match = j+1;
X }
X }
X else if (perfect==0) {
X if (match== -1) match=j+1;
X else match=0;
X }
X }
X }
X result[i] = match;
X }
X }
END_OF_FILE
if test 752 -ne `wc -c <'charmatch.c'`; then
echo shar: \"'charmatch.c'\" unpacked with wrong size!
fi
# end of 'charmatch.c'
fi
if test -f 'charmatch.d' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'charmatch.d'\"
else
echo shar: Extracting \"'charmatch.d'\" \(753 characters\)
sed "s/^X//" >'charmatch.d' <<'END_OF_FILE'
X.BG
X.FN charmatch
X.TL
Match based on Shortest Unique Identifier
X.CS
charmatch(inputs, target)
X.AG inputs
vector of strings to be matched against the target.
X.AG target
set of keywords over which to search.
X.RT
a vector of the same length as `inputs' giving, for each element of
X`inputs', the position in `target' of the matching keyword. If the
keyword was not found, then NA is returned. If there was an ambiguous
match, then 0 is returned as the position.
X.PP
This function is useful within other functions for matching an input
argument to a series of
alphabetic choices. If target were c("logit", "loglog", "probit"), for
instance, then an input of "p" would return a 3, one of "log" would return
a 0, and one of "arctan" would return an NA.
X.WR
END_OF_FILE
if test 753 -ne `wc -c <'charmatch.d'`; then
echo shar: \"'charmatch.d'\" unpacked with wrong size!
fi
# end of 'charmatch.d'
fi
if test -f 'charmatch.s' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'charmatch.s'\"
else
echo shar: Extracting \"'charmatch.s'\" \(773 characters\)
sed "s/^X//" >'charmatch.s' <<'END_OF_FILE'
X# SCCS: @(#)charmatch.s 1.1 1/26/90
X#
X# A "shortest unique identifier" match. If an input string matches
X# no output string, return NA. If it is an ambiguous match, return 0.
X# Otherwise return the index of the match.
X# Exception: if the target string contains "log" and "logistic", and the
X# user type "lo" it is ambiguous, but if he types "log" consider it a
X# perfect match.
X#
charmatch <- function( inputs, target) {
X if (!(is.character(inputs) && is.character(target)))
X stop ("Input must be character strings")
X
X if (len(inputs) <1) return(NULL)
X if (len(target) <1) return(rep(NA, len(inputs)))
X
X temp <- .C("charmatch", inputs, len(inputs), target, len(target),
X result=integer(len(inputs)))
X
X ifelse(temp$result<0, NA, temp$result)
X }
END_OF_FILE
if test 773 -ne `wc -c <'charmatch.s'`; then
echo shar: \"'charmatch.s'\" unpacked with wrong size!
fi
# end of 'charmatch.s'
fi
if test -f 'crossval.d' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'crossval.d'\"
else
echo shar: Extracting \"'crossval.d'\" \(1506 characters\)
sed "s/^X//" >'crossval.d' <<'END_OF_FILE'
X.BG
X.FN crossval
X.TL
crossval: Performs cross validation for windowing
X.CS
crossval(x,y,method="fixwidth",smparam=defpar,
X wtmethod="boxcar",estmethod="mean")
X
X.AG " x"
Independent variable
X.AG " y"
Dependent variable
X.AG " method"
X(optional) windowing method
X.br
Method by which decision is made to include or exclude data points
from the window.
X.br
Allowable values: "fixwidth","nearest neighbor"
X
Only enough of the argument to insure uniquness need be
given. Thus crossval(x,y,"n") works with the windowing method set to
nearest neighbor. This comment applies to other character arguments.
X.AG " smparam"
X(optional) vector of smoothing parameters
X.br If the windowing method is "fixwidth", smparam is the size of the
window. If ommitted, covers the range of x in increments of 10%.
If the windowing method is "nearest neighbor", it is the number of
neighbors to be included in the window. If ommitted, covers the number
of x values in increments of 10%.
X.AG " wtmethod"
X(optional) weighting method
X.br
Method by which weights are to be assigned to the values in the window.
Allowable values: "biquadratic", "boxcar".
X.AG "estmethod"
X(optional) estimation method
X.br
Method by which smoothed values of y elements in window are to be computed
If ommitted, set to "mean"
Allowable values: "linear fit", "mean", "quadratic fit".
X.RT
A vector with the cross validation term for each value in smparam
X.EX
X> crossval(x,y,"fixw",10,"box","linear")
X.KW cross validation
X.WR
X
END_OF_FILE
if test 1506 -ne `wc -c <'crossval.d'`; then
echo shar: \"'crossval.d'\" unpacked with wrong size!
fi
# end of 'crossval.d'
fi
if test -f 'crossval.s' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'crossval.s'\"
else
echo shar: Extracting \"'crossval.s'\" \(6494 characters\)
sed "s/^X//" >'crossval.s' <<'END_OF_FILE'
crossval <- function(x, y, method = "fixwidth", smparam = defpar,
X wtmethod = "boxcar",estmethod = "mean")
X{
X#########################################################################
X#
X# CROSSVAL
X#
X# FUNCTION
X#
X# Calls fortran windowing routines to perform cross validation
X#
X# ARGUMENTS
X#
X# x --> Independant variable for smoothing
X#
X# y --> Dependant variable for smoothing
X#
X# method --> Character string specifing the windowing method
X#
X# smparam --> Windowing parameter. If method is "fixwidth" then the
X# default covers the range of x in increments of 10%.
X# If method is "nearest neighbor" then the default covers
X# the length of x in increments of 10%.
X#
X# wtmethod --> Character string specifing the weighting method
X#
X# estmethod --> Character string specifing the estimation method
X#
X#########################################################################
X#
X# make sure x and y are there, are numeric, and have the same length
X#
X if(missing(x)) stop("Argument x is missing, with no default")
X if(missing(y)) stop("Argument y is missing, with no default")
X if(!is.numeric(x)) stop("x must be numeric")
X if(!is.numeric(y)) stop("y must be numeric")
X if(len(x) != len(y)) stop("Lengths of x and y must match")
X#
X# remove any NAs
X#
X xx <- x[!(is.na(x)) & !(is.na(y))]
X yy <- y[!(is.na(x)) & !(is.na(y))]
X nxy <- len(xx)
X if(nxy < 2) stop("x and y yield less than 2 points")
X#
X# sort them if we need to
X#
X if (any(c(-1.0E-35,xx) > c(xx,1.0E35))) {
X permut <- sort.list(xx)
X xx <- xx[permut]
X yy <- yy[permut]
X }
X#
X# get the various methods
X#
X if(!missing(method)) {
X if(!is.character(method))
X stop("method must be of mode character")
X target <- c("fixwidth", "nearest neighbor")
X imethod <- charmatch(method, target)
X if(is.na(imethod)) stop("windowing method unknown")
X else if(!(imethod > 0)) stop("windowing method ambiguous")
X }
X else imethod <- 1
X#
X if(!missing(wtmethod)) {
X if(!is.character(wtmethod))
X stop("wtmethod must be of mode character")
X target <- c("boxcar", "biquadratic")
X iwtmethod <- charmatch(wtmethod, target)
X if(is.na(iwtmethod)) stop("weighting method unknown")
X else if(!(iwtmethod > 0)) stop("weighting method ambiguous")
X }
X else iwtmethod <- 1
X#
X if(!missing(estmethod)) {
X if(!is.character(estmethod))
X stop("estmethod must be of mode character")
X target <- c("mean","linear fit", "quadratic fit")
X iestmethod <- charmatch(estmethod, target)
X if(is.na(iestmethod))
X stop("crossvalidation estimation method unknown")
X else if(!(iestmethod > 0)) stop("estimation method ambiguous")
X }
X else iestmethod <- 1
X#
X# If windowing method is fixwidth find the smallest window necessary
X#
X if (iestmethod == 1) npts <- 1
X else if(iestmethod == 2) npts <- 2
X else if(iestmethod == 3) npts <- 3
X if(imethod == 1) {
X smwind <- .Fortran("wdthmx",
X as.single(xx),
X as.integer(nxy),
X as.integer(npts),
X smlwnd = single(length = 1))
X smlwnd <- smwind$smlwnd
X }
X#
X# initialize the variables that need it
X#
X# if smparam wasn't given, define the default value for it
X#
X if(missing(smparam)) {
X if(imethod == 1){
X range_div10 <- (max(xx) - min(xx))/10
X from <- max(smlwnd,range_div10)
X defpar <- seq(from,to=max(xx),by=range_div10)
X }
X else {
X len_div10 <- nxy/10
X from <- max(npts,len_div10)
X defpar <- seq(from,to=nxy,by=len_div10)
X }
X }
X#
X# else since smparam was given make sure its elements are in
X# the right range
X#
X else {
X if(!is.numeric(smparam)) stop("smparam must be numeric")
X if(imethod == 1) {
X smparam[smparam < smlwnd] <- NA
X smparam[smparam > (max(xx)-min(xx))] <- NA
X }
X else {
X smparam[smparam < npts] <- NA
X smparam[smparam > nxy] <- NA
X }
X }
X#
X xdes <- xx
X ndes <- nxy
X sortar <- matrix(2, nxy)
X storage.mode(sortar) <- "single"
X xwgts <- vector(mode="numeric",length=nxy)
X xmiss <- vector(mode="numeric",length=nxy)
X qcross <- TRUE
X crsvalterm <- vector(mode = "numeric", length = len(smparam))
X#
X# loop over all the smoothing parameters
X#
X for(i in seq(along = smparam)) {
X#
X# if smparam isn't NA call the fortran to do the windowing
X#
X if(!(is.na(smparam[i]))) {
X wind <- .Fortran("window",
X as.single(xx),
X as.single(yy),
X as.integer(nxy),
X as.single(xdes),
X ysmo = single(length = ndes),
X as.integer(ndes),
X as.integer(imethod),
X as.single(smparam[i]),
X integer(length = 1),
X as.integer(iwtmethod),
X as.integer(iestmethod),
X single(length = 1),
X sortar,
X xwgts,
X as.integer(xmiss),
X as.logical(qcross),
X istat = integer(length = 1))
X if(wind$istat == 0)
X crsvalterm[i] <- sum((yy - wind$ysmo)^2)
X else crsvalterm[i] <- NA
X }
X else crsvalterm[i] <- NA
X }
X return(cbind(smparam,crsvalterm))
X}
END_OF_FILE
if test 6494 -ne `wc -c <'crossval.s'`; then
echo shar: \"'crossval.s'\" unpacked with wrong size!
fi
# end of 'crossval.s'
fi
if test -f 'windowing.d' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'windowing.d'\"
else
echo shar: Extracting \"'windowing.d'\" \(2214 characters\)
sed "s/^X//" >'windowing.d' <<'END_OF_FILE'
X.BG
X.FN windowing
X.TL
windowing: Window/Nonparametric Regression Smoothing/Estimation
X.CS
windowing(x,y,xdes=defxdes,method="fixwidth",smparam=defpar,
X wtmethod="boxcar",estmethod="mean",dermethod="linear",
X quantile=0.5)
X
X.AG " x"
Independent variable.
X.AG " y"
Dependent variable.
X.AG " xdes"
X(optional)Values at which smoothed y values are desired, i.e, window centers.
If omitted, set to 20 equally spaced values over the range of x.
X.AG " method"
X(optional) windowing method
X.br
Method by which decision is made to include or exclude data points
from the window.
Allowable values: "fixwidth", "nearest neighbor"
X
Note that only enough of the argument to insure uniquness need be
given. Thus windowing(x,y,"n") works with the windowing method set to
nearest neighbor. This comment applies to other character arguments.
X.AG " smparam"
X(optional) smoothing parameter
X.br
If the windowing method is "fixwidth" smparam is the size of the
window. If omitted, set to 1/10 the range of values.
If the windowing method is "nearest neighbor" it is the number of
neighbors to be include in the window. If omitted, set to 1/10 the
number of points.
X.AG " wtmethod"
X(optional) weighting method
X.br
Method by which weights are to be assigned to the values in the window.
Allowable values: "biquadratic", "boxcar".
X.AG "estmethod"
X(optional) estimation method
X.br
Method by which smoothed values of y elements in window are to be computed
X.br
Allowable values: "maximum","mean","minimum","number of points",
X"linear fit","quadratic fit","quantile","standard deviation","window width"
X.AG "dermethod"
X(optional) derivative method
X.br
Method by which the derivative of y as a function of x is estimated.
X.br
Allowable values:"linear fit", and "quadratic fit"
X.AG " quantile"
X(optional) Quantile (number in range 0-1) of y to be estimated
X.br
If omitted, set to 0.5
X.RT
A matrix with xdes in column one and the smoothed data in column two
X.EX
X> smodata <- windowing(x,y,estmethod="linear")
X.br
X # smooths using linear fit estimation
X.br
X> plot(x,y,type="l")
X.br
X> lines(smodata[,1],smodata[,2],lty=2)
X.br
X # plots the original data and the smoothed/estimated data
X.KW smooth
X.WR
END_OF_FILE
if test 2214 -ne `wc -c <'windowing.d'`; then
echo shar: \"'windowing.d'\" unpacked with wrong size!
fi
# end of 'windowing.d'
fi
if test -f 'windowing.f' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'windowing.f'\"
else
echo shar: Extracting \"'windowing.f'\" \(95080 characters\)
sed "s/^X//" >'windowing.f' <<'END_OF_FILE'
X*DECK DGEDI
X SUBROUTINE DGEDI(A,LDA,N,IPVT,DET,WORK,JOB)
X INTEGER LDA,N,IPVT(1),JOB
X DOUBLE PRECISION A(LDA,1),DET(2),WORK(1)
C
C DGEDI COMPUTES THE DETERMINANT AND INVERSE OF A MATRIX
C USING THE FACTORS COMPUTED BY DGECO OR DGEFA.
C
C ON ENTRY
C
C A DOUBLE PRECISION(LDA, N)
C THE OUTPUT FROM DGECO OR DGEFA.
C
C LDA INTEGER
C THE LEADING DIMENSION OF THE ARRAY A .
C
C N INTEGER
C THE ORDER OF THE MATRIX A .
C
C IPVT INTEGER(N)
C THE PIVOT VECTOR FROM DGECO OR DGEFA.
C
C WORK DOUBLE PRECISION(N)
C WORK VECTOR. CONTENTS DESTROYED.
C
C JOB INTEGER
C = 11 BOTH DETERMINANT AND INVERSE.
C = 01 INVERSE ONLY.
C = 10 DETERMINANT ONLY.
C
C ON RETURN
C
C A INVERSE OF ORIGINAL MATRIX IF REQUESTED.
C OTHERWISE UNCHANGED.
C
C DET DOUBLE PRECISION(2)
C DETERMINANT OF ORIGINAL MATRIX IF REQUESTED.
C OTHERWISE NOT REFERENCED.
C DETERMINANT = DET(1) * 10.0**DET(2)
C WITH 1.0 .LE. DABS(DET(1)) .LT. 10.0
C OR DET(1) .EQ. 0.0 .
C
C ERROR CONDITION
C
C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS
C A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED.
C IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY
C AND IF DGECO HAS SET RCOND .GT. 0.0 OR DGEFA HAS SET
C INFO .EQ. 0 .
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,DSCAL,DSWAP
C FORTRAN DABS,MOD
C
C INTERNAL VARIABLES
C
X DOUBLE PRECISION T
X DOUBLE PRECISION TEN
X INTEGER I,J,K,KB,KP1,L,NM1
C
C
C COMPUTE DETERMINANT
C
X IF (JOB/10 .EQ. 0) GO TO 70
X DET(1) = 1.0D0
X DET(2) = 0.0D0
X TEN = 10.0D0
X DO 50 I = 1, N
X IF (IPVT(I) .NE. I) DET(1) = -DET(1)
X DET(1) = A(I,I)*DET(1)
C ...EXIT
X IF (DET(1) .EQ. 0.0D0) GO TO 60
X 10 IF (DABS(DET(1)) .GE. 1.0D0) GO TO 20
X DET(1) = TEN*DET(1)
X DET(2) = DET(2) - 1.0D0
X GO TO 10
X 20 CONTINUE
X 30 IF (DABS(DET(1)) .LT. TEN) GO TO 40
X DET(1) = DET(1)/TEN
X DET(2) = DET(2) + 1.0D0
X GO TO 30
X 40 CONTINUE
X 50 CONTINUE
X 60 CONTINUE
X 70 CONTINUE
C
C COMPUTE INVERSE(U)
C
X IF (MOD(JOB,10) .EQ. 0) GO TO 150
X DO 100 K = 1, N
X A(K,K) = 1.0D0/A(K,K)
X T = -A(K,K)
X CALL DSCAL(K-1,T,A(1,K),1)
X KP1 = K + 1
X IF (N .LT. KP1) GO TO 90
X DO 80 J = KP1, N
X T = A(K,J)
X A(K,J) = 0.0D0
X CALL DAXPY(K,T,A(1,K),1,A(1,J),1)
X 80 CONTINUE
X 90 CONTINUE
X 100 CONTINUE
C
C FORM INVERSE(U)*INVERSE(L)
C
X NM1 = N - 1
X IF (NM1 .LT. 1) GO TO 140
X DO 130 KB = 1, NM1
X K = N - KB
X KP1 = K + 1
X DO 110 I = KP1, N
X WORK(I) = A(I,K)
X A(I,K) = 0.0D0
X 110 CONTINUE
X DO 120 J = KP1, N
X T = WORK(J)
X CALL DAXPY(N,T,A(1,J),1,A(1,K),1)
X 120 CONTINUE
X L = IPVT(K)
X IF (L .NE. K) CALL DSWAP(N,A(1,K),1,A(1,L),1)
X 130 CONTINUE
X 140 CONTINUE
X 150 CONTINUE
X RETURN
X END
X*DECK DGEFA
X SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO)
X INTEGER LDA,N,IPVT(1),INFO
X DOUBLE PRECISION A(LDA,1)
C
C DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION.
C
C DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED
C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) .
C
C ON ENTRY
C
C A DOUBLE PRECISION(LDA, N)
C THE MATRIX TO BE FACTORED.
C
C LDA INTEGER
C THE LEADING DIMENSION OF THE ARRAY A .
C
C N INTEGER
C THE ORDER OF THE MATRIX A .
C
C ON RETURN
C
C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
C WHICH WERE USED TO OBTAIN IT.
C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE
C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
C
C IPVT INTEGER(N)
C AN INTEGER VECTOR OF PIVOT INDICES.
C
C INFO INTEGER
C = 0 NORMAL VALUE.
C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR
C CONDITION FOR THIS SUBROUTINE, BUT IT DOES
C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO
C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE
C INDICATION OF SINGULARITY.
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,DSCAL,IDAMAX
C
C INTERNAL VARIABLES
C
X DOUBLE PRECISION T
X INTEGER IDAMAX,J,K,KP1,L,NM1
C
C
C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
X INFO = 0
X NM1 = N - 1
X IF (NM1 .LT. 1) GO TO 70
X DO 60 K = 1, NM1
X KP1 = K + 1
C
C FIND L = PIVOT INDEX
C
X L = IDAMAX(N-K+1,A(K,K),1) + K - 1
X IPVT(K) = L
C
C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
X IF (A(L,K) .EQ. 0.0D0) GO TO 40
C
C INTERCHANGE IF NECESSARY
C
X IF (L .EQ. K) GO TO 10
X T = A(L,K)
X A(L,K) = A(K,K)
X A(K,K) = T
X 10 CONTINUE
C
C COMPUTE MULTIPLIERS
C
X T = -1.0D0/A(K,K)
X CALL DSCAL(N-K,T,A(K+1,K),1)
C
C ROW ELIMINATION WITH COLUMN INDEXING
C
X DO 30 J = KP1, N
X T = A(L,J)
X IF (L .EQ. K) GO TO 20
X A(L,J) = A(K,J)
X A(K,J) = T
X 20 CONTINUE
X CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
X 30 CONTINUE
X GO TO 50
X 40 CONTINUE
X INFO = K
X 50 CONTINUE
X 60 CONTINUE
X 70 CONTINUE
X IPVT(N) = N
X IF (A(N,N) .EQ. 0.0D0) INFO = N
X RETURN
X END
X SUBROUTINE drvtv(x,y,weight,xval,nptwi,method,deriv,info)
C***********************************************************************
C
C SUBROUTINE DRVTV(X, Y, WEIGHT, XVAL, NPTWI, METHOD, DERIV, INFO)
C
C DeRiVaTiVe
C
C FUNCTION
C
C Estimates the derivative of the dependent variable with
C respect to the independent variable.
C
C ARGUMENTS
C
C X --> Vector of X values in window.
C REAL X(NPTWI)
C
C Y --> Vector of Y values in window.
C REAL Y(NPTWI)
C
C WEIGHT --> Weights associated with X and Y values.
C REAL WEIGHT(NPTWI)
C
C XVAL --> Value at window center at which derivative is estimate
C REAL XVAL
C
C NPTWI --> Number of points in the window.
C INTEGER NPTWI
C
C METHOD --> Method to estimate the derivative
C 1 - Linear
C 2 - Quadratic
C INTEGER METHOD
C
C DERIV --> The estimate of the derivative
C REAL DERIV
C
C INFO <-- Returned from LQBETA about the sucsess of finding
C the coeficients of regression.
C If BETAS are returned from LQBETA INFO = 0
C Else INFO > 0
C INTEGER INFO
C
C
C VARIABLES
C
C BETAS -- Coeficients of quadratic regresion
C REAL BETAS(NPARM)
C
C NPARAM -- Number of parameters of the regresion function.
C INTEGER NPARM
C
C
C METHOD
C
C***********************************************************************
X IMPLICIT LOGICAL (q)
X REAL x(nptwi),y(nptwi),weight(nptwi),xval,deriv
X INTEGER nptwi,method,info
X REAL betas(3)
X INTEGER nparam,ixval
C ..
C .. Main Body of Code
X
X IF (.NOT. (method.EQ.1)) GO TO 10
X nparam = 2
X GO TO 20
X
X 10 nparam = 3
C
C Call LQBETA to get the coeficients of linear regresion
C
X 20 ixval = -1
X CALL lqbeta(x,y,weight,nptwi,betas,nparam,ixval,info)
C
C If ok then use these BETAS to find the estimate for the derivative
C
X IF (.NOT. (info.EQ.0)) GO TO 50
X IF (.NOT. (method.EQ.1)) GO TO 30
X deriv = betas(2)
X GO TO 40
X
X 30 deriv = 2*betas(3)*xval + betas(2)
X 40 CONTINUE
X 50 RETURN
X
X END
X SUBROUTINE eslin(x,y,weight,nptwi,xval,yval,ixval,info)
C***********************************************************************
C SUBROUTINE ESLIN(X, Y, WEIGHT, NPTWI, XVAL, YVAL,IXVAL)
C
C EStimate by LINear fit
C
C FUNCTION
C
C Performs a weighted linear fit to X,Y data, then evaluates the
C linear function at XVAL.
C
C ARGUMENTS
C
C X --> Vector of X values in window.
C REAL X(NPTWI)
C
C Y --> Vector of Y values in window.
C REAL Y(NPTWI)
C
C WEIGHT --> Weights associated with X and Y values.
C REAL WEIGHT(NPTWI)
C
C XVAL --> Value at window center at which linear is evaluated.
C REAL XVAL
C
C NPTWI --> Number of points in the window.
C INTEGER NPTWI
C
C YVAL <-- Smoothed Y value for Y values in window.
C Returned only if INFO equals 0
C REAL YVAL
C
C IXVAL --> X(IXVAL) is not used in the fit when cross validating.
C If not cross validating IXVAL is negative.
C INTEGER IXVAL
C
C INFO <-- Returned from LQBETA about the sucsess of finding
C the coeficients of regression.
C If BETAS are returned from LQBETA INFO = 0
C Else INFO > 0
C INTEGER INFO
C
C VARIABLES
C
C BETAS -- The coeficients of regresion
C REAL BETAS(NPARAM)
C
C NPARAM -- The number of parameters of the regresion function.
C For linear regresion NPARM=2
C INTEGER NPARM
C
C
C METHOD
C
C Calls LQBETA to find the coeficients of linear regresion, then use
C these coeficients to find the estimate of Y at XVAL.
C
C***********************************************************************
C
C DECLARATIONS & SPECIFICATIONS
C
X REAL x(nptwi),y(nptwi),weight(nptwi),xval,yval
X INTEGER nptwi,ixval,nparam,info
X PARAMETER (nparam=2)
X REAL betas(nparam)
C
C Call LQBETA to get the coeficients of linear regresion
C
X CALL lqbeta(x,y,weight,nptwi,betas,nparam,ixval,info)
C
C If ok then use these BETAS to find the estimate for Y at XVAL
C
X IF (info.EQ.0) yval = betas(1) + betas(2)*xval
X RETURN
X
X END
X SUBROUTINE esmean(y,npt,weight,yval,iyval)
C***********************************************************************
C SUBROUTINE ESMEAN(Y, NPT, WEIGHT, YVAL, IYVAL)
C
C Weighted MEAN of array y
C
C FUNCTION
C
C This module returns the weighted mean of real vector Y.
C
C ARGUMENTS
C
C --> Y - REAL VECTOR - Vector for which weighted mean to be found.
C Y(NPT).
C --> NPT - INTEGER - Number of elements in Y.
C --> WEIGHT - REAL VECTOR - Weights for Y vector elements.
C WEIGHT(NPT).
C <-- YVAL - REAL - Weighted mean of Y values in window
C --> IYVAL - INTEGER - Y(IYVAL) is not used when Cross Validating
C
C METHODS
C
C The weighted mean is the sum of the products of each datum with
C its weight. The Y vector need not be sorted.
C*****************************************************************
C DECLARATIONS & SPECIFICATIONS
C
X IMPLICIT REAL (a-h),INTEGER (i-n),REAL (o-p),LOGICAL (q),
X + REAL (r-z)
C
X DIMENSION y(npt),weight(npt)
C---------------------------------------------------------------------
C
C SUM THE PRODUCTS OF WT*VALUE TO GET THE MEAN
X yval = 0.0
X DO 30,ii = 1,npt
X IF (.NOT. (ii.EQ.iyval)) GO TO 10
X yval = yval
X GO TO 20
X
X 10 yval = yval + y(ii)*weight(ii)
X 20 CONTINUE
X 30 CONTINUE
X RETURN
X
X END
X SUBROUTINE esquad(x,y,weight,nptwi,xval,yval,ixval,info)
C***********************************************************************
C SUBROUTINE ESQUAD(X, Y, WEIGHT, NPTWI, XVAL, YVAL, IXVAL)
C
C EStimate by QUADratic fit
C
C FUNCTION
C
C Performs a weighted quadratic fit to X,Y data, then evaluates
C the quadratic at XVAL.
C
C ARGUMENTS
C
C X --> Vector of X values in window.
C REAL X(NPTWI)
C
C Y --> Vector of Y values in window.
C REAL Y(NPTWI)
C
C WEIGHT --> Weights associated with X and Y values.
C REAL WEIGHT(NPTWI)
C
C XVAL --> Value at window center at which linear is evaluated.
C REAL XVAL
C
C NPTWI --> Number of points in the window.
C INTEGER NPTWI
C
C YVAL --> Smoothed Y value for Y values in window.
C REAL YVAL
C
C IXVAL --> X(IXVAL) is not used in the fit when cross validating.
C If not cross validating IXVAL is negative.
C INTEGER IXVAL
C
C INFO <-- Returned from LQBETA about the sucsess of finding
C the coeficients of regression.
C If BETAS are returned from LQBETA INFO = 0
C Else INFO > 0
C INTEGER INFO
C
C
C VARIABLES
C
C BETAS -- Coeficients of quadratic regresion
C REAL BETAS(NPARM)
C
C NPARAM -- Number of parameters of the regresion function.
C For quadratic regresion NPARM=3
C INTEGER NPARM
C
C
C METHOD
C
C Calls LQBETA to find the coeficients of quadratic regresion, then
C uses these coeficients to find the estimate of Y at XVAL.
C
C
C***********************************************************************
C
C DECLARATIONS & SPECIFICATIONS
C
X REAL x(nptwi),y(nptwi),weight(nptwi),xval,yval
X INTEGER nptwi,ixval,nparam,info
X PARAMETER (nparam=3)
X REAL betas(nparam)
C
C Call LQBETA to get the coeficients of quadratic regresion
C
X CALL lqbeta(x,y,weight,nptwi,betas,nparam,ixval,info)
C
C If ok then use these BETAS to find the estimate for Y at XVAL
C
X IF (.NOT. (info.EQ.0)) GO TO 10
X yval = betas(1) + betas(2)*xval + betas(3)*xval*xval
X 10 RETURN
X
X END
X SUBROUTINE esquan(y,weight,nptwi,quant,ysort,yval)
C***********************************************************************
C SUBROUTINE ESQUAN(Y, WEIGHT, NPTWI, QUANT, YSORT, YVAL)
C
C
C EStimate QUANtile
C
C
C Function
C
C
C Returns a value greater than QUANT percent of sorted Y values,
C i.e., the point at which the QUANT quantile occurs.
C
C
C**********************************************************************
C
C
C Method
C
C The Y vector and the X weights are sorted on the Y values.
C
C The WEIGHT vector is summed until the partial sum is greater
C than QUANT. The sorted Y values corresponding to the last and
C penultimate weights included in the partial sum are selected.
C The value is linearly interpolated between these YFIT values.
C
C***********************************************************************
C
C
C Arguments
C
C
C Y --> Vector of Y values in window.
C REAL Y(NPTWI)
C
C WEIGHT --> Weights of elements of Y vector.
C REAL WEIGHT(NPTWI)
C
C NPTWI --> Number of points in Y vector.
C INTEGER NPTWI
C
C QUANT --> Quantile for which value is to be found.
C REAL QUANT
C
C YSORT --> Array for sorting Y values and X weights.
C REAL YSORT(2,NPTWI)
C
C YVAL <-- Smoothed Y value for Y values in window.
C REAL YVAL
C
C
C***********************************************************************
C
C
X IMPLICIT LOGICAL (q)
C
X REAL y(nptwi),weight(nptwi),quant,ysort(2,nptwi),yval
X INTEGER nptwi
C
C Function used to sort Y values and X weights
C
X LOGICAL qgr
X EXTERNAL qgr
C
C
C*****************************************************************
C
C Main body of program
C
C
C Put Y values and X weights into sort array
C
X i = 1
X GO TO 20
X
X 10 i = i + 1
X 20 IF (i.GT.nptwi) GO TO 30
X ysort(1,i) = y(i)
X ysort(2,i) = weight(i)
X GO TO 10
X
X 30 CALL srmrl(qgr,ysort,2,2,nptwi,1)
C
C Sort on Y values
C
C
C Put sorted X weights back into WEIGHT vector
C
X i = 1
X GO TO 50
X
X 40 i = i + 1
X 50 IF (i.GT.nptwi) GO TO 60
X weight(i) = ysort(2,i)
X GO TO 40
C
C If QUANT is less than one-half of WEIGHT(1),
C extrapolate smoothed Y value for window.
C
X 60 IF (.NOT. (quant.LT.weight(1)/2.0)) GO TO 70
X yval = ysort(1,1) - (ysort(1,2)-ysort(1,1))*
X + (weight(1)/2.0-quant)/ ((weight(2)+weight(1))/2.0)
X GO TO 110
C
C
C Sum weights until sum is greater than QUANT or
C all weights are summed.
C
C SUM-WEIGHTS
X 70 ASSIGN 80 TO i99991
X GO TO 120
C
X 80 IF (.NOT. (qsmgtq)) GO TO 90
C
C Interpolate smoothed Y value
C
X yval = ysort(1,i) - (ysort(1,i)-ysort(1,i-1))* (sum-quant)/
X + ((weight(i)+weight(i-1))/2.0)
X GO TO 100
C
C All weights have been summed and sum is not greater than
C QUANT, therefore, extrapolate smoothed Y value.
C
X 90 yval = y(nptwi) + (quant-sum)/ ((weight(nptwi)+weight(nptwi-1))/
X + 2.0)* (ysort(1,nptwi)-ysort(1,nptwi-1))
X 100 CONTINUE
X 110 RETURN
C
C
X STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***'
C TO SUM-WEIGHTS
C
C Sum up weights until sum is greater than QUANT or
C all weights have been summed.
C
X 120 sum = 0.0
X qsmgtq = .FALSE.
C
X j = 1
X GO TO 140
X
X 130 j = j + 1
X 140 IF (j.GT.nptwi) GO TO 160
X i = j
X sum = sum + weight(j)/2.0
C
X IF (.NOT. (sum.GT.quant)) GO TO 150
X qsmgtq = .TRUE.
X GO TO i99991
X
X 150 sum = sum + weight(j)/2.0
C
X GO TO 130
X
X 160 GO TO i99991
C
X END
X SUBROUTINE est(x,y,nptwi,xval,estpar,dermth,estmet,xwgts,sortar,
X + yval,ixval,istats)
C*****************************************************************
C SUBROUTINE EST(X, Y, NPTWI, XVAL, ESTPAR, DERMTH, ESTMET, XWGTS,
C # SORTAR, YVAL, IXVAL, ISTATS)
C
C ESTimate value at window center
C
C
C Function
C
C
C Routine which calls specific estimation routines for each
C estimation method. Returns a smoothed Y value for Y values
C in window.
C
C**********************************************************************
C
C
C Arguments
C
C
C X --> Vector of X values in window. X(1) is lowest element
C in window.
C REAL X(NPTWI)
C
C Y --> Vector of Y values in window. Y(1) is Y value for
C lowest X element in window.
C REAL Y(NPTWI)
C
C NPTWI --> Number of points in current window.
C INTEGER NPTWI
C
C XVAL --> Window center
C REAL XVAL
C
C ESTPAR --> Parameter of the estimation method. Only needed and
C used for quantile estimation, in which case it is a
C quantile, i.e., a number between 0 and 1.
C REAL ESTPAR
C
C DERMTH --> Method to estimating derivative. Only needed if
C derivative estimation is to be done
C 1 - Linear
C 2 - Quadratic
C INTEGER DERMTH
C
C ESTMET --> Estimation method.
C 1 - Mean
C 2 - Linear
C 3 - Quadratic
C 4 - Maximum
C 5 - Minimum
C 6 - Quantile
C 7 - Standard deviation
C 8 - Derivative
C 9 - Window width
C 10 - Number of points
C INTEGER ESTMET
C
C XWGTS --> Weights of X values in current window (parallels X).
C REAL XWGTS(NPTWI)
C
C SORTAR --> Array for sorting Y values and X weights (parallels
C X). Used only for quantile estimation method.
C REAL SORTAR(2,NPTWI)
C
C YVAL <-- Smoothed Y value for Y values in window.
C REAL YVAL
C
C IXVAL --> X(IXVAL) will not be used in the estimation when cross
C validating. If not cross validating IXVAL is negative.
C INTEGER IXVAL
C
C ISTATS <-- 0 if an estimate for Y in this window was found.
C Greater than 0 if no estimate was found.
C Only of concern if the method of estimation is either
C linear or quadric fit.
C INTEGER ISTATS
C
C***********************************************************************
X REAL x(nptwi),y(nptwi),xval,estpar,xwgts(nptwi)
X REAL sortar(2,nptwi),yval
X INTEGER nptwi,dermth,estmet,ixval,istats
C
C**********************************************************************
C
C
C Parameters
C
C
C**********************************************************************
X INTEGER cmax,cmean,cmin,cnpt,cder
X INTEGER cquad,cquant,cstdev,cwnwid,clin
X PARAMETER (cmean=1,clin=2,cquad=3,cmax=4,cmin=5)
X PARAMETER (cquant=6,cstdev=7,cder=8,cwnwid=9,cnpt=10)
C
C***********************************************************************
C
C
C Main body of program
C
C
X istats = 0
X IF ((cmean).NE. (estmet)) GO TO 10
X CALL esmean(y,nptwi,xwgts,yval,ixval)
X GO TO 110
C
X 10 IF ((cquant).NE. (estmet)) GO TO 20
X CALL esquan(y,xwgts,nptwi,estpar,sortar,yval)
X GO TO 110
C
X 20 IF ((clin).NE. (estmet)) GO TO 30
X CALL eslin(x,y,xwgts,nptwi,xval,yval,ixval,istats)
X GO TO 110
C
X 30 IF ((cquad).NE. (estmet)) GO TO 40
X CALL esquad(x,y,xwgts,nptwi,xval,yval,ixval,istats)
X GO TO 110
C
X 40 IF ((cmax).NE. (estmet)) GO TO 50
X CALL estmax(y,nptwi,yval)
X GO TO 110
C
X 50 IF ((cmin).NE. (estmet)) GO TO 60
X CALL estmin(y,nptwi,yval)
X GO TO 110
C
X 60 IF ((cstdev).NE. (estmet)) GO TO 70
X CALL estdev(y,xwgts,nptwi,yval)
X GO TO 110
C
X 70 IF ((cnpt).NE. (estmet)) GO TO 80
X yval = nptwi
X GO TO 110
C
X 80 IF ((cwnwid).NE. (estmet)) GO TO 90
X CALL estwid(x,nptwi,yval)
X GO TO 110
C
X 90 IF ((cder).NE. (estmet)) GO TO 100
X CALL drvtv(x,y,xwgts,xval,nptwi,dermth,yval,istats)
X GO TO 110
C
X 100 CONTINUE
X 110 RETURN
C
X END
X SUBROUTINE estdev(y,weight,nptwi,yval)
C*****************************************************************
C SUBROUTINE ESTDEV(Y, WEIGHT, NPTWI, YVAL)
C
C Estimate STandard DEViation
C
C
C Function
C
C
C Calculates the weighted standard deviation of the elements
C of the Y vector.
C
C**********************************************************************
C
C
C Arguments
C
C
C Y --> Vector of Y values in window.
C REAL Y(NPTWI)
C
C WEIGHT --> Weights associated with elements of Y vector.
C REAL WEIGHT(NPTWI)
C
C NPTWI --> Number of points in Y vector.
C INTEGER NPTWI
C
C YVAL <-- Standard deviation of the elements of Y.
C REAL YVAL
C
C*****************************************************************
X REAL y(nptwi),weight(nptwi),yval
X INTEGER nptwi
C**********************************************************************
C
C
C Variables
C
C
C YMEAN --- Weighted mean of elements of Y.
C REAL YMEAN
C
C DIFF --- Difference from mean for each Y value
C REAL DIFF
C
C SGMASQ --- Sum of squares of differences from mean for
C each Y value
C REAL SGMASQ
C
C WTSUM --- Sum of weights associated with elements of Y vector
C REAL WTSUM
C
C**********************************************************************
X REAL ymean,diff,sgmasq,wtsum
C
C*****************************************************************
C
C
C Main body of program
C
C Get the weighted mean of elements of Y.
C
C and since we want all the points used IYVAL is set negative
C
X iyval = -1
X CALL esmean(y,nptwi,weight,ymean,iyval)
C
C Calculate differences from mean for each Y value
C
X sgmasq = 0.0
C
X i = 1
X GO TO 20
X
X 10 i = i + 1
X 20 IF (i.GT.nptwi) GO TO 30
X diff = y(i) - ymean
X sgmasq = sgmasq + weight(i)*diff*diff
X GO TO 10
X
X 30 wtsum = 0.0
C
C Calculate the sum of the weights
C
C
X i = 1
X GO TO 50
X
X 40 i = i + 1
X 50 IF (i.GT.nptwi) GO TO 60
X wtsum = wtsum + weight(i)
X GO TO 40
X
X 60 yval = sqrt(sgmasq/wtsum)
C
C Calculate standard deviation
C
C
X RETURN
X
X END
X SUBROUTINE estmax(y,nptwi,yval)
C*****************************************************************
C SUBROUTINE ESTMAX(Y, NPTWI, YVAL)
C
C ESTimate by MAXimum
C
C
C Function
C
C
C Returns the algebraic maximum of the elements of the Y vector.
C
C
C**********************************************************************
C
C
C Arguments
C
C Y --> Vector of Y values.
C REAL Y(NPTWI)
C
C NPTWI --> Number of points in current window.
C INTEGER NPTWI
C
C YVAL <-- Maximum Y value for Y values in window.
C REAL YVAL
C
C*****************************************************************
X REAL y(nptwi),yval
X INTEGER nptwi
C
C**********************************************************************
C
C
C Variables
C
C
C YMIN <-- Minimum Y value for Y values in window.
C REAL YMIN
C
C**********************************************************************
X REAL ymin
C
C**********************************************************************
C
C
C Main body of program
C
C
X CALL minmax(y,nptwi,ymin,yval)
C
X RETURN
X
X END
X SUBROUTINE estmin(y,nptwi,yval)
C*****************************************************************
C SUBROUTINE ESTMIN(Y, NPTWI, YVAL)
C
C ESTimate by MINimum
C
C
C Function
C
C
C Returns the algebraic minimum of the elements of the Y vector.
C
C
C**********************************************************************
C
C
C Arguments
C
C Y --> Vector of Y values.
C REAL Y(NPTWI)
C
C NPTWI --> Number of points in current window.
C INTEGER NPTWI
C
C YVAL <-- Minimum Y value for Y values in window.
C REAL YVAL
C
C*****************************************************************
X REAL y(nptwi),yval
X INTEGER nptwi
C
C**********************************************************************
C
C
C Variables
C
C
C YMAX <-- Maximum Y value for Y values in window.
C REAL YMAX
C
C**********************************************************************
X REAL ymin
C
C**********************************************************************
C
C
C Main body of program
C
C
X CALL minmax(y,nptwi,yval,ymax)
C
X RETURN
X
X END
X SUBROUTINE estwid(x,nptwi,wid)
C*****************************************************************
C SUBROUTINE ESTWID(X, NPTWI, WID)
C
C ESTimate WIDth of window
C
C
C Function
C
C
C Calculates width of a window containing a specified number
C of points.
C
C**********************************************************************
C
C
C Arguments
C
C
C X --> Vector of X values in window. Sorted into ascending
C order upon input to a routine calling this one.
C REAL X(NPTWI)
C
C NPTWI --> Number of points in current window.
C INTEGER NPTWI
C
C WID <-- Width of window.
C REAL YVAL
C
C*****************************************************************
X REAL x(nptwi),wid
X INTEGER nptwi
C
C
C Main body of program
C
C
C Since X vector is sorted, it is possible to caluculate width
C of window by subtracting first element from last.
C
X wid = x(nptwi) - x(1)
C
X RETURN
X
X END
X INTEGER FUNCTION IBSRL(VAL,TAB,NTAB,QINCR,QVALLO)
C
C**********************************************************************
C
C INTEGER FUNCTION IBSRL(VAL,TAB,NTAB,QINCR,QVALLO)
C I Binary Search of a ReaL table
C
C
C Function
C
C
C Performs a binary search of array TAB to find the value VAL.
C Returns the index of the element equal to VAL if there is one,
C else the element bounding the value above or below depending
C on QVALLO.
C
C
C Arguments
C
C
C IBSRL <-- Returns the index of the element of TAB that equals
C VAL if there is one (exact floating point equality
C is required for this to happen). If not, returns the
C index of the element that immediately bounds VAL
C below or above depending on the value of QVALLO. If
C VAL is outside the range of values of TAB, a 0 or
C NTAB + 1 value will be returned depending on which
C side of TAB the value VAL lies.
C IBSSRL is INTEGER
C
C VAL --> The value to be found by the binary search.
C VAL is REAL
C
C TAB --> The table to be searched for the value VAL. Note that
C the elements of TAB must be either non-increasing or
C non-decreasing according to the setting of QINCR;
C otherwise the behavior of the routine is unpredictable.
C TAB is REAL TAB(NTAB)
C
C NTAB --> The length of TAB.
C NTAB is INTEGER
C
C QINCR --> .TRUE. if the elements of TAB increase with their index
C (or at least do not decrease). .FALSE. if the elements
C of TAB decrease with their index (or at least do not
C increase).
C QINCR is LOGICAL
C
C QVALLO --> .TRUE. if the index returned by IBSRL should have the
C property that VAL .le. TAB(IBSRL); .FALSE. if the
C property is VAL .ge. TAB(IBSRL).
C QVALLO is LOGICAL
C
C**********************************************************************
C
C
C
C CASE IN WHICH TAB IS INCREASING (MONOTONE NON-DECREASING)
C
C
C .. Scalar Arguments ..
X REAL VAL
X INTEGER NTAB
X LOGICAL QINCR,QVALLO
C ..
C .. Array Arguments ..
X REAL TAB(NTAB)
C ..
C .. Local Scalars ..
X INTEGER IHI,ILO,IMID
C ..
C .. Executable Statements ..
X IF (.NOT.(QINCR)) GOTO 99998
X IF (QVALLO .AND. (VAL .GT. TAB(NTAB))) THEN
X IBSRL = NTAB + 1
X ELSE IF (.NOT.QVALLO .AND. (VAL .LT. TAB(1))) THEN
X IBSRL = 0
X ELSE
X ILO = 1
X IHI = NTAB
X GOTO 99997
X99996 IF ((IHI-ILO).LE.1) GOTO 99995
X99997 IMID = (ILO+IHI)/2
X IF (.NOT.(VAL.LE.TAB(IMID))) GOTO 99993
X IHI = IMID
X GOTO 99994
X99993 ILO = IMID
X99994 GOTO 99996
X99995 IF (.NOT.((QVALLO .AND. (VAL .LE. TAB(ILO))) .OR. (.NOT.QVALLO .AN
X 1D. (VAL .LT. TAB(IHI))))) GOTO 99991
X IBSRL = ILO
X GOTO 99992
X99991 IBSRL = IHI
X99992 ENDIF
X GOTO 99999
C
C
C CASE IN WHICH TAB IS DECREASING (MONOTONE NON-INCREASING)
C
C
X99998 IF (QVALLO .AND. (VAL .GT. TAB(1))) THEN
X IBSRL = 0
X ELSE IF (.NOT.QVALLO .AND. (VAL .LT. TAB(NTAB))) THEN
X IBSRL = NTAB + 1
X ELSE
X ILO = 1
X IHI = NTAB
X GOTO 99990
X99989 IF ((IHI-ILO).LE.1) GOTO 99988
X99990 IMID = (ILO+IHI)/2
X IF (.NOT.(VAL.LE.TAB(IMID))) GOTO 99986
X ILO = IMID
X GOTO 99987
X99986 IHI = IMID
X99987 GOTO 99989
X99988 IF (.NOT.((QVALLO .AND. (VAL .LE. TAB(IHI))) .OR. (.NOT.QVALLO .AN
X 1D. (VAL .LT. TAB(ILO))))) GOTO 99984
X IBSRL = IHI
X GOTO 99985
X99984 IBSRL = ILO
X99985 ENDIF
X99999 RETURN
X END
X INTEGER FUNCTION IDAMAX(N,DX,INCX) 00524
X DOUBLE PRECISION DX(1),DMAX 00529
X INTEGER I,INCX,IX,N 00530
X IDAMAX = 0 00532
X IF (N .LT. 1) RETURN 00533
X IDAMAX = 1 00534
X IF (N.EQ.1) RETURN 00535
X IF (INCX.EQ.1) GO TO 20 00536
X IX = 1 00540
X DMAX = DABS(DX(1)) 00541
X IX = IX + INCX 00542
X DO 10 I = 2,N 00543
X IF (DABS(DX(IX)).LE.DMAX) GO TO 5 00544
X IDAMAX = I 00545
X DMAX = DABS(DX(IX)) 00546
X 5 IX = IX + INCX 00547
X 10 CONTINUE 00548
X RETURN 00549
X 20 DMAX = DABS(DX(1)) 00553
X DO 30 I = 2,N 00554
X IF (DABS(DX(I)).LE.DMAX) GO TO 30 00555
X IDAMAX = I 00556
X DMAX = DABS(DX(I)) 00557
X 30 CONTINUE 00558
X RETURN 00559
X END 00560
X SUBROUTINE lqbeta(x,y,weight,n,betas,nbeta,ixval,info)
C***********************************************************************
C Function
C
C
C Does a weighted fit to X,Y data. Returns the coeficients
C of regression.
C
C
C**********************************************************************
C
C
C Method
C
C The sum of the squares of the deviations of the input data points
C from a fitted curve can be minimized by setting the partial
C derivative of that sum with respect to each parameter of the curve
C to zero. This yields as many linear equations as there are
C parameters. For weighted fit, the deviation of each data value
C fit, the deviation of each data value from the fitted curve is
C multiplied by the weight for that value. The equations can be
C manipulated to yield expressions for the curve parameters in terms
C of summations involving various powers and products of the X, Y
C and WEIGHT values. These expressions are coded in this routine.
C The summations are evaluated, then the curve parameters are
C calculated. no sort on X or Y is needed.
C
C**********************************************************************
C
C
C Arguments
C
C
C X --> Vector of X values.
C REAL X(N)
C
C Y --> Vector of Y values.
C REAL Y(N)
C
C WEIGHT --> Weights of X values.
C REAL WEIGHT(N)
C
C N --> Number of points.
C INTEGER N
C
C BETAS <-- Coeficients of regression.
C Returned only if INFO equals 0.
C REAL BETAS(NBETA)
C
C NBETA --> Number of parameters. 2 for linear regression
C 3 for quadractic regression
C INTEGER NBETA
C
C IXVAL --> All points but X(IXVAL) are used in the fit. If all
C points are to be used IXVAL is negative.
C INTEGER IXVAL
C
C INFO <-- Information about the sucsess of this routine.
C If sucsessful, INFO equals 0
C Otherwise INFO is greater than 0.
C INTEGER INFO
C***********************************************************************
X IMPLICIT INTEGER (i-n),DOUBLE PRECISION (a-h,o-p,r-z),LOGICAL (q)
C
X REAL x(n),y(n),weight(n),betas(nbeta)
X INTEGER n,nbeta,ixval,info
C
X INTEGER ipvt(10)
X INTEGER bmsize
X PARAMETER (bmsize=3)
X REAL small
X PARAMETER (small=1E-10)
X INTEGER maxn
X PARAMETER (maxn=500)
X DIMENSION det(2),xwx(bmsize,bmsize),xwy(bmsize),work(bmsize)
X DIMENSION dx(maxn),dy(maxn)
C
C Initialize sums to zero
C
X sumw = 0.0D0
X sumwy = 0.0D0
X sumwx = 0.0D0
X sumwxy = 0.0D0
X sumwx2 = 0.0D0
X smwx2y = 0.0D0
X sumwx3 = 0.0D0
X sumwx4 = 0.0D0
C
C DO WEIGHTED SUMS TO GET THE COEFFICIENTS FOR THE PARAMETERS
C OF THE FUNCTION.
C
X DO 30,i = 1,n
X IF (.NOT. ((i.NE.ixval).AND. (weight(i).GT.
X + small))) GO TO 20
X dx(i) = dble(x(i))
X dy(i) = dble(y(i))
X tempw = 1.0D0/dble(weight(i))
X tempwx = tempw*dx(i)
X sumw = sumw + tempw
X sumwy = sumwy + tempw*dy(i)
X sumwx = sumwx + tempwx
X sumwxy = sumwxy + tempwx*dy(i)
X sumwx2 = sumwx2 + tempwx*dx(i)
C
C If we are fitting a quadratic,
C then finish getting the weighted sums
C
X IF (.NOT. (nbeta.GT.2)) GO TO 10
X smwx2y = smwx2y + tempwx*dx(i)*dy(i)
X sumwx3 = sumwx3 + tempwx*dx(i)*dx(i)
X sumwx4 = sumwx4 + tempwx*dx(i)*dx(i)*dx(i)
C
X 10 CONTINUE
X 20 CONTINUE
X 30 CONTINUE
X xwx(1,1) = sumw
C
C Build matrix XWX (X-trans * W-inv * X)
C
X xwx(1,2) = sumwx
X xwx(1,3) = sumwx2
X xwx(2,1) = sumwx
X xwx(2,2) = sumwx2
X xwx(2,3) = sumwx3
X xwx(3,1) = sumwx2
X xwx(3,2) = sumwx3
X xwx(3,3) = sumwx4
C
C Build matrix XWY (X-trans * W-inv * Y)
C
X xwy(1) = sumwy
X xwy(2) = sumwxy
X xwy(3) = smwx2y
C
C
C Find inverse of (X-trans * w-inv * X)
C
X CALL dgefa(xwx,bmsize,nbeta,ipvt,info)
X job = 1
X IF (info.EQ.0) CALL dgedi(xwx,bmsize,nbeta,ipvt,det,work,job)
C
C Done Inverting -- Multiply XWX-inverse by XWY
C
X IF (.NOT. (info.EQ.0)) GO TO 60
X DO 50,i = 1,nbeta
X sum = 0.0D0
X DO 40,j = 1,nbeta
X sum = sum + xwx(i,j)*xwy(j)
X 40 CONTINUE
X betas(i) = real(sum)
X 50 CONTINUE
X 60 RETURN
X
X END
X SUBROUTINE minmax(y,ny,ymin,ymax)
C***********************************************************
C SUBROUTINE MINMAX(Y, NY, YMIN, YMAX)
C
C return MIN and MAX of real vector
C
C
C Function
C
C
C Finds minimum and maximum values of a vector.
C
C**********************************************************************
C
C
C Method
C
C The vector Y is scanned from beginning to end. Each element is
C compared to previous maximum and minimum.
C
C*****************************************************************
C
C
C Arguments
C
C Y --> Vector to be searched for minimum and maximum.
C REAL Y(NY)
C
C NY --> Number of points Y vector.
C INTEGER NY
C
C YMIN <-- Minimum value in Y.
C REAL YMIN
C
C YMAX <-- Maximum value in Y.
C REAL YMAX
C
C*****************************************************************
X REAL y(ny),ymin,ymax
X INTEGER ny
C
C**********************************************************************
C
C
C Main body of program
C
C
X ymax = y(1)
X ymin = y(1)
C
X i = 2
X GO TO 20
X
X 10 i = i + 1
X 20 IF (i.GT.ny) GO TO 30
X IF (y(i).GT.ymax) ymax = y(i)
X IF (y(i).LT.ymin) ymin = y(i)
X GO TO 10
X
X 30 RETURN
C
X END
X LOGICAL FUNCTION qgr(a,b,irow)
C*****************************************************************
C LOGICAL FUNCTION QGR(A, B, IROW)
C
C Q a GReater than b.
C
C FUNCTION
C
C Returns true value if A>B, else false. Called by COMPLIB
C routine SRMRL.
C
C ARGUMENTS
C
C A,B --> REAL - Values to be compared.
C IROW --> INTEGER - Where A and B are vectors, IROW is the
C element of the vectors which is to be compared. Not
C used here.
C*****************************************************************
C DECLARATIONS & SPECIFICATIONS
C
X IMPLICIT LOGICAL (q)
C
X IF (.NOT. (a.GT.b)) GO TO 10
X qgr = .TRUE.
X GO TO 20
X
X 10 qgr = .FALSE.
C
X 20 RETURN
X
X END
X SUBROUTINE SRMRL(QGR,A,NROWUS,NROWDM,NCOL,IROW)
C**********************************************************************
C
C SUBROUTINE SRMRL(QGR,A,NROWUS,NROWDM,NCOL,IROW)
C
C SoRt Matrix of ReaL values
C
C
C Function
C
C
C Sorts the columns of array A in ascending order acccording
C to the collating sequence implied by the user supplied
C function QGR.
C Note that each element of A must be less than or equal to
C 128 characters in length for the sort to work. Otherwise the
C routine will exit with an error indication.
C
C
C Arguments
C
C
C QGR --> The name of a logical function supplied by the user.
C It has three arguments, X, Y, and IROW, e.g.,
C QXGTY = QGR(X,Y,IROW)
C where X and Y are columns of array A and must be
C dimensioned accordingly. IROW is passed through to QGR
C unchanged from the value sent to this routine.
C QGR should return .TRUE. if X is greater than Y in the
C collating sequence desired and .FALSE. if X is less than
C or equal to Y.
C QGR must be declared to be EXTERNAL in the calling
C program.
C QGR is LOGICAL
C
C A --> The array whose columns are sorted.
C A is REAL A(NROWDM,NCOL)
C
C NROWUS --> The number of rows of A that are actually used. The
C sort routine works by permuting columns and only rows
C 1..NROWUS of each column are permuted. The remaining
C rows are not referenced by this routine.
C If a singly dimensioned array is being sorted
C then NROWUS should be set to 1.
C NROWUS is INTEGER
C
C NROWDM --> The number of rows declared in A in the program that
C dimensions it.
C If a singly dimensioned array is being sorted
C then NROWDM should be set to 1.
C NROWDM is INTEGER
C
C NCOL --> The number of columns of A to be sorted.
C Columns 1..NCOL are sorted. It is irrelevent
C to this routine what the number of columns declared
C in A is.
C NCOL is INTEGER
C
C IROW --> An integer argument passed unchanged to QGR. The intent
C is that it be the row on which the sort is performed in
C cases in which that makes sense. This allows one version
C of QGR to be used regardless of the row on which the
C sort is performed. Note that IROW can be used for
C any purpose that the user desires.
C IROW is INTEGER
C
C
C Note
C
C
C The algorithm used is based on the article by Robert Sedgewick,
C "Implementing Quicksort Programs", CACM, Oct 1978, 21:10,847-857.
C
C If QGR allows ties, the order of the tied columns is not
C determined by this routine.
C
C**********************************************************************
C***
C
C DECLARE ARGUMENTS
C
C***
C
C *** VARIABLES
C
C
C ** PARAMETERS AND CONSTANTS
C
C * CONSTANTS
C
C QFALSE -- MNEMONICALLY-NAMED LOGICAL CONSTANT
C
C
C QTRUE -- MNEMONICALLY-NAMED LOGICAL CONSTANT
C
C
C * PARAMETERS PECULIAR TO THIS APPLICATION
C
C CHGSRT -- THRESHOLD ABOVE WHICH QUICKSORT IS FASTER THAN
C INSERTION SORT
C
C
C STCKMX -- MAXIMUM NUMBER OF ENTRIES IN STACK
C
C
C
C ** OTHER SIGNIFICANT VARIABLES
C
C PART -- INDEX IN ARRAY OF ELEMENT CHOSEN AS PARTITION
C
C
C
C STCK -- STACK TO SAVE THE STARTING AND ENDING POSITIONS OF THE
C PARTITIONS
C
C
C STCKCT -- CURRENT LENGTH OF STACK
C
C
C SWAP* -- INDICES IN ARRAY OF ELEMENTS TO BE SWAPPED
C
C
C
C ** MISCELLANEOUS VARIABLES
C THE FOLLOWING VARIABLES ARE USED FOR SUCH PURPOSES AS INDICES AND
C TEMPORARY STORAGE. EACH ONE IS INTENDED TO BE SIGNIFICANT ONLY
C WITHIN A WELL-DEFINED BLOCK OF CODE (ALTHOUGH IT MAY APPEAR IN
C MORE THAN ONE BLOCK), AND ITS MEANING AND USE SHOULD BE APPARENT
C FROM ITS NAME AND CONTEXT.
C
C I --
C
C J --
C
C L --
C
C Q --
C
C QDONE --
C
C R --
C
C ICOL --
C
C II --
C
C JJ --
C
C IXROW --
C
C
C
C
C *** FUNCTIONS AND SUBROUTINES
C
C
C ** LIBRARY SUBPROGRAMS
C DETAILED INFORMATION ON EACH OF THE FOLLOWING SUBPROGRAMS CAN BE
C FOUND IN THE REFERENCE FOR ITS RESPECTIVE LIBRARY.
C
C * FORTRAN
C INFORMATION ON EACH OF THE FOLLOWING CAN BE FOUND IN THE FORTRAN
C LANGUAGE REFERENCE MANUAL.
C
C MAX -- MAXIMUM
C
C
C MIN -- MINIMUM
C
C
C LEN -- LENGTH OF A CHARACTER STRING
C
C
C
C
C
C *** QUICKSORT
C .. Parameters ..
X LOGICAL QFALSE
X PARAMETER (QFALSE=.FALSE.)
X LOGICAL QTRUE
X PARAMETER (QTRUE=.TRUE.)
X INTEGER CHGSRT
X PARAMETER (CHGSRT=10)
X INTEGER STCKMX
X PARAMETER (STCKMX=50)
C ..
C .. Scalar Arguments ..
X INTEGER IROW,NCOL,NROWDM,NROWUS
C ..
C .. Array Arguments ..
X REAL A(NROWDM,NCOL)
C ..
C .. Function Arguments ..
X LOGICAL QGR
X EXTERNAL QGR
C ..
C .. Local Scalars ..
X REAL TEMP
X INTEGER I,I99997,ICOL,IXROW,J,L,PART,R,STCKCT,SWAP1,SWAP2
X LOGICAL Q,QDONE
C ..
C .. Local Arrays ..
X INTEGER STCK(2,STCKMX)
C ..
C .. Intrinsic Functions ..
X INTRINSIC MAX,MIN
C ..
C .. Executable Statements ..
X STCKCT = 0
X L = 1
X R = NCOL
X QDONE = (R - L) .LT. CHGSRT
X99999 IF (QDONE) GOTO 99998
C
C
C *** 'MEDIAN-OF-THREE' MODIFICATION. SEE 'OVERVIEW' ABOVE
C
X SWAP1 = (L+R)/2
X SWAP2 = L + 1
C SWAP-ELEMENTS
X ASSIGN 99996 TO I99997
X GOTO 99997
C
X99996 IF (.NOT.(QGR(A(1,L+1),A(1,R),IROW))) GOTO 99995
X SWAP1 = L + 1
X SWAP2 = R
C SWAP-ELEMENTS
X ASSIGN 99994 TO I99997
X GOTO 99997
X99994 CONTINUE
X99995 IF (.NOT.(QGR(A(1,L),A(1,R),IROW))) GOTO 99993
X SWAP1 = L
X SWAP2 = R
C SWAP-ELEMENTS
X ASSIGN 99992 TO I99997
X GOTO 99997
X99992 CONTINUE
X99993 IF (.NOT.(QGR(A(1,L+1),A(1,L),IROW))) GOTO 99991
X SWAP1 = L + 1
X SWAP2 = L
C SWAP-ELEMENTS
X ASSIGN 99990 TO I99997
X GOTO 99997
X99990 CONTINUE
X99991 I = L + 1
C
C
C *** PARTITIONING OF SUBFILE
C
X J = R
X PART = L
C
X99988 GOTO 99987
X99986 IF (.NOT.(QGR(A(1,PART),A(1,I),IROW))) GOTO 99985
X99987 I = I + 1
X GOTO 99986
X99985 GOTO 99984
X99983 IF (.NOT.(QGR(A(1,J),A(1,PART),IROW))) GOTO 99982
X99984 J = J - 1
X GOTO 99983
C
X99982 IF (.NOT.(J .LT. I)) GOTO 99980
X GOTO 99989
X GOTO 99981
X99980 SWAP1 = I
X SWAP2 = J
C SWAP-ELEMENTS
X ASSIGN 99979 TO I99997
X GOTO 99997
X99979 CONTINUE
X99981 GOTO 99988
C
X99989 SWAP1 = L
C
X SWAP2 = J
C SWAP-ELEMENTS
X ASSIGN 99978 TO I99997
X GOTO 99997
C
C
C *** RECURSION STEP
C
X99978 IF (.NOT.(MAX((J-L),(R-I+1)) .LE. CHGSRT)) GOTO 99976
C *** DEFER FURTHER SORTING OF BOTH SUBFILES TO INSERTION SORT
C
X IF (.NOT.(STCKCT .EQ. 0)) GOTO 99974
X QDONE = QTRUE
X GOTO 99975
C *** POP STACK AND CONTINUE QUICKSORT
X99974 L = STCK(1,STCKCT)
X R = STCK(2,STCKCT)
X STCKCT = STCKCT - 1
X99975 GOTO 99977
C
C *** CONTINUE QUICKSORT ON AT LEAST ONE SUBFILE
C
X99976 IF (.NOT.(MIN((J-L),(R-I+1)) .LE. CHGSRT)) GOTO 99972
C *** DEFER SMALL TO INSERTION SORT, CONTINUE QUICKSORT ON LARGE
C
X IF (.NOT.((J-L) .GE. (R-I+1))) GOTO 99970
C *** LEFT SUBFILE IS LARGE ONE
C L = L
X R = J - 1
X GOTO 99971
C *** RIGHT SUBFILE IS LARGE ONE
X99970 L = I
C R = R
X99971 GOTO 99973
C
C *** CONTINUE QUICKSORT ON BOTH SUBFILES
C
X99972 IF (.NOT.(STCKCT .GE. STCKMX)) GOTO 99969
C *** STACK IS FULL
X STOP ' STACK OVERFLOW IN SRMRL'
X99969 STCKCT = STCKCT + 1
C
C *** PUSH LARGE SUBFILE ONTO STACK, CONTINUE QUICKSORT WITH SMALL
X IF (.NOT.((J-L) .GE. (R-I+1))) GOTO 99967
C *** LEFT SUBFILE IS LARGE ONE
C
X STCK(1,STCKCT) = L
X STCK(2,STCKCT) = J - 1
X L = I
C R = R
C
X GOTO 99968
C *** RIGHT SUBFILE IS LARGE ONE
C
X99967 STCK(1,STCKCT) = I
X STCK(2,STCKCT) = R
C L = L
X R = J - 1
C
C
X99968 CONTINUE
C
X99973 CONTINUE
X99977 GOTO 99999
C
C
X99998 CONTINUE
X1000 DO 99966, I=(NCOL-1),1,-1
C
C
C *** INSERTION SORT
C
C *** FOR EACH POSITION PRECEDING THE LAST POSITION ...
C
X IF (.NOT.(QGR(A(1,I),A(1,I+1),IROW))) GOTO 99964
C *** A(I) NEEDS TO BE MOVED FURTHER OUT
C
X J = I + 1
X GOTO 99963
X99962 IF (.NOT.(Q)) GOTO 99961
C
X99963 J = J + 1
C
X IF (.NOT.(J .GT. NCOL)) GOTO 99959
X Q = QFALSE
X GOTO 99960
X99959 Q = QGR(A(1,I),A(1,J),IROW)
C
X99960 GOTO 99962
X99961 DO 99958, IXROW=1,NROWUS
X TEMP = A(IXROW,I)
X DO 99956, ICOL=I,J-2
X A(IXROW,ICOL) = A(IXROW,ICOL+1)
X99956 CONTINUE
X A(IXROW,J-1) = TEMP
X99958 CONTINUE
C
X99964 CONTINUE
X99966 CONTINUE
C
X RETURN
X STOP '*** EXECUTION FLOWING INTO FLECS PROCEDURES ***'
C TO SWAP-ELEMENTS
X99997 DO 99954, IXROW=1,NROWUS
X TEMP = A(IXROW,SWAP1)
X A(IXROW,SWAP1) = A(IXROW,SWAP2)
X A(IXROW,SWAP2) = TEMP
X99954 CONTINUE
X GOTO I99997
C
C
X END
X SUBROUTINE window(x,y,nxy,xdes,ysmo,ndes,method,smparm,dermth,
X + wtmth,estmth,qtle,sortar,xwgts,xmiss,qcross,
X + istat)
C*********************************************************************
C SUBROUTINE WINDOW(X,Y,NXY,XDES,YSMO,NDES,METHOD,SMPARM,DERMTH
C # WTMTH,ESTMTH,QTLE,SORTAR,XWGTS,XMISS,QCROSS,ISTAT)
C
C WINDOW subroutine for S
C
C Function
C
C Returns smooth and continuous data representing the noisy or
C binary data input in the X and Y arrays.
C
C***********************************************************************
C
C Arguments
C
C X --> Vector of X values (parallels Y). Must be sorted
C in ascending order on input.
C REAL X(NXY)
C
C Y --> Vector of Y values (parallels X).
C REAL Y(NXY)
C
C NXY --> Length of X and Y arrays.
C INTEGER NXY
C
C XDES --> Vector of X values at which smoothed Y values are
C desired. (Parallels YSMO)
C REAL XDES(NDES)
C
C YSMO --> Vector of smoothed Y values (parallels Y).
C REAL YSMO(NDES)
C
C NDES --> Length of XDES and YSMO arrays.
C INTEGER NDES
C
C METHOD --> Windowing method.
C 1 - Fixed width
C 2 - Nearest neighbor
C INTEGER METHOD
C
C SMPARM --> Method-dependent windowing parameter.
C If METHOD = 1, SMPARM is the width of the fixed
C width window.
C If METHOD = 2, SMPARM is the number of nearest
C neighbors to be included in the window.
C REAL SMPARM
C
C DERMTH --> Method to estimating derivative. Only needed if
C derivative estimation is to be done
C 1 - Linear
C 2 - Quadratic
C INTEGER DERMTH
C
C WTMTH --> Weight method.
C 1 - Boxcar
C 2 - Biquadratic
C INTEGER WTMTH
C
C ESTMTH --> Estimation method.
C 1 - Mean
C 2 - Linear
C 3 - Quadratic
C 4 - Maximum
C 5 - Minimum
C 6 - Quantile
C 7 - Standard deviation
C 8 - Derivative
C 9 - Window width
C 10 - Number of points
C INTEGER ESTMTH
C
C QTLE --> Parameter of the estimation method. Only needed and
C used for quantile estimation, in which case it is a
C quantile, i.e. a number between 0 and 1.
C REAL QTLE
C
C SORTAR --> Array for sorting Y values by X weights.
C Used only for quantile estimation method.
C REAL SORTAR(2,NXY)
C
C XMISS <-- Vector which contains the indicies ofthose X values
C that yield an empty window.
C INTEGER XMISS(NXY)
C
C QCROSS --> .TRUE. if cross valdiation is to be done
C LOGICAL QCROSS
C
C ISTAT <-- 0 if everything went ok, else positive
C INTEGER ISTAT
C
C**********************************************************************
C
X REAL x(nxy),xwgts(nxy),y(nxy),xdes(ndes),ysmo(ndes),xmiss(nxy)
X REAL sortar(2,nxy)
X INTEGER nxy,ndes
X REAL smparm,qtle
X INTEGER method,dermth,wtmth,estmth
X LOGICAL qcross
C***********************************************************************
C
C
C Variables
C
C
C IXLO -- Index of first element in X array in current window.
C INTEGER IXLO
C
C NPTWI -- Nunber of points in current window.
C INTEGER NPTWI
C
C NTIELO -- Number of X elements on left boundary of window.
C INTEGER NTIELO
C
C NTIEHI -- Number of X elements on right boundary of window.
C INTEGER NTIEHI
C
C QEMPTY -- .TRUE. if current window is empty, otherwise .FALSE.
C LOGICAL QEMPTY
C
C***********************************************************************
X INTEGER ixlo,nptwi,ntielo,ntiehi
X LOGICAL qempty
X REAL sum
C
C***********************************************************************
C
C Parameters
C
C
C***********************************************************************
X INTEGER cmax,cmean,cmin,cnpt,cder
X INTEGER cquad,cquant,cstdev,cwnwid,clin
X PARAMETER (cmean=1,clin=2,cquad=3,cmax=4,cmin=5)
X PARAMETER (cquant=6,cstdev=7,cder=8,cwnwid=9,cnpt=10)
C
C***********************************************************************
C
C
C Main body of program
C
C
X nempty = 0
X i = 1
X GO TO 20
X
X 10 i = i + 1
X 20 IF (i.GT.ndes) GO TO 230
C
C Iterate over windows
C
C Get highest and lowest indices of X array such that all elements
C of X with indices between these bounds are in window
C
X CALL windo(x,nxy,xdes(i),method,smparm,ixlo,nptwi,ntielo,ntiehi)
C
C Check for empty window
C
X qempty = .FALSE.
X IF ((cmean).NE. (estmth)) GO TO 30
X IF (nptwi.LT.1) qempty = .TRUE.
X GO TO 130
X
X 30 IF ((clin).NE. (estmth)) GO TO 40
X IF (nptwi.LT.2) qempty = .TRUE.
X GO TO 130
X
X 40 IF ((cquad).NE. (estmth)) GO TO 50
X IF (nptwi.LT.3) qempty = .TRUE.
X GO TO 130
X
X 50 IF ((cmax).NE. (estmth)) GO TO 60
X IF (nptwi.LT.1) qempty = .TRUE.
X GO TO 130
X
X 60 IF ((cmin).NE. (estmth)) GO TO 70
X IF (nptwi.LT.1) qempty = .TRUE.
X GO TO 130
X
X 70 IF ((cquant).NE. (estmth)) GO TO 80
X IF (nptwi.LT.2) qempty = .TRUE.
X GO TO 130
X
X 80 IF ((cstdev).NE. (estmth)) GO TO 90
X IF (nptwi.LT.1) qempty = .TRUE.
X GO TO 130
X
X 90 IF ((cder).NE. (estmth)) GO TO 100
X IF (dermth.EQ.1 .AND. nptwi.LT.2) qempty = .TRUE.
X IF (dermth.EQ.2 .AND. nptwi.LT.3) qempty = .TRUE.
X GO TO 130
X
X 100 IF ((cwnwid).NE. (estmth)) GO TO 110
X IF (nptwi.LT.2) qempty = .TRUE.
X GO TO 130
X
X 110 IF ((cnpt).NE. (estmth)) GO TO 120
X IF (nptwi.LT.2) qempty = .TRUE.
X GO TO 130
X
X 120 CONTINUE
X 130 IF (.NOT. (qempty)) GO TO 140
C
C Count empty windows and record indices of empty windows.
C
X nempty = nempty + 1
X xmiss(nempty) = i
X GO TO 220
C
C If window is non-empty calculate weights of points in window
C and estimate Y value for window.
C
X 140 CALL wgt(x(ixlo),nptwi,wtmth,method,xdes(i),smparm,ntielo,ntiehi,
X + xwgts(ixlo))
C
X IF (.NOT. (qcross)) GO TO 170
X ixval = i - ixlo + 1
C
C Adjust the weights so they will sum to one without X(IXVAL)
C
X sum = 0.0
X DO 150,ii = ixlo,ixlo + nptwi - 1
X IF (ii.NE.ixlo+ixval-1) sum = sum + xwgts(ii)
X 150 CONTINUE
X DO 160,ii = ixlo,ixlo + nptwi - 1
X xwgts(ii) = xwgts(ii)/sum
X 160 CONTINUE
X GO TO 180
X
X 170 ixval = -1.0
X 180 CALL est(x(ixlo),y(ixlo),nptwi,xdes(i),qtle,dermth,estmth,
X + xwgts(ixlo),sortar,ysmo(i),ixval,istat)
X IF (.NOT. (istat.NE.0)) GO TO 210
X IF (.NOT. (qcross)) GO TO 190
X RETURN
X
X GO TO 200
X
X 190 nempty = nempty + 1
X xmiss(nempty) = i
X 200 CONTINUE
X 210 CONTINUE
X 220 GO TO 10
X
X 230 RETURN
C
X END
X SUBROUTINE wdthmx(x,n,npts,smlwnd)
C***********************************************************************
C
C This subroutine returns the smallest window size that includes
C NPTS points.
C
C ARGUEMENTS
C
C X --> The real vector of data points.
C REAL X(N)
C
C N --> The number of values in X.
C INTEGER N
C
C NPTS --> The minimum number of points needed in a window for the
C particular evaluation method.
C INTEGER NPTS
C
C SMLWND <-- The smallest window width in X units containing
C NPTS points.
C
C NOTE
C
C This subroutine assumes X is sorted in ascending order.
C
C***********************************************************************
X IMPLICIT INTEGER (i-n),REAL (a-h,o-p,r-z),LOGICAL (q)
X REAL x(n),smlwnd
X INTEGER n,npts
C
X smlwnd = 0.0
C
C Loop over all Values of X
C
X DO 90,i = 1,n
X ixlo = i
X ixhi = i
X xval = x(i)
X ndiff = 0
C
C Loop until number of distinct points is equal to NPTS
C
X 10 IF (ndiff.EQ.npts) GO TO 80
C
C Execute if not on a boundary
C
X IF (.NOT. (ixlo.NE.1.AND.ixhi.NE.n)) GO TO 40
X xlo = xval - x(ixlo-1)
X xhi = x(ixhi+1) - xval
C
X IF (.NOT. (xlo.LT.xhi)) GO TO 20
X ixlo = ixlo - 1
X xdiff = xlo - x(ixlo+1)
X GO TO 30
X
X 20 ixhi = ixhi + 1
X xdiff = xhi - x(ixhi-1)
X 30 GO TO 70
C
C Execute if on either boundary
C
X 40 IF (.NOT. (ixlo.EQ.1)) GO TO 50
X ixhi = ixhi + 1
X xdiff = x(ixhi) - x(ixhi-1)
X GO TO 60
X
X 50 ixlo = ixlo - 1
X xdiff = x(ixlo) - x(ixlo+1)
X 60 CONTINUE
X 70 IF (xdiff.NE.0.000) ndiff = ndiff + 1
C
X GO TO 10
X
X 80 xmin = 2.0*max(xval-x(ixlo),x(ixhi)-xval)
C
C End of Until Loop
C
X IF (xmin.GT.smlwnd) smlwnd = xmin
X 90 CONTINUE
X RETURN
C
C End of Do Loop
C
X END
X SUBROUTINE wgt(x,nptwi,weight,method,xval,parmet,ntielo,ntiehi,
X + xwgts)
C**********************************************************************
C SUBROUTINE WGT(X, NPTWI, WEIGHT, METHOD, XVAL, PARMET,
C # NTIELO, NTIEHI, XWGTS)
C
C WeiGhT
C
C Function
C
C Returns vector of weights of data points in window.
C Methods are documented in lower level routines.
C
C
C**********************************************************************
C
C
C Arguments
C
C
C X --> Vector of X values in window. X(1) is lowest element
C in window.
C REAL X(NPTWI)
C
C NPTWI --> Number of points in current window.
C INTEGER NPTWI
C
C WEIGHT --> Weighting method.
C 1 - Boxcar
C 2 - Biquartic
C INTEGER WEIGHT
C
C METHOD --> Windowing method.
C 1 - Fixed width
C 2 - Nearest neighbor
C INTEGER METHOD
C
C XVAL --> Window center
C REAL XVAL
C
C PARMET --> Method-dependent windowing parameter.
C If METHOD = 1, PARMET is the width of the fixed
C width window.
C If METHOD = 2, PARMET is the number of nearest
C neighbors to be included in the window.
C REAL PARMET
C
C NTIELO --> Number of X elements on left boundary of window.
C INTEGER NTIELO
C
C NTIEHI --> Number of X elements on right boundary of window.
C INTEGER NTIEHI
C
C XWGTS --> Weights of X values in current window (parallels X).
C REAL XWGTS(NPTWI)
C
C**********************************************************************
X REAL x(nptwi),xval,parmet,xwgts(nptwi)
X INTEGER ntielo,ntiehi
X INTEGER nptwi,weight,method
C
C**********************************************************************
C
C
C Parameters
C
C
C**********************************************************************
X INTEGER wfix,wnear
X PARAMETER (wfix=1,wnear=2)
C
X INTEGER wbox,wbi
X PARAMETER (wbox=1,wbi=2)
C
C**********************************************************************
C
C
C Main body of program
C
C
X IF (.NOT. (method.EQ.wnear)) GO TO 40
C Nearest neighbor
X iparmt = int(parmet)
X IF ((wbi).NE. (weight)) GO TO 10
C Biquartic weighting
X CALL wtbqnr(x,nptwi,iparmt,ntielo,ntiehi,xval,xwgts)
X GO TO 30
X
X 10 IF ((wbox).NE. (weight)) GO TO 20
C
C Boxcar weighting
X CALL wtbcnr(x,nptwi,iparmt,ntielo,ntiehi,xwgts)
X GO TO 30
C
X 20 CONTINUE
C
X 30 CONTINUE
X 40 IF (.NOT. (method.EQ.wfix)) GO TO 80
C
C Fixed width
X IF ((wbox).NE. (weight)) GO TO 50
C Boxcar weighting
X CALL wtbcwd(x,nptwi,xwgts)
X GO TO 70
X
X 50 IF ((wbi).NE. (weight)) GO TO 60
C
C Biquartic weighting
X CALL wtbqwd(x,nptwi,xval,parmet,xwgts)
X GO TO 70
C
X 60 CONTINUE
C
X 70 CONTINUE
X 80 RETURN
C
X END
X SUBROUTINE windo(x,nx,xval,method,parmet,ixlo,nptwi,ntielo,ntiehi)
C*****************************************************************
C SUBROUTINE WINDO(X, NX, XVAL, METHOD, PARMET, IXLO, NPTWI,
C # NTIELO, NTIEHI)
C
C performs WINDOwing
C
C
C Function
C
C
C Returns the lowest index of the X element which is in the
C window, and the number of X vector elements which are in the
C window. When nearest neighbor windowing is to be done, number
C of ties at the left and right boundaries are also returned.
C
C
C Arguments
C
C
C X --> Vector of X values (parallels Y). Must be sorted
C in ascending order on input..
C REAL X(NX)
C
C NX --> Length of X array.
C INTEGER NX
C
C XVAL --> X value for which window is to be found.
C REAL NX
C
C METHOD --> Windowing method.
C 1 - Fixed width
C 2 - Nearest neighbor
C INTEGER METHOD
C
C PARMET --> Method-dependent windowing parameter.
C If METHOD = 1, PARMET is the width of the fixed
C width window.
C If METHOD = 2, PARMET is the number of nearest
C neighbors to be included in the window.
C REAL PARMET
C
C IXLO <-- Index of first element in x array in current window.
C INTEGER IXLO
C
C NPTWI <-- Number of points in current window. Ties (points
C lying on boundaries) are considered to be in the
C window.
C INTEGER NPTWI
C
C NTIELO <-- Number of X elements on left boundary of window.
C INTEGER NTIELO
C
C NTIEHI <-- Number of X elements on right boundary of window.
C INTEGER NTIEHI
C
C**********************************************************************
C
X REAL x(nx),xval,parmet
X INTEGER nx,method,ixlo,nptwi,ntielo,ntiehi
C
C**********************************************************************
C
C
C Parameters
C
C
C**********************************************************************
X INTEGER wfix,wnear
X PARAMETER (wfix=1,wnear=2)
C
C**********************************************************************
C
C Main body of program
C
C
C Call proper windowing routine
C
X IF ((wnear).NE. (method)) GO TO 10
C Nearest neighbor
X iparmt = int(parmet)
X CALL wnnrn(x,nx,xval,iparmt,ixlo,nptwi,ntielo,ntiehi)
X GO TO 30
X
X 10 IF ((wfix).NE. (method)) GO TO 20
C
C Fixed width
X CALL wnfwd(x,nx,xval,parmet,ixlo,nptwi)
X GO TO 30
C
X 20 CONTINUE
X 30 RETURN
C
X END
X SUBROUTINE wnfwd(x,nx,xcen,width,ixlo,nptwi)
C*****************************************************************
C SUBROUTINE WNFWD(X, NX, XCEN, WIDTH, IXLO, NPTWI)
C
C WiNdow, Fixed WiDth
C
C
C Function
C
C
C Determines which elements of the X vector lie within window of
C specified width (WIDTH) and center (XCEN).
C
C**********************************************************************
C
C
C Method
C
C The X limits of the window are calculated. The elements of the X
C vector lying closest to these limits are found by binary search.
C The number of X vector elements in the window is calculated from
C the indices of the values at the limits.
C
C*****************************************************************
C
C
C Arguments
C
C
C X --> Vector of X values, sorted in ascending order.
C REAL X(NX)
C
C NX --> Number of points X vector.
C INTEGER NX
C
C XCEN --> X value about which window is to be centered.
C REAL XCEN
C
C WIDTH --> Width of window.
C REAL WIDTH
C
C IXLO <-- Index of lowest X element in window.
C INTEGER IXLO
C
C NPTWI <-- Number of points in window.
C INTEGER NPTWI
C
C*****************************************************************
X REAL x(nx),xcen,width
X INTEGER nx,ixlo,nptwi
C
C*****************************************************************
C
C
C Variables
C
C
C IXHI --- Index of highest X element in window.
C INTEGER IXHI
C
C XLO --- Lowest X element in window.
C REAL XLO
C
C XHI --- Highest X element in window.
C REAL XHI
C
C QINCR --- .TRUE. if the elements of X are sorted in ascending
C order, otherwise .FALSE.
C LOGICAL QINCR
C
C QVALLO --- .TRUE. if the index returned by IBSRL should have the
C property that the value searched for is less than
C the value whose index is returned, otherwise .FALSE.
C LOGICAL QVALLO
C
C**********************************************************************
X REAL xlo,xhi
X INTEGER ixhi
X LOGICAL qincr,qvallo
C
C**********************************************************************
C
C
C Main body of program
C
C
C Find X values at limits of window
C
X xlo = xcen - width/2.0
X IF (xlo.LT.x(1)) xlo = x(1)
C
X xhi = xcen + width/2.0
X IF (xhi.GT.x(nx)) xhi = x(nx)
C
C Get index of lowest element in window
C
X qincr = .TRUE.
X qvallo = .TRUE.
C
X ixlo = ibsrl(xlo,x,nx,qincr,qvallo)
C
C Get index of highest element in window
C
X qvallo = .FALSE.
C
X ixhi = ibsrl(xhi,x,nx,qincr,qvallo)
C
C Set number of points in window
C
X nptwi = ixhi - ixlo + 1
C
X RETURN
X
X END
X SUBROUTINE wnnrn(x,nx,xcen,knbr,ixlo,nptwi,ntielo,ntiehi)
C*****************************************************************
C SUBROUTINE WNNRN(X,NX,XCEN,KNBR,IXLO,NPTWI,NTIELO,NTIEHI)
C
C WiNdowing - NeaRest Neighbor
C
C
C Function
C
C
C Selects specified number of nearest neighbors of a given X
C value.
C
C**********************************************************************
C
C
C Method
C
C Beginning at X(1), a binary search of the X vector is done,
C using COMPLIB routine IBSRL, to determine the element closest
C to XCEN. This index of this element is IXCENT.
C
C X is scanned left and right from IXCENT to find the KNBR values
C nearest to XCEN. If the left limit of the vector, element 1,
C is reached, scanning continues to the right only. Reaching
C the right limit is handled in like manner. The distances from
C XCEN to the right and left elements are compared. If the right
C element is closer, the array index denoting the right window
C boundary is incremented. If the left element is closer, the
C left boundary is decremented.
C
C X elements adjacent to the left and right window boundary
C elements are examined to determine if ties exist at the
C boundaries. If multiple elements of the same value are detected
C at a boundary, their number is determined, NTIELO being the
C number of tied elements on left boundary, and NTIEHI being
C the number of tied elements on the right boundary. The value
C IXLO is decremented to include all ties on left boundary in
C the window, and the value NPTWI is incremented to include all
C ties in the window.
C
C**********************************************************************
X IMPLICIT LOGICAL (q)
C**********************************************************************
C
C
C Arguments
C
C
C X --> Vector of X values, sorted in ascending order.
C REAL X(NX)
C
C NX --> Number of points X vector.
C INTEGER NX
C
C XCEN --> X value about which window is to be centered.
C REAL XCEN
C
C KNBR --> Number of nearest neighbors to be included in
C window.
C INTEGER KNBR
C
C IXLO <-- Index of lowest X element in window.
C INTEGER IXLO
C
C NPTWI <-- Number of points in window (including values on
C window boundaries).
C INTEGER NPTWI
C
C NTIELO <-- Number of X elements on left boundary of window.
C INTEGER NTIELO
C
C NTIEHI <-- Number of X elements on right boundary of window.
C INTEGER NTIEHI
C
C*****************************************************************
X REAL x(nx),xcen
X INTEGER nx,knbr,ixlo,nptwi,ntielo,ntiehi
C
C*****************************************************************
C
C
C Variables
C
C
C DL --- Distance between center of window and current
C candidate for inclusion on left side of window.
C REAL DL
C
C DR --- Distance between center of window and current
C candidate for inclusion on right side of window.
C REAL DR
C
C IXHI --- Index of highest X element in window.
C INTEGER IXHI
C
C IXCENT --- Index of X element closest to XCEN.
C INTEGER IXCENT
C
C IBLLO --- Index of first tied value on left boundary.
C INTEGER IBLLO
C
C IBLHI --- Index of last tied value on left boundary.
C INTEGER IBLHI
C
C IBRLO --- Index of first tied value on right boundary.
C INTEGER IBRLO
C
C IBRHI --- Index of last tied value on right boundary.
C INTEGER IBRHI
C
C QINCR --- .TRUE. if the elements of X are sorted in ascending
C order, otherwise .FALSE.
C LOGICAL QINCR
C
C QVALLO --- .TRUE. if the index returned by IBSRL should have the
C property that the value searched for is less than
C the value whose index is returned, otherwise .FALSE.
C LOGICAL QVALLO
C
C QNOTIE --- used in testing for tied points at window boundaries
C LOGICAL QNOTIE
C
C*****************************************************************
X REAL dl,dr
X INTEGER ixhi,ixcent,ibllo,iblhi,ibrlo,ibrhi
X LOGICAL qincr,qvallo,qnotie
C
C*****************************************************************
C
C Logical function to see if two values are approximately equal
C
X qeqrl(vx,vy) = (abs(vx-vy)/amax1(abs(vx)+abs(vy),1.0E-20)) .LT.
X + (1.0E-5)
C
C*****************************************************************
C
C
C Main body of program
C
C
C Find two X elements bracketing XCEN on either side
C X(NDX-1) XCEN X(NDX)
C
X qincr = .TRUE.
X qvallo = .TRUE.
C
X ndx = ibsrl(xcen,x,nx,qincr,qvallo)
C
C Choose element closer to XCEN
C
X IF (.NOT. (ndx.GT.1)) GO TO 30
X IF (.NOT. (abs(xcen-x(ndx)).LT.abs(xcen-x(ndx-1)))) GO TO 10
X ixcent = ndx
X GO TO 20
X
X 10 ixcent = ndx - 1
X 20 GO TO 40
X
X 30 ixcent = ndx
C
C Set initial conditions for search for KNBR nearest neighbors
C
X 40 ilo = ixcent
X ihi = ixcent + 1
X qllim = .FALSE.
X qrlim = .FALSE.
C
X IF (.NOT. (ilo.LE.1)) GO TO 50
X qllim = .TRUE.
X ilo = 1
X 50 IF (.NOT. (ihi.GE.nx)) GO TO 60
C
X qrlim = .TRUE.
X ihi = nx
X 60 IF (.NOT. ((ihi-ilo+1).LT.knbr)) GO TO 130
C
C Add point closest to center of window to window until KNBR
C points are in window. If either limit of X vector is reached,
C add points from other side of center only.
C
X IF (.NOT. (qrlim)) GO TO 70
X ilo = ilo - 1
X GO TO 120
X
X 70 IF (.NOT. (qllim)) GO TO 80
X ihi = ihi + 1
X GO TO 120
X
X 80 dl = abs(xcen-x(ilo))
X dr = abs(xcen-x(ihi))
X IF (.NOT. (dl.GT.dr)) GO TO 90
X ihi = ihi + 1
X GO TO 110
X
X 90 IF (.NOT. (dr.GT.dl)) GO TO 100
X ilo = ilo - 1
X GO TO 110
X
X 100 ihi = ihi + 1
X ilo = ilo - 1
X 110 CONTINUE
X 120 IF (ilo.EQ.1) qllim = .TRUE.
X IF (ihi.EQ.nx) qrlim = .TRUE.
X GO TO 60
X
X 130 ixlo = ilo
C
C Save indices of upper and lower bounds of window
C
X ixhi = ihi
C
C If we haven't hit the left side of the data set then
C check for tied values on left outside of window
C
X IF (.NOT. (.NOT.qllim)) GO TO 180
X qnotie = .FALSE.
X ilo = ixlo
X ibllo = ixlo
C
X 140 IF (qllim .OR. qnotie) GO TO 170
X ilo = ilo - 1
X IF (.NOT. (qeqrl(x(ilo),x(ixlo)))) GO TO 150
X ibllo = ilo
X IF (ilo.EQ.1) qllim = .TRUE.
X GO TO 160
X
X 150 qnotie = .TRUE.
X 160 GO TO 140
X
X 170 ixlo = ibllo
C
C Set new lower bound for window to first tied value
C
C
X GO TO 190
X
X 180 ibllo = ixlo
C
C Check for tied values on left inside of window
C
X 190 IF (.NOT. (x(ixcent).EQ.x(ixlo))) GO TO 200
X iblhi = ixcent
X GO TO 230
X
X 200 ihi = ixlo + 1
X iblhi = ixlo
C
X 210 IF (.NOT. (x(ihi).EQ.x(ixlo))) GO TO 220
X iblhi = ihi
X ihi = ihi + 1
X GO TO 210
X
X 220 CONTINUE
X 230 IF (.NOT. (.NOT.qrlim)) GO TO 280
C
C If we have'nt hit the right side of the data set then
C check for tied values on right outside window
C
X qnotie = .FALSE.
X ihi = ixhi
X ibrhi = ixhi
C
X 240 IF (qrlim .OR. qnotie) GO TO 270
X ihi = ihi + 1
X IF (.NOT. (qeqrl(x(ihi),x(ixhi)))) GO TO 250
X ibrhi = ihi
X IF (ihi.EQ.nx) qrlim = .TRUE.
X ihi = ihi + 1
X GO TO 260
X
X 250 qnotie = .TRUE.
X 260 GO TO 240
X
X 270 GO TO 290
X
X 280 ibrhi = ixhi
C
C Check for tied values on right inside window
C
X 290 IF (.NOT. (x(ixcent).EQ.x(ixhi))) GO TO 300
X ibrlo = ixcent
X GO TO 330
X
X 300 ilo = ixhi - 1
X ibrlo = ixhi
C
X 310 IF (.NOT. (x(ilo).EQ.x(ixhi))) GO TO 320
X ibrlo = ilo
X ilo = ilo - 1
X GO TO 310
X
X 320 CONTINUE
X 330 ntielo = iblhi - ibllo + 1
C
C Calculate number of tied values on each side of window
C
X ntiehi = ibrhi - ibrlo + 1
C
C Calculate number of elements in window, including tied values
C
X nptwi = ibrhi - ibllo + 1
C
X RETURN
X
X END
X SUBROUTINE wtbcnr(x,nptwi,knbr,ntielo,ntiehi,xwgts)
C*****************************************************************
C SUBROUTINE WTBCNR(X, NPTWI, KNBR, NTIELO, NTIEHI, XWGTS)
C
C WeighT BoxCar nearest NeighboR
C
C
C Function
C
C
C Calculates equal weights for all elements in window except
C those on boundary (boxcar weighting).
C
C**********************************************************************
C
C
C Method
C
C The weight of each element in the window is set to 1/(KNBR).
C
C The weights for elements tied in value with a boundary element
C are reduced by the factor
C
C (KNBR - (NPTWI-NTIELO-NTIEHI))/(NTIELO+NTIEHI).
C
C The weights are summed, and each individual weight is divided by
C that sum to normalize the sum of weights to 1.
C
C*****************************************************************
C
C
C Arguments
C
C
C X --> Vector of X values in window.
C REAL X(NPTWI)
C
C NPTWI --> Number of points in current window.
C INTEGER NPTWI
C
C KNBR --> Number of data nearest neighbors to be included
C in the window.
C INTEGER KNBR
C
C NTIELO --> Number of X elements on left boundary of window.
C INTEGER NTIELO
C
C NTIEHI --> Number of X elements on right boundary of window.
C INTEGER NTIEHI
C
C XWGTS <-- Weights of X values in current window (parallels X).
C REAL XWGTS(NPTWI)
C
C**********************************************************************
X REAL x(nptwi),xwgts(nptwi)
X INTEGER nptwi,knbr,ntielo,ntiehi
C
C**********************************************************************
C
C
C Variables
C
C
C WEIGHT --- Weight assigned to data points in window.
C REAL WEIGHT
C
C BFACT --- Factor used in accounting for points tied on
C boundary of window.
C REAL BFACT
C
C WGTSUM --- Sum of weights of points.
C REAL WGTSUM
C
C**********************************************************************
X REAL weight,bfact,wgtsum
C
C Main body of program
C
C
C Calculate weight, and fudge factor for boundary points
C
X weight = 1.0/knbr
X bfact = (knbr-float(nptwi-ntielo-ntiehi))/float(ntielo+ntiehi)
C
C Initialize sum of weights
C
X wgtsum = 0.0
C
C Accumulate weights of points on left boundary
C
X i = 1
X GO TO 20
X
X 10 i = i + 1
X 20 IF (i.GT.ntielo) GO TO 30
X xwgts(i) = weight*bfact
X wgtsum = wgtsum + xwgts(i)
X GO TO 10
X
X 30 i = ntielo + 1
X GO TO 50
X
X 40 i = i + 1
X 50 IF (i.GT.nptwi-ntiehi) GO TO 60
C
C Accumulate weights of points not on a boundary
C
X xwgts(i) = weight
X wgtsum = wgtsum + xwgts(i)
X GO TO 40
X
X 60 i = nptwi - ntiehi + 1
X GO TO 80
X
X 70 i = i + 1
X 80 IF (i.GT.nptwi) GO TO 90
C
C Accumulate weights of points on right boundary
C
X xwgts(i) = weight*bfact
X wgtsum = wgtsum + xwgts(i)
X GO TO 70
X
X 90 i = 1
X GO TO 110
X
X 100 i = i + 1
X 110 IF (i.GT.nptwi) GO TO 120
C
C Normalize weights so sum of weights is 1
C
X xwgts(i) = xwgts(i)/wgtsum
X GO TO 100
X
X 120 RETURN
C
X END
X SUBROUTINE wtbcwd(x,nptwi,xwgts)
C*****************************************************************
C SUBROUTINE WTBCWD(X, NPTWI, XWGTS)
C
C WeighT BoxCar fixWiDth
C
C
C Function
C
C
C Boxcar weighting. Calculates weights of points in fixed
C width window by assigning all points a weight of 1/NPTWI.
C
C**********************************************************************
C
C
C Method
C
C The weight of each element in the window is set to 1/(NPTWI).
C
C*****************************************************************
C
C
C Arguments
C
C
C X --> Vector of X values in window.
C REAL X(NPTWI)
C
C NPTWI --> Number of points in current window.
C REAL X(NPTWI)
C
C XWGTS <-- Weights of X values in current window (parallels X).
C REAL XWGTS(NPTWI)
C
C*****************************************************************
X REAL x(nptwi),xwgts(nptwi)
X INTEGER nptwi
C**********************************************************************
C
C
C Variables
C
C
C WEIGHT --- Weight of each point
C REAL WEIGHT
C
C**********************************************************************
X REAL weight
C
C**********************************************************************
C
C
C Main body of program
C
C
X weight = 1.0/float(nptwi)
C
X i = 1
X GO TO 20
X
X 10 i = i + 1
X 20 IF (i.GT.nptwi) GO TO 30
X xwgts(i) = weight
X GO TO 10
C
X 30 RETURN
X
X END
X SUBROUTINE wtbqnr(x,nptwi,knbr,ntielo,ntiehi,xcen,xwgts)
C*****************************************************************
C SUBROUTINE WTBQNR(X, NPTWI, KNBR, NTIELO, NTIEHI, XCEN, XWGTS)
C
C WeighT BiQuartic nearest NeighboR
C
C
C Function
C
C
C Calculates weights for all elements of X vector in window
C including those on boundary. Weights based on biquartic
C function of distance from window center.
C
C**********************************************************************
C
C
C Method
C
C The weights for all elements in the window are calculated to be
C
C (1 - ((X(I)-XCEN)/WIDTH)**2)**2.
C
C The weights for elements tied in value with a boundary element
C are reduced by the factor
X
C (KNBR - (NPTWI-NTIELO-NTIEHI))/(NTIELO+NTIEHI).
C
C The weights are summed, and each individual weight is divided by
C that sum to normalize the sum of weights to 1.
C
C*****************************************************************
C
C
C Arguments
C
C X --> Vector of X values in window.
C REAL X(NPTWI)
C
C NPTWI --> Number of points in current window.
C INTEGER NPTWI
C
C KNBR --> Number of data nearest neighbors to be included
C in the window.
C INTEGER KNBR
C
C NTIELO --> Number of X elements on left boundary of window.
C INTEGER NTIELO
C
C NTIEHI --> Number of X elements on right boundary of window.
C INTEGER NTIEHI
C
C XCEN --> X value of window center
C REAL XCEN
C
C XWGTS <-- Weights of X values in current window (parallels X).
C REAL XWGTS(NPTWI)
C
C*****************************************************************
X REAL x(nptwi),xcen,xwgts(nptwi)
X INTEGER nptwi,knbr,ntielo,ntiehi
C
C**********************************************************************
C
C Variables
C
C
C WIDTH --- Range of X values in window.
C REAL WIDTH
C
C BFACT --- Factor used in accounting for points tied on
C boundary of window.
C REAL BFACT
C
C WEIGHT --- Weight assigned to data points in window.
C REAL WEIGHT
C
C WGTSUM --- Sum of weights of points.
C REAL WGTSUM
C
C**********************************************************************
X REAL width,bfact,weight,wgtsum
C
C Main body of program
C
C
C Calculate fudge factor for boundary points and window width
C
X IF (.NOT. ((ntielo+ntiehi).GT.2)) GO TO 10
X bfact = float(knbr- (nptwi-ntielo-ntiehi))/float(ntielo+ntiehi)
X GO TO 20
X
X 10 bfact = 1.0
C
X 20 width = abs(x(nptwi)-x(1))
C
C Initialize sum of weights
C
X wgtsum = 0.0
C
C Accumulate weights of points on left boundary
C
X i = 1
X GO TO 40
X
X 30 i = i + 1
X 40 IF (i.GT.ntielo) GO TO 50
X xwgts(i) = bfact* (1.0- ((x(i)-xcen)/width)**2)**2
X wgtsum = wgtsum + xwgts(i)
X GO TO 30
X
X 50 i = ntielo + 1
X GO TO 70
X
X 60 i = i + 1
X 70 IF (i.GT.nptwi-ntiehi) GO TO 80
C
C Accumulate weights of points not on a boundary
C
X xwgts(i) = (1.0- ((x(i)-xcen)/width)**2)**2
X wgtsum = wgtsum + xwgts(i)
X GO TO 60
X
X 80 i = nptwi - ntiehi + 1
X GO TO 100
X
X 90 i = i + 1
X 100 IF (i.GT.nptwi) GO TO 110
C
C Accumulate weights of points on right boundary
C
X xwgts(i) = bfact* (1.0- ((x(i)-xcen)/width)**2)**2
X wgtsum = wgtsum + xwgts(i)
X GO TO 90
X
X 110 i = 1
X GO TO 130
X
X 120 i = i + 1
X 130 IF (i.GT.nptwi) GO TO 140
C
C Normalize weights so sum of weights is 1
X xwgts(i) = xwgts(i)/wgtsum
X GO TO 120
X
X 140 RETURN
C
X END
X SUBROUTINE wtbqwd(x,nptwi,xcen,width,xwgts)
C*****************************************************************
C SUBROUTINE WTBQWD(X, NPTWI, XCEN, WIDTH, XWGTS)
C
C WeighT BiQuartic fixWiDth
C
C
C Function
C
C
C Biquartic weighting. The distance from XCEN for each point is
C calculated as a function of the window width. The weights are
C normalized to sum to 1.
C
C**********************************************************************
C
C
C Method
C
C The distance from XCEN is calculated as a fraction of the
C window width, squared and added to one. This quantity is then
C squared again. The weights are normalized to sum to one by
C dividing each by the sum of them all.
C
C*****************************************************************
C
C
C Arguments
C
C X --> Vector of X values in window.
C REAL X(NPTWI)
C
C NPTWI --> Number of points in current window.
C INTEGER NPTWI
C
C XCEN --> X value of window center
C REAL XCEN
C
C WIDTH --- Width of window.
C REAL WIDTH
C
C XWGTS <-- Weights of X values in current window (parallels X).
C REAL XWGTS(NPTWI)
C
C*****************************************************************
X REAL x(nptwi),xcen,width,xwgts(nptwi)
X INTEGER nptwi
C
C*****************************************************************
C
C Variables
C
C
C WGTSUM --- Sum of weights of points.
C REAL WGTSUM
C
C**********************************************************************
X REAL wgtsum
C
C**********************************************************************
C
C
C Main body of program
C
C
X wgtsum = 0.0
C
C Calculate weight of each point
C
X i = 1
X GO TO 20
X
X 10 i = i + 1
X 20 IF (i.GT.nptwi) GO TO 30
X xwgts(i) = (1.0- ((x(i)-xcen)/width)**2)**2
X wgtsum = wgtsum + xwgts(i)
X GO TO 10
X
X 30 i = 1
X GO TO 50
X
X 40 i = i + 1
X 50 IF (i.GT.nptwi) GO TO 60
X xwgts(i) = xwgts(i)/wgtsum
X GO TO 40
C
C Normalize weights to sum to 1
C
C
X 60 RETURN
X
X END
END_OF_FILE
if test 95080 -ne `wc -c <'windowing.f'`; then
echo shar: \"'windowing.f'\" unpacked with wrong size!
fi
# end of 'windowing.f'
fi
if test -f 'windowing.s' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'windowing.s'\"
else
echo shar: Extracting \"'windowing.s'\" \(6691 characters\)
sed "s/^X//" >'windowing.s' <<'END_OF_FILE'
windowing <- function(x,y,xdes=defxdes,method="fixwidth",smparam=defpar,
X wtmethod="boxcar",estmethod="mean",
X dermethod="linear",quantile=0.5)
X{
X########################################################################
X#
X# WINDOWING
X#
X# FUNCTION
X#
X# Calls the fortran windowing routines to perform smoothing
X#
X# ARGUMENTS
X#
X# x --> Independant variable for smoothing
X#
X# y --> Dependant variable for smoothing
X#
X# xdes --> Values at which smoothed y values are desired, i.e. window
X# center values. Default is 20 equally spaced values over the
X# range of x
X#
X# method --> Character string specifing the windowing method
X#
X# smparam --> Windowing parameter. If method is "fixwidth" it's the width
X# of the window, the default is 1/10 the range of x. If method
X# is "nearest neighbor" it's the number of neighbors to be
X# included in the window, the default is max(10,len(x)/10)
X#
X# wtmethod --> Character string specifing the weighting method
X#
X# estmethod --> Character string specifing the estimation method
X#
X# dermethod --> Character string specifing the method to calculate
X# the derivative
X#
X# quantile --> Quantile to be used in the quantile estimation method
X#
X#########################################################################
X#
X# make sure x and y are there, are numeric, and the same length
X#
X if(missing(x)) stop("Argument x is missing, with no default")
X if(missing(y)) stop("Argument y is missing, with no default")
X if(!is.numeric(x)) stop("x must be numeric")
X if(!is.numeric(y)) stop("y must be numeric")
X if(len(x) != len(y)) stop("Lengths of x and y must match")
X#
X# remove any NAs
X#
X xx <- x[!(is.na(x)) & !(is.na(y))]
X yy <- y[!(is.na(x)) & !(is.na(y))]
X nxy <- len(xx)
X if(nxy < 2) stop("x and y yield less than 2 points")
X#
X# sort them if we need to
X#
X if (any(c(-1.0E-35,xx) > c(xx,1.0E35))) {
X permut <- sort.list(xx)
X xx <- xx[permut]
X yy <- yy[permut]
X }
X#
X# get the various methods
X#
X if(!missing(method)) {
X if(!is.character(method))
X stop("method must be of mode character")
X target <- c("fixwidth","nearest neighbor")
X imethod <- charmatch(method,target)
X if (is.na(imethod)) stop("windowing method unknown")
X }
X else imethod <- 1
X#
X if(!missing(wtmethod)) {
X if(!is.character(wtmethod))
X stop("wtmethod must be of mode character")
X target <- c("boxcar","biquadratic")
X iwtmethod <- charmatch(wtmethod,target)
X if (is.na(iwtmethod)) stop("weighting method unknown")
X else if (!(iwtmethod > 0)) stop("weighting method ambiguous")
X }
X else iwtmethod <- 1
X#
X if(!missing(estmethod)) {
X if(!is.character(estmethod))
X stop("estmethod must be of mode character")
X target <- c("mean","linear fit","quadratic fit","maximum",
X "minimum","quantile","standard deviation",
X "derivatives","window width","number of points")
X iestmethod <- charmatch(estmethod,target)
X if (is.na(iestmethod)) stop("estimation method unknown")
X else if (!(iestmethod > 0)) stop("estimation method ambiguous")
X }
X else iestmethod <- 1
X#
X if(!missing(dermethod)) {
X if(!is.character(dermethod))
X stop("dermethod must be of mode character")
X target <- c("linear fit","quadratic fit")
X idermethod <- charmatch(dermethod,target)
X if (is.na(idermethod)) stop("Derivative method unknown")
X }
X else idermethod <- 1
X#
X# initialize the variables that need it
X#
X xlo <- min(xx)
X xhi <- max(xx)
X#
X# if xdes wasn't given, define the default value for it
X#
X if(missing(xdes)) {
X xd <- c(0:19)
X delx <- (xhi-xlo)/19
X defxdes <- (xd*delx)+xlo
X }
X#
X# else since xdes was given make sure it is ok
X#
X else {
X if(!is.numeric(xdes)) stop("xdes must be numeric")
X xdes <- xdes[!is.na(xdes)]
X if(len(xdes) < 1) stop("xdes is all NA")
X }
X ndes <- len(xdes)
X#
X# if smparam wasn't given, define the default value for it
X#
X if(missing(smparam)) {
X if(imethod==1) defpar <- max(1e-10,(xhi-xlo)/10)
X else defpar <- max(2,min(10,len(x)/10))
X }
X#
X# else since smparam was given make sure it's ok
X#
X else {
X if(len(smparam)>1) {
X smparam <- smparam[1]
X warning("Only the first element of smparam used")
X }
X if(is.na(smparam)) stop("Missing value not allowed")
X if(imethod==1) {
X if(smparam < 1e-10) stop("window width is too small")
X }
X else {
X if(smparam < 2) stop("number of neighbors is too few")
X }
X }
X#
X# if quantile was given make sure it is ok
X#
X if((iestmethod == 7) & (!missing(quantile))) {
X if(len(quantile) > 1) {
X quantile <- quantile[1]
X warning("Only the first element of quantile used")
X }
X if(is.na(quantile)) stop("Missing value not allowed")
X if((quantile >= 1.0) & (quantile <= 0.0))
X stop("quantile must be between 0.0 and 1.0")
X }
X#
X sortar <- matrix(2,nxy)
X storage.mode(sortar) <- "single"
X xwgts <- vector(mode="numeric",length=nxy)
X storage.mode(xwgts) <- "single"
X xmiss <- rep(0,nxy)
X qcross <- FALSE
X#
X# call the fortran to do the windowing
X#
X wind <- .Fortran("window",
X as.single(xx),
X as.single(yy),
X as.integer(nxy),
X as.single(xdes),
X ysmo = single(length = ndes),
X as.integer(ndes),
X as.integer(imethod),
X as.single(smparam),
X as.integer(idermethod),
X as.integer(iwtmethod),
X as.integer(iestmethod),
X as.single(quantile),
X sortar,
X xwgts,
X xmiss=as.integer(xmiss),
X as.logical(qcross),
X integer(length = 1))
X ysmo <- wind$ysmo
X if(sum(wind$xmiss) != 0) {
X xdes[wind$xmiss] <- NA
X ysmo[wind$xmiss] <- NA
X }
X return(cbind(xdes,ysmo))
X}
END_OF_FILE
if test 6691 -ne `wc -c <'windowing.s'`; then
echo shar: \"'windowing.s'\" unpacked with wrong size!
fi
# end of 'windowing.s'
fi
if test -f 'windows.doc' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'windows.doc'\"
else
echo shar: Extracting \"'windows.doc'\" \(32562 characters\)
sed "s/^X//" >'windows.doc' <<'END_OF_FILE'
X
X
X
X
X
X
X
X
X
X
X WINDOWS
X Version 2.0 August 1990
X
X
X
X
X
X Windowed / Kernel / Moving Average Smoothing and Estimation
X
X
X
X
X
X
X
X Barry W. Brown
X David Gutierrez
X James Lovato
X Marty Spears
X Patrick Tibbits
X
X
X
X
X
X This work was supported by grant CA-16672 from the National Cancer
X Institute.
X
X
X
X
X
X Copyright 1988 for
X The University of Texas, M.D. Anderson Cancer Center
X Department of Biomathematics, Box 237
X 1515 Holcome Blvd.
X Houston, TX 77030
X
X
X Contact: BWB at above address or bwb@mdaali.cancer.utexas.edu
X
X
X
X WINDOWS can be freely copied and (noncommercially) distributed.
X
X
X
X
X
X
X
X
X
X
X Part I. Introduction to
X WINDOW and KERNEL Estimation
X
X
X A. Overview
X
X
X The primary use of window and kernel estimation techniques is to
X display the form of the mean of a dependent variable, Y, as a function
X of an independent variable, X, when Y is subject to noise. In
X addition to the mean of Y, other properties (e.g., the variance) can
X be estimated as a function of X. Windowing methods can also be used
X to determine whether parameters of equations fitting Y to X vary with
X X.
X
X A window is a set of X values, for example, those values within a
X particular distance from the specified center of the window. Window
X estimation performs a summary operation, e.g., calculation of the mean
X of Y, on data within a series of windows, the centers of which span
X the values of X. The statistic generated by the summary is associated
X with the value of X at the center of each window.
X
X In performing the summary operation, the set of points in each
X window may be weighted. If weights are used, the points near the
X center given greater weights than points at the periphery. The
X function of the distance from the center of the window that determines
X the weights given points is termed the kernel of the estimator.
X
X Properties of Y other than the mean can also be estimated. Those
X provided by the WINDOWS program include minimum and maximum values,
X percentiles, standard deviations and the derivative of Y with respect
X to X.
X
X Windowing techniques can be used to examine the appropriateness of
X a model that fits Y as a function of X. The techniques are used to
X determine whether the parameters of the model change systematically
X with X. If so then the model is not a good summary of the data. The
X only fitting capability in the current version of WINDOWS is the
X estimation of the derivative of Y with respect to X. The derivative
X provides a test of the adequacy of a linear model. Additional
X enhancements in examining fits of models are planned.
X
X
X B. Example
X
X Figure 1 shows the results of a simulated study of the effect of
X the dose of some agent on the response rates of subjects. Each of 101
X subjects received a dose of the agent; the doses given are 0, 1, ...,
X 100 units. The outcome is 0 if the subject does not respond to the
X dose received, and 1, if the subject does respond. The results
X graphed in Figure 1 show that more subjects respond at high doses than
X at low doses.
X
X
X
X
X Figure 1
X
X ............................................................
X .
X 1.0.. o o o o o o o oooo ooooooooooooooooooooooo
X .
X .
X .
X .
X .
X .
X .
X .
X 0.8..
X .
X .
X .
X .
X .
X .
X .
X .
X 0.6..
X R .
X e .
X s .
X p .
X o .
X n .
X s .
X e .
X .
X 0.4..
X .
X .
X .
X .
X .
X .
X .
X .
X 0.2..
X .
X .
X .
X .
X .
X .
X .
X .
X 0.0.. ooooooooooooooooooooooooooooooooooooooooo oooo o o
X .
X ............................................................
X 0 20 40 60 80 100
X
X Dose
X
X
X
X
X The dose administered is the independent variable (X) as it is
X controlled by the experimenter; the outcome variable, response, (Y) is
X the dependent variable because it varies systematically with X. In
X many cases, the difference between independent and dependent variable
X is difficult to distinguish, and this difference may only reflect the
X viewpoint of the observer. The X and Y terminology results from the
X usual graphic representation of the independent variable on the x-axis
X and the dependent variable on the y-axis.
X
X The purpose of the experiment was to obtain information on the
X probability of response as a function of dose. Figure 1 is not
X particularly useful for this purpose. It would be desirable to use
X the data to obtain an estimate of the probability of response as a
X function of the dose. Figures 2, 3, and 4 show such estimates using
X different window widths. The estimates are shown a s's with periods
X connecting them. Since the data was generated for a known function
X relating response rate to dose, the true relation is shown as
X asterisks.
X
X To obtain the figures, an estimate of the probability of response
X was performed at doses of 0, 2, 4, ..., 96, 98, 100. To estimate the
X response rate at each dose, X, all of the data obtained at doses
X within a window centered at X is examined. The proportion of
X responders in this data is used as the estimate of the response rate
X at X.
X
X The window centered at X contains all data with doses within half
X the window width of the dose X, since the window extends from its
X center to half its width on the right and half of its width on the
X left. As X is moved, the window moves with it.
X
X
X
X
X Figure 2
X Window Width = 5
X ............................................................
X .
X 1.0.. ssss
X . .
X . . *
X . . ***
X . .**
X . .s ss*
X . .. *
X . ..**
X . .**
X 0.8.. s**s
X . *
X . *
X . .*
X . *
X . *s
X . *.
X . * .
X . * .
X 0.6.. s .
X R . ** .
X e . *. .
X s . **. .
X p . * ...
X o . s. **. ss
X n . ..** .
X s . . .* .
X e . . ** .
X . . * .
X 0.4.. s**ssss
X . .*
X . .*
X . *
X . *s
X . *
X . ** .
X . ** .
X . ** .
X 0.2.. s ssss s
X . s . ** . .
X . . .** . .
X . . *** ..
X . s sss sss ss sssss ss
X . . . *****
X . . ********* . .
X . ******** . . ..
X . . . .. .
X 0.0.. sss s ss s
X .
X ............................................................
X 0 20 40 60 80 100
X
X Dose
X
X
X
X Figure 3
X Window Width = 25
X ............................................................
X .
X 1.0..
X .
X . s
X . .ss*
X . s*
X . s*
X . *
X . **
X . *ss
X 0.8.. **
X . *s
X . *.
X . *.
X . * .
X . * s
X . * .
X . * s
X . *s.
X 0.6.. *
X R . **.
X e . * s
X s . **s.
X p . *
X o . **s
X n . ** .
X s . * .
X e . ** .
X . * s
X 0.4.. ** .
X . * .s
X . *s s
X . * .
X . **s.
X . ss.
X . *s.
X . *s
X . **
X 0.2.. *ss
X . **
X . **s.
X . sss s
X . ssss*
X . s *ss**
X . ss s*ssss*s*
X . s**ss*s*s
X .
X 0.0..
X .
X ............................................................
X 0 20 40 60 80 100
X
X Dose
X
X
X Figure 4
X Window Width = 50
X ............................................................
X .
X 1.0..
X .
X . *
X . ***
X . **
X . **
X . *
X . **
X . ** s
X 0.8.. **
X . * s
X . * .
X . * ss
X . * ..
X . * ss
X . * .s
X . * .s
X . * s.
X 0.6.. * s.
X R . ** s.
X e . *ss
X s . **
X p . *s
X o . *s.
X n . **
X s . *s
X e . **
X . *s
X 0.4.. **
X . ss
X . *
X . s.
X . s*
X . .*
X . *s
X . .ss
X . .s*
X 0.2.. ss*
X . s.**
X . s**
X . .ss**
X . ssss**
X . sss*s***
X . s sss ssssss****
X . *s******
X .
X 0.0..
X .
X ............................................................
X 0 20 40 60 80 100
X
X Dose
X
X
X
X C. Choosing a Window Size
X
X As the window size increases, the amount of noise in the estimate
X decreases. Unfortunately, as the window size gets quite large, the
X estimate may be so smooth that significant features in the data are
X missed. Extreme cases illustrate this. If the window width is so
X small that there is only one point in each window, then there is no
X smoothing of the raw data at all. If the window width is so large
X that all points are included in every window, then the estimated curve
X is a horizontal line regardless of the data.
X
X In Figure 2, the window width is 5 and the smoothed curve jumps
X around a great deal and frequently is not close to the true value.
X When the window width is increased to 25 (Figure 3), the fit of the
X smoothed curve is much better. Increasing the window width to 50
X (Figure 4) smooths the curve still more, but the smoothed curve
X deviates from the true value appreciably at the right end.
X
X A rule of thumb for choosing the window width is to decrease the
X window size until there is appreciable noise in the estimate, and then
X increase it a bit. Fortunately, the estimate usually changes very
X slowly with window size making this heuristic guideline feasible.
X
X Theoretical work is ongoing on methods to choose window widths to
X minimize some measure of the total error in window estimation; this
X error has components of noise and oversmoothing. The method for
X choosing window widths incorporated into the WINDOWS program is cross
X validation. For this method, a list of possible window widths is
X used. For each width, an estimate of Y is made separately using a
X window centered at each data point. The estimate of Y, however, does
X not use the data point at which it is centered. The sum of squares of
X deviations of the estimated Y values from the actual Y values is the
X cross validation measure. The smaller this value, the better the
X choice of window width.
X
X Cross validation is only available for methods that estimate Y as
X a function of X, not for the estimates of other properties of Y.
X
X
X D. Known Difficulties with Windowing
X
X Windowing becomes less accurate at the edges of the data because
X there are fewer points in the window as the center moves past the
X range of the data. The problem is worse in cases in which the data is
X on only one side of the window center. Windowing also tends to smooth
X out rapid changes in Y with X. Both problems can be ameliorated
X somewhat by the use of a smoothing function more robust than the
X weighted mean of our example. WINDOWS offers linear and quadratic
X regression as alternatives to the weighted mean to provide estimates
X of the value of Y as a function of X. These methods fit either a
X linear or quadratic equation in Y as a function of X within each
X window and return the Y value at the center. Either regression should
X perform better at the extremes of the data than does the weighted
X mean; this is particularly the case when Y is monotone with X in these
X regions. The quadratic fit should track rapid changes in Y as a
X function of X better than the weighted mean or linear regression.
X
X
X
X
X
X
X E. Comparison with Stratification and Parametric Estimation
X
X Window estimation should be contrasted with stratification in
X which the independent variable is categorized into a number of
X mutually exclusive groups so that each point belongs to only one
X group. The problem with stratification is that the number of groups
X must typically be small or else there will be few points in each
X group. For example, if we have decided that there should be 25 points
X in each group corresponding to the window size of 25, then there could
X be only four dose categories. The resulting estimate would look like
X that obtained from the window method on every twelfth (50/4) point.
X The graph of response rate versus dose would contain only four points,
X and increasing the number of points would greatly reduce the size in
X each group and thus increase the noise associated with each point.
X
X In windowing methods using each data point is typically used in
X each of several windows. This allows the number of points displayed
X on the curve to be independent of the amount of noise at each point,
X but there is a corresponding cost. Statistical techniques for testing
X for changes in categorized data rely heavily on the fact that each
X data points occur in only one category. There are no comparable
X techniques for windowed estimates. Such estimates should be
X considered exploratory, designed to show the form of the data, rather
X than confirmatory, intended to demonstrate differences.
X
X An alternative to windowing is the use of parametric methods in
X which a form for the relation of response to dose is assumed. The
X most commonly used form for the shape of a dose response curve is that
X provided by the logistic model. Figure 5 shows the windowed estimate
X with a window width of 25 on a transformed scale. If the logistic
X model is correct, the curve should be a straight line. A quite
X noticeable curvature is present indicating that an unmodified logistic
X model is probably not a good fit to the data. Since the data was
X generated from a model that is quadratic, not linear, on the logistic
X scale, this finding is not a surprise.
X
X
X
X Figure 5
X
X ......................................................................
X .
X 3..
X . .
X . .
X . .
X .
X .
X 2.. .
X .
X . .
X . .
X . .
X .
X 1..
X . .
X R . .
X e . .
X s .
X p . . .
X o 0.. .
X n .
X s . .
X e . .
X . ..
X . .
X -1.. ...
X . .
X .
X . ..
X . .
X . . .. .
X -2.. . .
X . . .
X .
X . . ..
X . .
X . . . ... . .
X -3.. . .
X . .
X ......................................................................
X 0 20 40 60 80 100
X
X Dose
X
X
X
X
X This examination illustrates our usual use of windowing techniques
X as a method of exploratory data analysis and as an aid to model
X identification. It also demonstrates why WINDOWS is usually used in
X conjunction with other computer packages. WINDOWS provides only crude
X character graphics, it does not provide transformation capabilities,
X nor does it perform the fitting of parametric models. Data
X transformation, model fitting, and quality graphics must be found
X elsewhere. The interface to other packages that WINDOWS provides is
X the ability to write smoothed data to a file. There is a version of
X WINDOWS that runs within the S statistical package (copyright AT&T),
X and it is freely available from statlib.
X
X E. Windows, Kernels, and Moving Averages
X
X These three techniques are minor variants on the same method. The
X basic window technique has been explained. A kernel or weighting
X function is sometimes used to weight points according to their
X distance from the center of the window. This variation on the
X estimation technique is termed kernel estimation. Generally, points
X near the center of the window are given more weight than points far
X from the center in computing the estimate. This form of weighting
X results in a greater degree of smoothing than does a simple window
X estimate using a fixed window width because the weight associated with
X a point smoothly decreases as the window center moves away from the
X point. In the simple windowing technique, the weight is one if the
X point is in the window and the weight drops abruptly to zero when the
X window moves away from the point.
X
X The ordinary window estimate is a special case of kernel
X estimation in which the weighting function is identically one inside
X the window. This function is referred to as the boxcar weighting
X function because its graph resembles a boxcar (the function is zero to
X the left of the window, one inside the window, and zero outside the
X window). In the current version of the WINDOW package, only the
X biquartic weight is provided in addition to the boxcar. The graph of
X the biquartic kernel looks like a bell shaped curve that decreases to
X zero at the window boundary.
X
X Moving average techniques use an average of the data within a
X fixed number of points of the point at which the estimate will be
X applied. Frequently, the points are weighted according to their
X distance from the point of interest.
X
X Window and kernel techniques are quite similar to the moving
X average method. Indeed, the moving average and windowing technique
X are identical if the moving average includes points on both sides of
X the position at which the estimate is desired. Moving average
X techniques are frequently one-sided, particularly in time series work
X in which only past values are included. The WINDOW routines include
X only two-sided estimates. As will be noted below, estimates other
X than averages are available with windowing or kernel techniques and so
X do not technically fit into the framework of moving averages, although
X they are there in spirit.
X
X
X
X
X
X
X F. Fixed Width and Nearest Neighbor Estimation
X
X If the values of the independent variable are close to equally
X spaced, then fixed width window or kernel estimation works quite well.
X Many studies have an observational component that does not determines
X the spacing of x-values. For example, consider the relation of the
X distance of stars from the earth to their average size. If fixed
X window width estimation is used on data of widely varying spacing,
X some windows will contain a large number of observations while others
X will contain few. Consequently, some estimates will be well founded
X while others are based on very little data.
X
X In this circumstance, it may be wise to adjust the size of each
X window so that all windows contain the same number of observations.
X The window width is increased or decreased until it is just sufficient
X to contain a specified number of observations. This form of
X estimation is termed nearest neighbor estimation for obvious reasons.
X When using nearest neighbor estimates, the number of points to appear
X in each window is specified; when using fixed width windows, the
X window width is specified.
X
X A technical difficulty of nearest neighbor windowing occurs when
X data can be tied on the value of the independent variable. A clump of
X tied data can occur right at a point which would otherwise be a window
X endpoint. Without the clump, there are not enough points in the
X window; with it, there are too many. The solution to this dilemma in
X the WINDOW program is to include the clump, but to weight the points
X in the clump so that the total weight is precisely as needed to round
X out the number of points in the window. This weighting is applied
X before the weighting from a specified kernel function; that is the
X final weight assigned to a point depends both on the kernel and on
X whether the point belongs to a clump of points with tied values at the
X boundary of the window.
X
X G. Estimators and Fits
X More Uses for Window and Kernel Techniques
X
X Windowing and kernel techniques can perform either estimation or
X fitting of data. Estimation requires only the values of the dependent
X or y-values of the data in the window plus the weight to be applied to
X each observation. There are a large number of possible estimators.
X Here are a few examples:
X
X (1) Mean of the y-values. This is the estimate most commonly used
X and the one used in the example above.
X
X (2) Standard deviation of the y-values.
X
X (3) Number of points in the fixed width window. This gives a rough
X estimate of the density of the y-values.
X
X (4) Quantile of the y-values. A popular choice for the quantile
X is the median.
X
X (5) Kaplan-Meier estimates of the survival distribution of the
X y-values in the window. The y-values are interpreted as survival
X times. (An indicator of whether each time is a time to death or a
X time to last contact is needed in addition to the times themselves.)
X
X Fits differ from estimators as defined above in that fits require
X values of both the independent and dependent variables (x- and
X y-values) as well as the weight to be applied to each point. At least
X two useful results can be obtained from fitting simple functions to
X the data. An example of a simple function is a low order polynomial,
X say a quadratic or cubic.
X
X (1) The value of the polynomial at the center of the window can be
X used as the central value instead of the mean. This generally
X increases the degree of smoothing.
X
X (2) The derivative of the polynomial with respect to the X-value
X can be used as an estimate of the derivative of the estimated
X function.
X
X Window or kernel estimators can also be used to assess the
X adequacy of a model for data. The assessment is obtained by fitting
X the model within each window and plotting the coefficients obtained
X from these fits versus the independent variable. If the model is
X adequate over the range of the independent variable then the
X coefficients should be constant and the plot should reveal no
X systematic change. This capability is not now in the WINDOW computer
X code.
X
X H. Other Methods
X
X There are other methods for smoothing data. One uses smoothing
X splines instead of fixed width or nearest neighbor estimation.
X Another method uses a different window width for different window
X centers; the widths are chosen to optimize some goodness of fit
X criterion for each window separately. Finally, kernels other than the
X boxcar and biquartic are in common use.
X
X I. References
X
X I. References
X
X Hardle, W. Smoothing Techniques with Implementation in S.
X Springer-Verlag, New York. (1991)
X
X T. J. Hastie and R. J. Tibshirani. Generalized Additive Models.
X Chapman and Hall, London. (1990)
X
END_OF_FILE
if test 32562 -ne `wc -c <'windows.doc'`; then
echo shar: \"'windows.doc'\" unpacked with wrong size!
fi
# end of 'windows.doc'
fi
echo shar: End of shell archive.
exit 0