info       MACRO
) This file contains MacAnova macros used in multivariate analysis
) Copyright (C) 2004 by Christopher Bingham and Gary W. Oehlert
)
) Many use features available only in MacAnova4.10 or later
) Some use features available only in MacAnova 4.13 or later
)) File version 041216
)) All macros are the "new" format (no line count, macro body terminated
)) by %macroname%)
) It contains the following macros.
) Descriptive statistics
)        covar          distcomp       groupcovar     standardize
) Mean testing related macros
)        hotellval      hotell2val
) MANOVA hypothesis testing related macros
)        cumpillai      cumtrace       cumwilks       seqF
) Factor analysis related macros
)        facanal        glscrit *      mlcrit *       ulscrit *
)        glsfactor      glsresids +    stepuls **     stepgls **
)        stepml **      ulsfactor      ulsresids +
) Obsolete factor analysis related macro
)        goodfit ++ 
) Discriminant analysis related macros
)        discrim        jackknife      discrimquad    probsquad
) Subset discriminant analysis related macros
)        dastepsetup    dastepstatus   dasteplook     daentervar
)        daremovevar    davarid        dastepselect   dascreen
) Obsolete subset discriminant analysis related macros
)        compf +        forstep +      backstep +
) Plotting macros
)        chiqqplot
) Other miscellaneous macros
)        rmvnorm        mulvarhelp
) *  Used by facanal()
) ** Obsolete except for use in getting starting values in facanal()
) +  Used by glsfactor() or ulsfactor()
) ++ Obsolete
)) 001210 added subtopics to help
)) 030114 moved standardize() here from MacAnova.mac and wrote help
))        added macros facanal() and supporting macros glscrit(),
))        mlcrit() and ulscrit() and wrote help for them
))        renamed mvngen() to rmvnorm() (left in a stub)
))        added keys for most topics
)) 030814 added subtopic titles
)) 030922 minor correction to help
)) 031224 exact match of names on dasteplook() is no longer required
))        new macro davarid(), used by daentervar() and daremovevar()
)) 031226 dastepsetup() uses pushmodel() before invoking manova()
))        New macros dastepselect(), black box driver for dastepsetup(),
))        daentervar() and daremovevar() and dascreen(), a macro
))        for screening all subsets
)) 040106 added distcomp() keywords 'center' and 'quadform'
)) 041012 added keyword 'silent' and replaced 'verbose' by 'quiet' in
))        dastepselect(); changed default for 'mbest', and shape
))        and labelling of output in dascreen()
)) 041018 Reorganized cumtrace() and changed default operation; added
))        4-th term and method to choose how many terms to use when
))        not specified
)) 041026 Fixed bug in cumtrace() introduced on 041018; removed extra ')'
))        in distcomp()
)) 041111 Added keyword phrase 'nsteps:n' to stepuls(), stepgls(),
))        stepml(); use eigen() instead of releigen() in stepgls()
))        and stepml()
)) 041115 Added keyword phrase 'presteps:n' to facanal() specifying
))        starting value obtained by n steps of stepuls(), stepgls()
))        or stepml(), default is presteps:1; illegal with start:psi0
)) 041201 Added usage lines to dastepselect() help.
)) 041202 Cosmetic change to dasteplook()
)) 041216 Changed dastepselect() to use F-statistics to choose candidates
))        for entering or removing a variable.  Previously it used P-values
))        which caused indeterminacy when there were several values of
))        P = 0.
%info%

===> standardize <===
standardize     MACRO DOLLARS
) usage:
)   ynew <- standardize(y [,locations [,scales]])
)    y            REAL vector of length n or REAL n by p matrix
)    locations    REAL scalar or vector or row vector with
)                 length(locations)= ncols(y) (default describe(y,means:T))
)                 scalar is equivalent to rep(locations, p)
)                 MISSING elements replaced by mean
)    scales       REAL scalar vector or row vector with length(scales)
)                 = ncols(y) (default describe(y,stddev:T))
)                 scalar is equivalent to rep(scales, p)
)                 MISSING elements replaced by standard deviation
)     ynew        REAL vector or matrix the same size and shape as y
)
)) Version 991219 removed $$
)) 030114 allow scalar locations and/or scales
))        allow MISSING elements in locations and/or scales
# usage: ynew <- $S(y [,locations [,scales]])
if($v > 3 || $v == 0){
	error("usage: $S(y [,locs [,scales])",macroname:F)
}
@y <- argvalue($1,"$1","real matrix")
@p <- ncols(@y)
if($v < 3){
	@stats <- describe(@y,mean:T, stddev:T)
}
@loc <- if($v >= 2){
	matrix(argvalue($2 ,"argument 2","real matrix"))
}else{
	matrix(@stats$mean)
}

@scale <- if($v == 3){
	matrix(argvalue($3,"argument 3","real matrix"))
}else{
	matrix(@stats$stddev)
}
if (min(dim(@scale)) > 1){
	error("argument 2 not REAL vector or row vector")
}
if (min(dim(@loc)) > 1){
	error("argument 2 not REAL vector or row vector")
}
if (isscalar(@loc)){
	@loc <- rep(@loc,@p)
} elseif(length(@loc) != @p){
	error("length($3) != ncols($1)")
}
if (isscalar(@scale)){
	@scale <- rep(@scale,@p)
} elseif(length(@scale) != @p){
	error("length($3) != ncols($1)")
}
if (anymissing(@loc)){
	if (!isdefined(@stats)){
		@stats <- describe(@y, mean:T, stddev:T)
	}
	@J <- ismissing(@loc)
	@loc[@J] <- @stats$mean[@J]
	delete(@J)
}	
if (anymissing(@scale)){
	if (!isdefined(@stats)){
		@stats <- describe(@y, mean:T, stddev:T)
	}
	@J <- ismissing(@scale)
	@scale[@J] <- @stats$stddev[@J]
	delete(@J)
}	
if(isdefined(@stats)){
	delete(@stats)
}
(delete(@y,return:T) - vector(delete(@loc,return:T))')/\
	vector(delete(@scale,return:T))'
%standardize%

===> covar <===
covar              MACRO DOLLARS
) Macro to compute sample covariance matrix and mean vector
) Usage:
)    stats <- covar(x)
)    x             REAL N by p matrix with no MISSING values
)    stats         structure(n:N, mean:ybar, covariance:S)\
)     N            sample size = nrows(x)
)     ybar         1 by p REAL row vector of sample variable means
)     S            p by p REAL sample variance matrix
)
)) Version of 000531, uses tabs(y,covar:T,n:T,mean:T)
# usage: d <- $S(x), where x is a REAL matrix with no MISSING values
@R <- tabs(argvalue($1,"argument to $S","real matrix nonmissing"),\
	mean:T,count:T,covar:T)
@R <- structure(n:@R$count[1], mean:@R$mean', covariance:@R$covar)
delete(@R, return:T)
%covar%

===> groupcovar <===
groupcovar            MACRO DOLLARS
) macro to compute group means and pooled sample covariance matrix
) Usage:
)  result <- groupcovar(groups, y)
)   groups        factor or vector of positive integers
)   y             REAL matrixwith nrows(y) = length(groups)
)
)   result        structure(n:N, mean:ybar, covariance:S)
)    N            integer vector of group sizes, length(N) = g =
)                 max(groups)
)    ybar         g by ncols(y) matrix of group mean vectors, with
)                 MISSING values for any empty groups
)    S            pooled sample variance-covariance matrix
)) 000601 new version using tabs(y,means:T,count:T,covar:T)
)) 011024 fixed problem with factor groups with max(groups) < nfactors
@groups <- vector(argvalue($1, "group numbers", "positive integer vector"))
@n <- tabs(,@groups)
if (max(@n) < 2){
	error("no groups have 2 or more cases")
}
@R <- tabs(argvalue($2,"argument to $S","real matrix nonmissing"),\
	delete(@groups, return:T), mean:T, covar:T)
@ng <- length(@n)
@cov <- 0
for(@i,1,@ng){
	if (@n[@i] > 1){
		@cov <-+ (@n[@i]-1)*matrix(@R$covar[@i,,])
	}
}
@cov <-/ sum(@n[@n > 0] - 1)
delete(@ng,@i)
@R <- structure(n:delete(@n,return:T),means:@R$mean, \
	covariance:delete(@cov,return:T))
delete(@R,return:T)
%groupcovar%

===> distcomp <===
distcomp           MACRO DOLLARS
) Macro to compute vector of generalized distances
) Usage:
)   d <- distcomp(y [,center:cent] [,quadform:Q])
)     y            REAL N by p data matrix with no MISSING values
)     cent         REAL vector or row vector of length p with no
)                  MISSING values; default: cent = ybar = sum(y)/N
)     Q            REAL positive definite symmetric p by p matrix,
)                  default: Q = S = ((y-ybar)' %*% (y-ybar))/(N-1)
)     d            REAL vector with length N
)
) d[i] is (y[i,] - cent)' %*% solve(Q) %*% (y[i,] - cent),
) When needed, tabs() or describe() are used to compute ybar and S
)
)) Version of 000826, uses tabs() instead of macro covar()
)) 040105 added keywords 'center' and 'quadform'
# usage: d <- $S(y) where y is a REAL matrix with no MISSING values
@y <- matrix(argvalue($1, "data matrix", "real matrix nonmissing"))
@cent <- keyvalue($K, "cent*","nonmissing real matrix")
@Q <- keyvalue($K,"quad*","nonmissing real square")
@p <- ncols(@y)
if (!isnull(@cent)){
	if (!isvector(@cent)) {
		if (nrows(@cent) != 1) {
			error("center not column or row vector")
		}
		@cent <- @cent'
	}
	if (length(@cent) != @p){
		error("length(center) != ncols(y)")
	}
}
if (isnull(@Q)){
	@Q <- if (isnull(@cent)){
		@stats <- tabs(@y,covar:T,mean:T)
		@cent <- @stats$mean
		delete(@stats,return:T)$covar
	} else {
		tabs(@y,covar:T)
	}	
} else {
	if (nrows(@Q) != @p){
		error("quadform matrix has wrong dimensions")
	}
	@testval <- max(abs(vector(@Q - @Q')))
	if (@testval > 1e-10*max(abs(vector(@Q)))) {
		error("quadratic form matrix not symmetric")
	}
	if (delete(@testval, return:T) != 0){
		@Q <- (@Q + @Q')/2
	}
	@vals <- eigenvals(@Q)
	if (anytrue(@vals[1] <= 0, @vals[@p]/@vals[1] <= 1e-12)) {
		error("quadratic form matrix not positive definite")
	}
	delete(@vals)
	if (isnull(@cent)) {
		@cent <- describe(y, mean:T)
	}
}
@y <-- delete(@cent,return:T)'

@dsq <- vector((@y * (@y %/% delete(@Q,return:T))) %*% rep (1,@p))
if (haslabels(@y)){
	setlabels(@dsq, getlabels(@y,1))
}
delete(@p, @y)
delete(@dsq, return:T)
%distcomp%

===> chiqqplot <===
chiqqplot       MACRO DOLLARS
) Macro to make a chisquare or sqrt(chisquare) QQ plot of generalized
) distance of residuals from 0
) usage:
)   chiqqplot(dsq, df [,sqrt:T] [, symbol:sym])
)     dsq               REAL vector or matrix of squared distances >= 0
)     df                REAL scalar > 0
) or
)   chiqqplot(y  [,sqrt:T] [, symbol:sym])
)     y                 REAL data matrix with no MISSING values from which
)                       a vector dsq of squared generalized distances from
)                       the mean vector and df = ncols(y).
)
) Ordered values of dsq are plotted against equally spaced quantiles
) of chi-squared on df degrees of freedom.  When dsq is a matrix, each
) column is plotted
)
) With 'sqrt:T', sqrt(dsq) is plotted against sqrt(chi-squared quantiles)
)
) With symbol:0, when ncols(d) = 1, the symbols are the case numbers, and
) when ncols(d) > 1, the symbols are the column numbers.  Otherwise,
) symbol:sym works as for plot().
)
) The usual graphics keywords such as 'xlab', 'ylab' and 'title' may also
) be used
))version 000601
#$S(d, p [, sqrt:T] [,symbols:T] [,graphics keyword phrases])
if ($v > 2){
	error("usage: $S(y [keyword phrases]) or $S(d, df [keyword phrases]")
}
if ($v == 1){ # need to compute vector of distances
	if (!ismacro(covar)){
		getmacros(covar,quiet:T,printname:F)
	}
	if (!ismacro(distcomp)){
		getmacros(distcomp,quiet:T,printname:F)
	}
	@d <- matrix(argvalue($1,"data","nonmissing real matrix"))
	@df <- ncols(@d)
	@d <- distcomp(@d)
}else{ # arg 1 is matrix of distance; arg2 is d.f.
	@d <- matrix(argvalue($1,"distances","nonmissing real matrix"))
	@df <- argvalue($2,"degrees of freedom","positive number")
}
@sqrt <- keyvalue($K,"sqrt","TF",default:F)

@n <- nrows(@d)
@m <- ncols(@d)
@J <- grade(@d)

@ch <- vector("\1","\2","\3","\4","\5","\6","\7","\10")[1 + run(0,@m-1)%%8]
@symbols <- keyvalue($K,"symbol*","matrix",default:@ch)
if (!anytrue(ischar(@symbols),\
	alltrue(isreal(@symbols), sum(vector(@symbols!=floor(@symbols))) == 0))){
	error("Value of 'symbols' not CHARACTER or integer")
}
if (alltrue(isscalar(@symbols,real:T),@symbols == 0)){
	@symbols <- if (@m == 1){
		@J
	}else{
		run(@m)'
	}
}
@d <- sort(@d)
@vals <- invchi((run(@n) - .5)/@n,@df)
if (@sqrt){
	@d <- sqrt(@d)
	@vals <- sqrt(@vals)
	@what <- "Chi"
	@ylab <- "Distances"
}else{
	@what <- "Chi-squared"
	@ylab <- "Squared distances"
}
@xlab <- paste(@what,"(",@df,") probability points",sep:"")
@title <- paste(@what, "Q-Q plot of $1")

plot(delete(@vals,return:T), delete(@d,return:T),\
	symbols:delete(@symbols,return:T), $K,\
	xlab:delete(@xlab,return:T),ylab:delete(@ylab,return:T),\
	title:delete(@title,return:T),xmin:0, ymin:0)
%chiqqplot%

===> rmvnorm <===
rmvnorm  macro  dollars
) Macro to compute multivariate normal sample
) Usage:
)  y = rmvnorm(n , p [,mu])
)   n               REAL integer > 0
)   p               REAL integer > 0
)   mu              nonMISSING REAL row or column vector of length p,
)                   default is mu = rep(0,p)
)   y               random REAL n by p matrix with rows i.i.d
)                   MVN(mu, I_p)
)  y = rmvnorm(n, sigma [,mu])
)   sigma           positive p by p definite symmetrix matrix with p > 1
)   y               random REAL n by p matrix with rows i.i.d
)                   MVN(mu, sigma)
)
) Use mu + sigma*rnorm(n) to generate a univariate sample
)
)) Version 991120, C. Bingham, kb@stat.umn.edu
)) 030406 fixed bug when mu not an argument
#$S(n,sigma [,mu])
@n <- argvalue($1,"sample size","positive count")
@sigma <- argvalue($2,"covariance matrix","nonmissing real square")
if (isscalar(@sigma)){
	@p <- argvalue(@sigma,"dimension","positive count")
	@sigma <- NULL
}else{
	if(max(abs(vector(@sigma-@sigma')/vector(@sigma))) > 1e-12){
		error("Covariance matrix not symmetric")
	}
	@sigma <- cholesky(@sigma, nonposok:T)
	if (isnull(@sigma)){
		error("covariance not positive definite")
	}
	@p <- nrows(@sigma)
}
@mu <- if ($v > 2){
	matrix(argvalue($3,"mean vector","nonmissing real matrix"))
}else{
	rep(0,@p)'
}

if (min(dim(@mu)) > 1){
	error("mu not row or column vector")
}

if(length(@mu) != @p){
	error("mean vector length != nrows(sigma)")
}
@y <- matrix(rnorm(@n*@p),@n)
if (!isnull(@sigma)){
	@y <- @y %*% @sigma
}
delete(@n,@sigma,@p)
delete(@y,return:T) + vector(delete(@mu, return:T))'
%rmvnorm%

===> mvngen <===
mvngen       MACRO
) alias for rmvnorm() to maintain backward compatibility
if (!ismacro(rmvnorm)){
	getmacros(rmvnorm, silent:T)
}
rmvnorm($0)
%mvngen%

===> cumpillai <===
cumpillai   macro
) alias for cumtrace(tr,fh,fe,p,pillai:T ...) to compute
) probabilities for Pillai's trace test
)) 031230 written by C. Bingham, kb@umn.edu
#$S(tr, fh, fe, p [,upper:T] [,nterms:m] [,all:T])
if (!ismacro(cumtrace)){
	getmacros(cumtrace, silent:T)
}
cumtrace($0,pillai:T)
%cumpillai%

===> cumtrace <===
cumtrace   macro    dollars
) Macro to compute cumulative distribution of MANOVA trace tests
) T_0^2/fe = trace(solve(E,H)) = sum(vals) (Hotelling)
) V/(fe+fh) = trace(solve(E+H,H)) = sum(vals/(1+vals)) (Pillai)
) where vals is the vector of eigenvalues of H relative to E.
)
) Usage:
)  cumP <- cumtrace(tr, fh, fe, p [,upper:T] [,nterms:m or useF:T])
)  cumP <- cumtrace(tr, fh, fe, p, pillai:T [,upper:T] [,nterms:m or useF:T])
)  allF <- cumtrace(tr, fh, fe, p [,pillai:T] [,upper:T], useF:T, all:T)
)  allP <- cumtrace(tr, fh, fe, p [,pillai:T] [,upper:T], all:T)
)   tr          Positive REAL scalar, vector, matrix or array with
)               no MISSING values; see below.
)   fh, fe, p   Positive integer scalars, vectors, matrices or arrays
)   pillai:T    Argument 1 is trace(solve(H+E,H)) = V/(fe+fh) where V =
)               Pillai statistic; otherwise it's trace(solve(E,H)) =
)               T_0^2/fe, T_0^2 = Hotelling statistic
)   m           1, 2, 3 or 4, number of terms of series used; ignored
)               with all:T; illegal with useF:T
)   useF:T      Use approximate F-distribution; illegal with nterms:m
)   cumP        REAL scalar or vector, matrix or array of probabilities
)               with dimensions matching those any non-scalar tr
)               fh, fe, or p.
)   allF        structure(P:cumP, f:fstat, df1:dfnum, df2:dfdenom)
)   allP        structure(P:cumPStr, terms:termStr)
)    cumPStr    structure(P_1:cumP1, P_2:cumP2, P_3:cumP3, P_4:cumP4)
)    termStr    structure(Term_1:t1, Term_2:t2, Term_3:t3, Term_4:t4)
)
)  All non scalar tr, fh, fe or p must have identical dimensions
)
)  Without pillai:T (or with pillai:F),
)       tr = T_0^2/fe = trace(solve(E) %*% H)) = sum(vals)
)  where T_0^2 is Hotelling's generalized T^2
)
)  With pillai:T,
)       tr = V/(fe+fh) = trace(solve(E+H) %*% H) = sum(vals/(1 + vals))
)  where V = Pillai's trace statistic
)
)  H and E are p by p MANOVA hypothesis and error matrices with degrees
)  of freedom fh and fe and vals are the eigenvalues of H relative to E
)  (eigenvalues of solve(E,H))
)
)  With upper:T, upper tail probability is computed (usually what is
)  wanted)
)
)  With neither useF:T or nterms:m, 1 to 4 terms of Fujikoshi's series
)  is used, the choice depending on the size of the terms, except
)  when min(fh,p) = 1, when the exact value based on F is computed.
)
)  The first term of Fujikoshi's series is the usual chi-squared
)  approximation, cumchi(mj*tr, fh*p [,upper:T]), where
)    mj = m2 = fe - p - 1, without pillai:T
)    mj = m3 = fh + fe, with pillai:T
)  The j-th term of the series is O(1/mj^(j-1)).
)
)  With nterms:m, probabilities are computed using m terms of Fujikoshi's
)  series, even when min(p,fh) = 1
)
)  With all:T, but without useF:T, the component P of structure ansP
)  contains the probabilities computed with 1, 2, 3 and 4 terms of
)  Fujikoshi's series; component terms contains the actual terms
)  themselves.
)
)  With useF:T, probabilities are computed from approximations based
)  on the F-distribution which are used in SAS and SPSS which are
)  exact when min(p,fh) = 1.  With all:T, the returned structure
)  contains the F-statistics and their degrees of freedom in addition
)  to the probabilities.  Except when min(p,fh) = 1, the probabilities
)  appear to be less accurate than those from Fujikoshi's series.
)
) Examples
) Cmd> list(Y,groups) # response matrix and factor
) groups          REAL   50    FACTOR with 4 levels
) Y               REAL   50    4
)
) Cmd> manova("Y = groups",silent:T)
)
) Cmd> h <- matrix(SS[2,,]);fh <- DF[2]
)
) Cmd> e <- matrix(SS[3,,]);fe <- DF[3]
)
) Cmd> p <- ncols(Y); vector(fh,fe,p)
) (1)           3          46           4
)
) Cmd> trhotelling <- trace(solve(e,h)); trhotelling
) (1)        0.77
)
) Cmd> trpillai <- trace(solve(h+e,h)); trpillai
) (1)     0.48115
)
) Cmd> cumtrace(trhotelling,fh,fe,p,upper:T) # default
) (1)   0.0047263
)
) Cmd> cumtrace(trpillai,fh,fe,p,pillai:T,upper:T)
) (1)    0.013864
)
) Cmd> cumtrace(trhotelling,fh,fe,p,upper:T,all:T)
) component: P
)   component: P_1
) (1)   0.0016117        [usual chi-squared approximaition]
)   component: P_2
) (1)    0.003814
)   component: P_3
) (1)   0.0047398
)   component: P_4
) (1)   0.0047263        [value returned by default usage above]
) component: terms
)   component: Term_1
) (1)   0.0016117
)   component: Term_2
) (1)   0.0022023
)   component: Term_3
) (1)  0.00092578
)   component: Term_4
) (1) -1.3487e-05
)
) Cmd> cumtrace(trpillai,fh,fe,p,upper:T,pillai:T,nterms:1)
) (1)    0.023213        [usual chi-squared approximation]
)
) Cmd> cumtrace(trhotelling,fh,fe,p,upper:T,useF:T) # F-approximation
) (1)   0.0031321        [F-approximation probability]
)
) Cmd> cumtrace(trpillai,fh,fe,p,upper:T,useF:T,all:T)
) component: P
) (1)    0.081057        [F-approximation probability]
) component: f
) (1)      1.6707
) component: df1
) (1)          12
) component: df2
) (1)         125
)
) From a simulation using 100,000 samples, the P-value for
) trhotelling is estimated as 0.00451 and for trpillai
) is estimated to be 0.01431.  These are quite close to
) the values computed by the default usage of cumtrace() and
) less close to the value computed with useF:T
)
)   See Fujikoshi, Yasunori (1973)  Asymptotic formulas for the
)     distributions of three statistics for multivariate linear
)     hypothesis, Ann. Inst. Statist. Math. 25 423-437
)) 031015 written by C. Bingham, kb@umn.edu
)) 031017 corrected bug
)) 031229 added all:T, other minor changes
)) 031230 added pillai:T to handle Pillai's trace
)) 041018 added order m^-3 term
))        always uses n terms with nterms:n
))        returns terms and partial sums with all:T and without useF:T
))        Otherwise, without useF:T, for items with s > 1,
))        finds appropriate truncation of Fujikoshi's series
)) 041026 Fixed bug that forced using just 1 term
#$S(tr, fh, fe, p [,pillai:T] [,upper:T] [,nterms:m] [,all:T])
@trace <- argvalue($1, "$1", "positive")
@fh <- argvalue($2, "fh", "positive integer")
@fe <- argvalue($3, "fe", "positive integer")
@p <- argvalue($4, "p", "positive integer")

@upper <- keyvalue($K,"up*","TF",default:F)
@nterms <- keyvalue($K,"nterm*","positive count")
@useF <- keyvalue($K,"useF","TF",default:F)
@pillai <- keyvalue($K,"pillai*","TF",default:F)
@all <- keyvalue($K,"all","TF",default:F)

if (@useF) {
	if (!isnull(@nterms)){
		error("'useF:T' illegal with 'nterms:m'")
	}
} else {
	if (isnull(@nterms)){
		@nterms <- -4
	} elseif (@all) {
		print("WARNING: 'nterms' ignored with 'all:T'",macroname:T)
	} elseif (@nterms > 4) {
		error("nterms > 4 not implemented")
	}
}

@names <- vector("trace","fh","fe","p")
@args <- structure(@trace,@fh,@fe,@p)
for (@i,1,3) {
	for(@j,@i+1,4) {
		@ai <- @args[@i]
		@aj <- @args[@j]
		if (!isscalar(@ai) && !isscalar(@aj) &&\
			anytrue(ndims(@ai) != ndims(@aj),\
				    sum(dim(@ai) != dim(@aj)) != 0)) {
			error(paste(@names[@i],"and",@names[@j],"dimension mismatch"))
		}
	}
}

@dim <- 1
for (@i,1,4) {
	if (!isscalar(@args[@i])) {
		@dim <- dim(@args[@i])
		break
	}
}
delete(@i,@j,@ai,@aj,@names,@args)

@N <- prod(@dim)
if (prod(@dim) > 1) {
	if (isscalar(@trace)) {
		@trace <- array(rep(@trace,@N),@dim)
	}
	if (isscalar(@fh)) {
		@fh <- array(rep(@fh,@N),@dim)
	}
	if (isscalar(@fh)) {
		@fe <- array(rep(@fe,@N),@dim)
	}
	if (isscalar(@p)) {
		@p <- array(rep(@p,@N),@dim)
	}
}
delete(@N,@dim)

@df1 <- @p*@fh
@s <- (@fh + @p - abs(@fh - @p))/2 # min(p,fh)
@JF <- if (@useF){
	rep(T,length(@trace))
} elseif (@all || @nterms < 0) {
	rep(F,length(@trace))
} else {
	vector(@s == 1)
}
@JG <- !@JF
									
if (max(@JF) > 0) { # use F approximation like SAS and SPSS
	#@m <- (abs(@fh - @p) - 1)/2
	@n <- (@fe - @p - 1)[@JF]/2
	if (!@pillai) {
		# 2*(s*n+1)*trace(solve(e,h)) / (s^2*(2*m+s+1))
		#         ~= F(s*(2*m+s+1), 2*(s*n+1))
		# and s*(2*m+s+1) = fh*p
		@df2 <- 2*(@s[@JF]*@n + 1)
		@f <- ((@df2/@df1[@JF])*(@trace/@s)[@JF])
	} else {
		@df2 <- @s[@JF]*(2*@n + @s[@JF] + 1)
		@f <- (@df2/@df1[@JF])*(@trace/(@s - @trace))[@JF]
	}
	@probsF <- cumF(@f,@df1[@JF],@df2,upper:@upper)
	if (@all) {
		@result <- structure(P:@probsF, f:@f, df1:@df1, df2:@df2)
	}
	delete(@n,@s,@f,@df2)
}

if (max(@JG) > 0) {
	# use Fujikoshi series
	@mj <- if (!@pillai) {
		(@fe - @p - 1)[@JG] # m2
	} else {
		(@fh + @fe)[@JG]    # m3
	}
	@x <- @mj * @trace[@JG]
	@f <- @df1[@JG]
	@P0 <- @probsG <- cumchi(@x, @f,upper:@upper) # term 1
	
	if (@all || abs(@nterms) > 1) {
		if (@all) {
			@P <- structure(P_1:@P0, P_2:NULL, P_3:NULL, P_4:NULL)
			@terms <- structure(Term_1:@P0, Term_2:NULL, Term_3:NULL, \
								Term_4:NULL)
		}

		@gam <- (@p + @fh + 1)[@JG]
		# Compute term 2
		@P2 <- cumchi(@x, @f + 2, upper:@upper)
		@P4 <- cumchi(@x, @f + 4, upper:@upper)
		@term <- (@f*@gam)*(@P4 - 2*@P2 + @P0)/(4*@mj)
		if (@pillai) {
			@term <- -@term
		}

		if (@all) {
			@P[2] <- @P[1] + @term
			@terms[2] <- @term
		} elseif (@nterms > 1){
			@probsG <-+ @term
		} elseif (@nterms < 0) {
			@term2 <- @term
		}

		if (@all || abs(@nterms) > 2) {
			# Compute term 3
			#Fujikoshi has following in the form (p+fh+1) + p*fh + 2
			@g <- ((@p+1)*(@fh+1))[@JG] + 2

			@h0 <- (3*@f - 8)*@gam^2 + 4*@g
			@h1 <- 12*@f*@gam^2
			@h2 <- 6*(3*@f + 8)*@gam^2
			@h3 <- 4*((3*@f + 16)*@gam^2 + 4*@g)
			@h4 <- (3*@f + 24)*@gam^2 + 12*@g
			# h's should sum to 0

			@P6 <- cumchi(@x, @f + 6,upper:@upper)
			@P8 <- cumchi(@x, @f + 8,upper:@upper)

			@term <- @f*(@h0*@P0 - @h1*@P2 + @h2*@P4 - \
				@h3*@P6 + @h4*@P8)/(96*@mj^2)

			if (@all) {
				@P[3] <- @P[2] + @term
				@terms[3] <- @term
			} elseif (@nterms > 2) {
				@probsG <-+ @term
			} elseif (@nterms < 0) {
				@term3 <- @term
			}
			if (@all || abs(@nterms) > 3) {
				# Compute term 4
				# order fe^-3 terms
				@g0 <- @gam*((@f-4)^2*@gam^2 + \
						 4*(@f-4)*@gam + 4*(@f-4)*(@f+2))
				@g1 <- 2*@f*@gam*@h0
				@g2 <- @f*@gam*((15*@f+40)*@gam^2 + 4*@gam + 4*(@f+2))
				@g3 <- 4*@gam*(5*(@f+4)^2*@gam^2 + \
							    4*(@f+4)*@gam + \
								4*(@f+4)*(@f+2))
				@g4 <- 5*rational(@f,vector(144,40,3))*@gam^3 + \
					   (44*@f+432)*@gam^2 + \
					   4*rational(@f,vector(288,130,11))*@gam + 96*(@f + 2)
				@g5 <- 2*(rational(@f,vector(288,56,3))*@gam^3 + \
					      4*(5*@f + 72)*@gam^2 + \
						  4*rational(@f,vector(216,82,5))*@gam +\
						  96*(@f+2))
				@g6 <- rational(@f,vector(160,24,1))*@gam^3 + \
					   4*(3*@f + 56)*@gam^2 + \
					   4*rational(@f,vector(184,62,3))*@gam +\
					   96*(@f+2)
				# @g0 - @g1 + @g2 - @g3 + @g4 - @g5 + @g6 = 0
				@P10 <- cumchi(@x, @f + 10,upper:@upper)
				@P12 <- cumchi(@x, @f + 12,upper:@upper)

				@term <- @f*(@g0*@P0 - @g1*@P2 + @g2*@P4 - @g3*@P6 +\
						     @g4*@P8 - @g5*@P10 + @g6*@P12)/(384*@mj^3)
				if (@pillai) {
					@term <- -@term
				}

				if (@all) {
					@P[4] <- @P[3] + @term
					@terms[4] <- @term
					@result <- structure(P:delete(@P,return:T),\
										 terms:delete(@terms, return:T))
				} elseif (@nterms > 3) {
					@probsG <-+ @term
				}  elseif (@nterms < 0) {
					# Work backwards, adding terms and earlier terms
					# if a term is smaller than preceding term
					
					@L <- vector(abs(@term) <= abs(@term3))
					if (max(@L) > 0){ # include term 4
						@probsG[@L] <- @probsG[@L] + @term[@L]
					}

					@L <- @L || vector(abs(@term3) <= abs(@term2))
					if (max(@L) > 0){ # include term 3
						@probsG[@L] <- @probsG[@L] + @term3[@L]
					}

					@L <- @L || vector(abs(@term2) <= (abs(1-2*@probsG) + 1)/2)
					if (max(@L) > 0){ # include term 2
						@probsG[@L] <- @probsG[@L] + @term2[@L]
					}

					@L <- vector(@probsG < 0 || @probsG > 1)
					if (max(@L) > 0) {#if use term 1 if out of range
						@probsG[@L] <- @P0[@L]
					}
					delete(@term2,@term3,@L)
				}
				delete(@g0,@g1,@g2,@g3,@g4,@g5,@g6,@P10,@P12)
			}
			delete(@h0,@h1,@h2,@h3,@h4,@g)
			delete(@P6,@P8)
		}
		delete(@gam, @term, @P2, @P4)
	}
	delete(@P0,@x,@mj)
}
if (!@all) {
	@result <- @trace # ensures dimensions are OK
	if (max(@JG) > 0) {
		@result[@JG] <- delete(@probsG, return:T)
	}
	if (max(@JF) > 0) {
		@result[@JF] <- delete(@probsF, return:T)
	}
}
delete(@trace,@fh,@fe,@upper,@nterms, @useF, @pillai, @df1)
delete(@result,return:T)
%cumtrace%

===> cumwilks <===
cumwilks   macro  dollars
) Macro to compute cumulative distribution of MANOVA LR
) statistic lambda* = det(E)/det(H+E) = 1/prod(1+vals)
) It uses either an F-approximation (exact for min(p,fh) <=2) or
) a series in powers of 1/m1^2, m1 = fe - (p - fh + 1)/2.
) Usage:
)  cumP <- cumwilks(lambdastar, fh, fe, p [,upper:T] [,useF:F] [,nterms:m])
)  allF <- cumwilks(lambdastar, fh, fe, p, all:T, [,upper:T])
)  allP <- cumwilks(lambdastar, fh, fe, p, all:T, useF:F [,upper:T])
)   lambdastar  Positive REAL scalar, vector, matrix or array with
)               no MISSING values whose elements should be
)               det(E)/det(H+E) = 1/prod(1+vals) = (LR)^(2/N), where
)               vals are the relative eigenvalues of H relative to E
)               and LR is the unadjusted likelihood ratio statistic
)               -N*log(lambdastar) is the usual -2*log(lambda) form
)   fh, fe, p   Positive integer scalars, vectors, matrices or arrays,
)               designating hypothesis and error degrees of freedom and
)               the dimension
)   useF:F      Probabilities are computed from a series involving
)               chi-square rather than from F
)   nterms:m    m = 1, 2, 3, 4 = number of terms of series involving
)               chi-square; implies useF:F; not legal with useF:T
)   all:T       Output is structure with additional information
)   cumP        REAL scalar or vector, matrix or array of lower (default)
)               pr upper tail probabilities with dimensions matching
)               those of any non-scalar lambdastar, fh, fe, or p.
)   allF        structure(P:cumP, f:fstats, df1:numdf, df2:denomdf)
)   allP        structure(P:cumPStr, terms:termStr)
)    cumPStr    structure(P_1:cumP1, P_2:cumP2, P_3:cumP3, P_4:cumP4)
)    termStr    structure(Term_1:t1, Term_2:t2, Term_3:t3, Term_4:t4)
)
)  All non scalar lambdastar, fh, fe or p must have identical dimensions
)
)  With upper:T, upper tail probability is computed (usually not
)  what's wanted since small lamdastar is evidence against H_0)
)
)  With useF:T or without either useF:F or nterms:m
)   P is computed from Rao's formula approximating the distribution
)   of (1 - lambdastar^(1/t))/lambdastar^(1/t) by that of a multiple of F,
)   where t is computed from fh, fe and p.  This is exact when
)   s = min(p,fh) <= 2.
)
)   With all:T, ans = structure(P:P,f:fstat,df1:df1,df2:df2)
)
)  With nterms:m
)   P is computed from m terms of a series computing probabilities
)   for x = -m1*log(lambdastar), m1 = fe - (p - fh + 1)/2.
)
)   With m = 1, x is treated as chi-squared(p*fh), the standard
)   approximation.  With m = 2, 3 or 4, the chi-squared probability is
)   adjusted by m-1 additional terms of magnitudes O(1/m1^2), O(1/m1^4),
)   and O(1/m1^6).
)
)   The series is always used with nterms:m, even when s = min(p,fh) <= 2
)
)  With useF:F but without nterms:m
)   Probablilities are computed from F when s = min(p,fh) <= 2; otherwise
)   they are computed from the series, using from 1 to 4 terms, the number
)   chosen depending on the relative sizes of the terms
)
)   With all:T, allP = structure(P:cumPStr, terms:termStr) where cumPstr =
)   structure(P_1:cumP1, P_2:cumP2, P_3:cumP3, P_4:cumP4) and termStr =
)   structure(Term_1:t1, Term_2:t2, Term_3:t3, Term_4:t4), where cumP1,...,
)   cumP4 are the 1, ...,4 term approximations and t1,...,t4 are the
)   terms themselves.
)
) Example,
)  Cmd> manova("Y=groups",silent:T)
)
)  Cmd> fh <- DF[2]; fe <- DF[3], p <- ncols(Y)
)
)  Cmd> vector(fh,fe,p)
)  (1)           3          46           4
)
)  Cmd> h <- matrix(SS[2,,]); e <- matrix(SS[3,,])
)
)  Cmd> lambdastar <- det(e)/det(e + h);lambdastar
)  (1)      0.54836
)
)  Cmd> cumwilks(lambdastar,fh,fe,p) # Rao's F-approximation
)  (1)    0.0077151
)
)  Cmd> cumwilks(1/prod(1+releigenvals(h,e)),fh,fe,p) # same
)  (1)    0.0077151
)
)  Cmd> cumwilks(lambdastar,fh,fe,p,all:T)
)  component: P
)  (1)   0.0077151
)  component: f
)  (1)      2.4232
)  component: df1
)  (1)          12
)  component: df2
)  (1)      114.06
)
)  Cmd> cumwilks(lambdastar,fh,fe,p,useF:F) # Rao/Box series
)  (1)    0.0077152
)
)  Cmd> cumwilks(lambdastar,fh,fe,p,nterms:1) # use 1st term (chisq)
)  (1)    0.0076322
)
)  Cmd> cumwilks(lambdastar,fh,fe,p,nterms:2) # use two terms
)  (1)    0.0077147
)
)  Cmd> cumwilks(lambdastar,fh,fe,p,useF:F,all:T)
)  component: P
)    component: P_1
)  (1)   0.0076322
)    component: P_2
)  (1)   0.0077147
)    component: P_3
)  (1)   0.0077152
)    component: P_4
)  (1)   0.0077152
)  component: terms
)    component: Term_1
)  (1)   0.0076322
)    component: Term_2
)  (1)  8.2563e-05
)    component: Term_3
)  (1)   4.444e-07
)    component: Term_4
)  (1)  2.1917e-09
)
)   See the following references
)    Fujikoshi, Yasunori (1973)  Asymptotic formulas for the
)     distributions of three statistics for multivariate linear
)     hypothesis, Ann. Inst. Statist. Math. 25 423-437
)    Box, G.E.P (1949) A general distribution theory for a class
)     of likelihood criteria, Biometrika 36, 317-346
)    Rao, C.R. (1948) Tests of significance in multivariate analysis
)     Biometrika 35, 467-480
)) 031015 written by C. Bingham, kb@umn.edu
)) 031229 minor modifications and change of defaults; it no longer
))        recognizes 'series:T' but does recognize 'useF:F'
))        With nterms:m and all:T, result has all three terms
)) 031230 Except with all:T, always uses exact distribution based on
))        F when min(p,fh) = 1
)) 041019 Corrected formulas involved in O(1/m2^4) term from Box
)) 041020 Added O(1/m2^3) term from Box and modified variable names
))        to be closer to Box; @a2, @a4, @a6 are his alpha_2, alpha_4,
))        alpha_6; @a2 was previously @b1 and @a4 + @a2^2/2 was @b2
))        With all:T and useF:F, result is structure with probabilities
))        with O(1/m2^2), O(1/m2^4), O(1/m2^6), and O(1/m2^8)
))        accuracy and the O(1), O(1/m2^2), O(1/m2^4), O(1/m2^6) terms
))        With nterms:m, always uses series even for min(p.fh) <= 2
#$S(lambdastar, fh, fe, p [,upper:T] [,chisq:T [,nterms:m]])
@lamstar <- argvalue($1, "lambdaStar", "positive")
@fh <- argvalue($2, "fh", "positive integer")
@fe <- argvalue($3, "fe", "positive integer")
@p <- argvalue($4, "p", "positive integer")
@upper <- keyvalue($K,"upper","TF",default:F)
@useF <- keyvalue($K,"useF","TF")
@nterms <- keyvalue($K,"nterm*","positive count")
@all <- keyvalue($K,"all","TF",default:F)
if (!isnull(@nterms)){
	if (@all) {
		print("WARNING: 'nterms' ignored with 'all:T'",macroname:T)
	}
	if (isnull(@useF)){
		@useF <- F
	} elseif (@useF) {
		error("useF:T illegal with nterms:m")
	}
	if (@nterms > 4) {
		error("nterms > 4 not implemented")
	}
} elseif (isnull(@useF)){
	@useF <- T
} elseif (!@useF) {
	@nterms <- (2*@all - 1)*4
}

@names <- vector("lambda*","fh","fe","p")
@args <- structure(@lamstar,@fh,@fe,@p)
for (@i,1,3) { # Check dimensions of 1st 4 arguments
	for(@j,@i+1,4) {
		@ai <- @args[@i]
		@aj <- @args[@j]
		if (!isscalar(@ai) && !isscalar(@aj) &&\
			anytrue(ndims(@ai) != ndims(@aj),\
					 sum(dim(@ai) != dim(@aj)) != 0)) {
			error(paste(@names[@i],"and",@names[@j],"dimension mismatch"))
		}
	}
}
@dim <- 1
for (@i,1,4) {
	if (!isscalar(@args[@i])) {
		@dim <- dim(@args[@i])
		break
	}
}
delete(@i,@j,@ai,@aj,@names,@args)
@N <- prod(@dim)
if (prod(@dim) > 1) { # replicated any scalar args
	@lamstar <- if (isscalar(@lamstar)) {
		rep(@lamstar,@N)
	} else {
		vector(@lamstar)
	}
	@fh <- if (isscalar(@fh)) {
		rep(@fh,@N)
	} else {
		vector(@fh)
	}
	@fe <- if (isscalar(@fh)) {
		rep(@fe,@N)
	} else {
		vector(@fe)
	}
	@p <- if (isscalar(@p)) {
		rep(@p,@N)
	} else {
		vector(@p)
	}
}
@Jbad <- @fe < @p
@N1 <- @N - sum(@Jbad)
if (@N1 != @N) {
	if (@N1 == 0) {
		error("all values have fe < p")
	}
	print(paste("WARNING:",@N - @N1,"values have fe < p"),macroname:T)
	@lamstar <- @lamstar[!@Jbad]
	@fh <- @fh[!@Jbad]
	@fe <- @fe[!@Jbad]
	@p <- @p[!@Jbad]
}

@df1 <- @p*@fh
@m1 <- @fe - (@p - @fh + 1)/2

@JF <- if (@useF) {
	rep(T,@N1)
} elseif (@all || @nterms > 0) {
	rep(F,@N1)
} else { # useF:F without nterms:m
	vector(@p <= 2 || @fh <= 2) # cases where F is exact
}
@JG <- !@JF
 # @JF is selector for values where F will be used
 # @JG is selector for values where chi-square or series will be used

if (max(@JF) > 0) { # use Rao's formula
	@df1F <- @df1[@JF]
	@t <- rep(1,@N1)

	@J <- (@p != 1 && @fh != 1)[@JF]
	if (max(@J) > 0) {
		@t[@J] <- sqrt((@df1F^2 - 4)[@J]/(@p^2 + @fh^2 - 5)[@JF][@J])
	}

	@df2F <- (@m1[@JF]*@t - @df1F/2 + 1)

	@fstats <- (@df2F/@df1F)*(@lamstar[@JF]^(-1/@t) - 1)
	@probsF <- cumF(@fstats,@df1F,@df2F,upper:!@upper)
	delete(@t)

	if (@all) {
		if (@N == @N1) {
			@fstats <- array(@fstats, @dim)
			@probsF <- array(@probsF, @dim)
			@df1F <- array(@df1F, @dim)
			@df2F <- array(@df2F, @dim)
		} else {
			@tmp <- @fstats
			@fstats <- array(rep(?,@N),@dim)
			@fstats[!@Jbad] <- @tmp

			@tmp <- @probsF
			@probsF <- @fstats
			@probsF[!@Jbad] <- @tmp

			@tmp <- @df1F
			@df1F <- @fstats
			@df1F[!@Jbad] <- @tmp

			@tmp <- @df2F
			@df2F <- @fstats
			@df2F[!@Jbad] <- delete(@tmp,return:T)
		}			
		@result <- structure(P:@probsF,f:@fstats, df1:@df1F,df2:@df2F)
	}
	delete(@fstats,@df1F,@df2F,@J)
}

if (max(@JG) > 0) { # use chi-squared approximation or Fujikoshi series
	@f <- @df1[@JG]
	@m1 <- @m1[@JG]
	@fh <- @fh[@JG]
	@p <- @p[@JG]
	@x <- -@m1 * log(@lamstar[@JG])
	@probsG <- cumchi(@x, @f,upper:!@upper)

	if (@all){
		@tmp <- array(rep(?,@N),@dim)
		@tmp[!@Jbad] <- @probsG
		@P <- structure(P_1:@tmp, P_2:NULL, P_3:NULL, P_4:NULL)
		@terms <- structure(Term_1:@tmp, Term_2:NULL, Term_3:NULL, \
							Term_4:NULL)
	}

	if (abs(@nterms) > 1 || @all) {
		#@Pk is cumchi(@x, @f + k, upper:!@upper),k=0,2,4,...
		@P0 <- @probsG
		@P4 <- cumchi(@x, @f + 4,upper:!@upper)
		@a2 <- @f*(@p^2 + @fh^2 - 5)/48
		@term <- @a2*(@P4 - @P0)/@m1^2
		if (@nterms < 0) {
			@term2 <- @term
		} else {
			@probsG <-+ @term
		}
		if (@all) {
			@tmp[!@Jbad] <- @probsG
			@P[2] <- @tmp
			@tmp[!@Jbad] <- @term
			@terms[2] <- @tmp
		}
		if (@all || abs(@nterms) > 2) {
			@P8 <- cumchi(@x, @f + 8,upper:!@upper)
			# Rao 1948 incorrectly had 150 in place of 159
			@a4 <- @f*(3*(@p^4+@fh^4) + 10*@f^2 - 50*(@p^2+@fh^2) + 159)/1920

			@term <- ((@a4 + @a2^2/2)*(@P8-@P0) - @a2^2*(@P4-@P0))/@m1^4
			if (@nterms < 0) {
				@term3 <- @term
			} else {
				@probsG <-+ @term
			}

			if (@all) {
				@tmp[!@Jbad] <- @probsG
				@P[3] <- @tmp
				@tmp[!@Jbad] <- @term
				@terms[3] <- @tmp
			}
			if (@all || abs(@nterms) > 3) {
				@P12 <- cumchi(@x, @f + 12,upper:!@upper)
				@a6 <- @f*(3*(@p^6 + @fh^6) - \
						 105*(@p^4 + @fh^4) + \
						1113*(@p^2 + @fh^2) + \
							 (21*@p^2 - 350 + 21*@fh^2)*@f^2 - \
						2995)/16128
				@term <- (@a6*(@P12 - @P0) + \
						 2*@a2*@a4*(@P12 - @P8 - @P4 + @P0) + \
						 @a2^3*(@P12 - 3*@P8 + 3*@P4 - @P0)/6)/@m1^6
				if (@nterms < 0) {
					# use from 1 to 4 terms
					# Trim off terms that are larger than their predecessor
					@term4 <- @term
					@L <- abs(@term4) <= abs(@term3)
					if (max(@L) > 0) {
						@probsG[@L] <- @probsG[@L] + @term4[@L]
					}
					@L <- @L || abs(@term3) <= abs(@term2)
					if (max(@L) > 0) {
						@probsG[@L] <- @probsG[@L] + @term3[@L]
					}
					@L <- @L || abs(@term2) <= (abs(1-2*@P0) + 1)/2
					if (max(@L) > 0) {
						@probsG[@L] <- @probsG[@L] + @term2[@L]
					}
					@L <- @probsG < 0 || @probsG > 1
					if (max(@L) > 0) {#if use term 1 if out of range
						@probsG[@L] <- @P0[@L]
					}
					delete(@term2, @term3, @term4)
				} else {
					@probsG <-+ @term
				}
				if (@all) {
					@tmp[!@Jbad] <- @probsG
					@P[4] <- @tmp
					@tmp[!@Jbad] <- @term
					@terms[4] <- delete(@tmp, return:T)
				}
				delete(@P12, @a6)
			}
			delete(@P8, @a4)
		}
		delete(@P0, @P4, @a2, @term)		
	}

	if (@all) {
		@result <- structure(P:delete(@P, return:T), \
							 terms:delete(@terms, return:T))
	}
	delete(@x, @f)
}

delete(@fh,@fe,@p,@upper,@df1,@m1,@useF,@nterms)

if (!delete(@all,return:T)) {
	@tmp <- rep(0,@N1)
	if (max(@JF) > 0) {
		@tmp[@JF] <- delete(@probsF, return:T)
	}
	if (max(@JG) > 0) {
		@tmp[@JG] <- delete(@probsG, return:T)
	}
	@result <- array(rep(?,@N),@dim)
	@result[!@Jbad] <- delete(@tmp, return:T)
}
delete(@lamstar,@JF,@JG,@N,@N1,@dim)
delete(@result,return:T)
%cumwilks%

===> seqF <===
seqF         MACRO DOLLARS
) Macro to compute sequential F-statistics for MANOVA hypothesis
) Usage:
)   result <- seqF(hypTerm [,errTerm ][,model:Model] [,order:Jorder])
)    hypTerm  scalar CHARACTER term name or integer term number
)    errTerm  scalar CHARACTER term name or integer term number;
)             default = number of last term in model
)    Model    CHARACTER scalar specifying a manova() model;
)             default is STRMODEL
)    Jorder   Length p vector of positive integers containing a
)             permutation of 1, 2, ..., p = dimension specifying
)             variable order for computing sequential F-statistics;
)             default is Jorder = run(p)
)    result   structure(f:fstats,fh:fh,fe:fe), f, fh and fe
)             vectors of F-statistics, hypothesis DF and error
)             DF,each of length ncols(y), where y is the matrix
)             of response variables
)
)) Written by C. Bingham (kb@umn.edu)
)) 011022: modified so errTerm could be omitted
)) 031021: added keyword 'order'
)) 031022: use labels on response matrix
)) 031228: minor changes
)) 041011: more minor changes
# $S(hypterm [, errterm] [,model:Model] [,order:Jorder])
if ($v == 0 || $v > 2) {
	error("Usage is $S(hypterm [, errterm] [,model:Model])")
}
@model <- keyvalue($K,"model","string")
if (isnull(@model)){
	if (!isscalar(STRMODEL,char:T)){
		error("No model for $S and STRMODEL doesn't exist")
	}
	if (!isvector(DF, real:T) || !isreal(SS) || !isvector(TERMNAMES, char:T)){
		@model <- STRMODEL
	}
}
@Jorder <- keyvalue($K,"order","positive integer vector")

@hypterm <- argvalue($1,"hypothesis term")
if (!isscalar(@hypterm,real:T) && !isscalar(@hypterm,char:T)){
	error("argument 1 to $S must be positive integer or term name")
}

if ($v > 1){
	@errterm <- argvalue($2,"error term")
	if (!isscalar(@errterm,real:T) && !isscalar(@errterm,char:T)){
		error("argument 2 must be positive integer or error term name")
	}
} else {
	@errterm <- 0
}

@pushed <- !isnull(@model) && isfunction(pushmodel)
if (@pushed) {
	@pushed <-  pushmodel(canpush:T)
}

if (!isnull(@model)){
	if (@pushed){
		pushmodel()
	}
	manova(@model,silent:T)
}

@nterms <- length(TERMNAMES)

if (ischar(@hypterm)) {
	@hypterm <- match(@hypterm,TERMNAMES,0)
	if (@hypterm == 0){
		error(paste("$1 is not the name of a term in the model"))
	}
} elseif (!isreal(@hypterm, positive:T, integer:T) || @hypterm > @nterms){
	error(paste("$1 not positive integer between 1 and",@nterms))
}

if (ischar(@errterm)){
	@errterm <- match(@errterm,TERMNAMES,0)
	if (@errterm == 0){
		error(paste("$2 is not the name of a term in the model"))
	}
} elseif (@errterm < 0 || !isreal(@errterm,integer:T) || @errterm > @nterms){
	error(paste("$2 not an integer between 1 and",@nterms))
} elseif (@errterm == 0){
	@errterm <- @nterms
} else {
	argvalue(@errterm, "positive count")
}

@e <- matrix(SS[@errterm,,]) # e
@hpluse  <- matrix(SS[@hypterm,,]) + @e # h + e
@p <- nrows(@e)
@fh <- rep(DF[@hypterm], @p)
@fe <- DF[@errterm] - run(@p) + 1
if (@pushed) {
	popmodel()
}
delete(@model)

@labs <- if (@p > 1 && haslabels(@e)) {
	getlabels(@e,1)
} else {
	NULL
}

@fstats <- rep(0, @p)
if (!isnull(@Jorder)) {
	if(anytrue(length(@Jorder) != @p,sum(sort(@Jorder) != run(@p)) != 0)) {
		error("value of 'order' not permutation of run(dimension)")
	}
	@hpluse <- @hpluse[@Jorder,@Jorder]
	@e <- @e[@Jorder,@Jorder]
	if (isnull(@labs)) {
		@labs <- getlabels(vector(rep(0,@p),labels:"y"))
	}
	@labs <- @labs[@Jorder]
}
if (!isnull(@labs)) {
	setlabels(@fstats, @labs)
	setlabels(@fh, @labs)
	setlabels(@fe, @labs)
}

delete(@hypterm,@errterm,@nterms,@labs)

for(@i,1,@p){
	if (@i > 1){
		@hpluse <- swp(@hpluse, @i-1)
		@e <-  swp(@e, @i-1)
	}
	@fstats[@i] <- (@fe[@i]/@fh[@i])*(@hpluse[@i,@i] - @e[@i,@i])/@e[@i,@i]
}

delete(@hpluse,@e,@p,@i)

structure(f:delete(@fstats,return:T),fh:delete(@fh,return:T),\
	fe:delete(@fe,return:T))
%seqF%


===> discrim <===
discrim      MACRO DOLLARS
) Macro to compute coefficients and offsets for multi-group linear
) classification, assuming Multivariate normal populations with equal
) variances matrices.
) Usage:
)  result <- discrim(groups, y)
)    groups          vector of positive integers defining membership
)                    in g non-empty groups
)    y               N by p REAL matrix with no MISSING values, with
)                    N = nrows(groups)
)    result          structure(coefs:L, addcon:C), L p by g matrix,
)                    C 1 by g vector where linear discriminant function for
)                    group j is y' %*% L[,j] + C[j]
)
)) 000825 uses tabs() instead of macro covar
# usage: $S(groups, y), N by 1 vector groups, N by p matrix y
if($v != 2){
	error("usage is $S(groups, y)",macroname:F)
}
@groups <- argvalue($1,"group numbers","positive integer vector")
@y <- argvalue($2, "data", "nonmissing matrix")

@g <- max(@groups) # number of groups
@p <- ncols(@y) # number of variables
@n <- nrows(@y) # number of cases
if (length(@groups) != @n){
	error("length of group index differs from number of cases")
}
if (@n <= @p){
	error("number of cases not larger than dimension")
}

@varlabs <- if (haslabels(@y)){
	getlabels(@y,2)
}else{
	"@"
}
@facname <- "$1"
if (!isname(@facname) || match("@*",@facname,0,exact:F) != 0){
	@facname <- "Group "
}

@ybar <- matrix(rep(0,@p*@g),@p,labels:structure(@varlabs,@facname))
@spooled <- matrix(rep(0,@p*@p),@p,labels:structure(@varlabs,@varlabs))

@stuff <- tabs(@y,@groups,covar:T,count:T,mean:T)
@N <- vector(@stuff$count[,1])
@ybar <- matrix(@stuff$mean',labels:structure(@varlabs,@facname))
@spooled <-\
	matrix(sum((@N-1)*matrix(@stuff$covar,length(@N))),@p,\
		labels:structure(@varlabs,@varlabs))/sum(@N-1)
delete(@stuff,@N,@y,@p)

@coefs <- solve(delete(@spooled, return:T), @ybar) # p x g
@addcon <- -0.5*sum(@coefs * delete(@ybar, return:T)) # 1 x g
structure(coefs:delete(@coefs,return:T), addcon:delete(@addcon,return:T))
%discrim%

===> jackknife <===
jackknife          MACRO   dollars
) Macro to do jackknife validation of linear discrimination
) Usage
)  probs <- jackknife(groups,y [,prior:P])
)   groups             vector of positive integers defining membership
)                      in g = max(groups) groups, each with at least two
)                      cases
)   y                  n by p REAL data matrix with n = length(groups)
)   P                  REAL positive vector of g "prior probabilities",
)                      default = rep(1/g,g)
)
)   probs              n by g+1 matrix
)
)   Column j of probs, 1 <= j <= g contains the estimated posterior
)   probabilities of group j for each case
)   Column g+1 of probs contains the leave-one-out classification of
)   each case
)
)) 991120 use down dating formula for computing coefficients
)) Version of 991120
# usage: $S(groups, y [,prior:P]), factor groups, matrix y, positive vector P
if ($v != 2 || $k > 1){
	error("Usage: $S(groups, y [,prior:P]), factor groups, matrix y, positive vector P",\
		macroname:F)
}
@groups <- argvalue($1, "group numbers","positive integer vector")
@y <- argvalue($2, "data","nonmissing matrix")

@g <- max(@groups)
@P <- keyvalue($K,"prior","positive vector", default:rep(1,@g)/@g)
if (length(@P) != @g){
	error("number of prior probabilities != number of groups")
}

@n <- tabs(,@groups)
if(min(@n) < 2){
	error("min(group sizes) < 2")
}
@N <- nrows(@y)
@p <- ncols(@y)
@fe1 <- @N - 1 - @g

 #compute mean vectors and pooled error matrix
@ybar <- matrix(rep(0,@p*@g),@p)#,labels:structure(@varlabs,@facname))
@E <- matrix(rep(0,@p*@p),@p)#,labels:structure(@varlabs,@varlabs))

for(@i,run(@g)){
	@yi <- @y[@groups == @i,]	
	@ybar1 <- sum(@yi)/@n[@i]
	@ybar[,@i] <- @ybar1
	@yi <-- @ybar1
	@E <-+ @yi %c% @yi
}

@n1 <- @n - 1
@nn1 <- (@n*@n1) # sqrt(n(n-1))

@probs <- matrix(rep(0,@N*@g),@g) #g by N
@coefs <- solve(delete(@E, return:T),hconcat(@ybar,dmat(@p,1)))
@Einv <- @coefs[,-run(@g)] # p by p
@coefs <- @coefs[,run(@g)] # p by g
for(@i,run(@N)){
	@ybar1 <- @ybar
	@j <- @groups[@i]
	@yi <- @y[@i,]
	@d <- (@yi' - @ybar[,@j])/@n1[@j]
	@ybar1[,@j] <- @ybar[,@j] - @d
	@tmp <-  @Einv %*% @d
	@c <- 1 - @nn1[@j] * @tmp %c% @d
	@coefs1 <- @coefs + @nn1[@j] * @tmp * (@d %c% @coefs)/@c
	@coefs1[,@j] <- @coefs1[,@j] - @tmp/@c
	@probs[,@i] <- @fe1*vector(@yi %*% @coefs1  - .5*sum(@ybar1 * @coefs1))
}

@probs <- @P * exp(@probs - @probs[1,])
@probs <- vconcat(@probs/sum(@probs),grade(@probs)[@g,])'
delete(@groups, @y, @n, @g, @N, @n1, @nn1, @fe1, @tmp, @i, @j,\
	@Einv,@P,@ybar,@coefs,@ybar1,@coefs1)
delete(@probs, return:T)
%jackknife%

discrimquad macro dollars
) Macro to compute coefficients for quadratic discrimination,
) assuming multivariate normal.
)
) Usage:
)  quadfn <- discrimquad(groups, y)
)   groups        length N vector of positive integers defining
)                 samples from g populations of p-dimensional vectors.
)                 It is an error if the smallest sample size <= p
)   y             N by p REAL data matrix with no MISSING elements
)   quadfn        structure(Q:q, L:l, addcon:c, grandmean:ybar)
)     q           structure of g p by p REAL matrices,
)     l           structure of g length p REAL vectors,
)     c           length g REAL vector
)     ybar        length p REAL vector of grand means
)
) If x is a length p vector to be classified, the quadratic score for
) group j is
)   qs[j] = (x-ybar)' %*% q[,j] %*% (x-ybar) + (x-ybar)' %*% l[,j] + c[j]
) If P is a length g vector of prior probabilities and qs is the vector
) of scores, the vector of posterior probabilities can be computed as
) P*exp(qs-K)/sum(P*exp(qs-K)), where K is a scalar such as max(qs)
) chosen to the exponentials are computable.
)
)) Macro written 991119 by C. Bingham, kb@stat.umn.edu
)) Version of 000826 uses tabs() instead of covar
#$S(groups,y)
@groups <- argvalue($1,"group index","positive integer vector")
@y <- argvalue($2, "data", "nonmissing matrix")
@p <- ncols(@y)
@g <- max(@groups)
@N <- nrows(@y)

if (length(@groups) != @N){
	error("length of group index differs from number of cases")
}

@stats <- tabs(@y, @groups,count:T,mean:T,covar:T)
@n <- vector(@stats$count[,1]) # group sizes
@means <- @stats$mean # row i is mean of group i
@covars <- @stats$covar # @covars[i,,] is s for group i
delete(@stats)

if (min(@n) <= @p){
	error("All groups must have at least p+1 cases")
}

@varlabs <- if (haslabels(@y)){
	getlabels(@y,2)
}else{
	"@"
}
delete(@y)
@facname <- "$1"
if (!isname(@facname) || match("@*",@facname,0,exact:F) != 0){
	@facname <- "Group"
}
@facname <- paste(@facname," ",sep:"")

@grandmean <- vector(sum(@n*@means)/@N,labels:@varlabs)

 #make structures and row vector to fill
@l <- @q <- strconcat(split(rep(0,@g)'),labels:@facname)
@addcon <- matrix(rep(0,@g)',labels:structure("@",@facname))
for(@i,run(@g)){
	@s <- matrix(@covars[@i,,],labels:structure(@varlabs,@varlabs))
	@offset <- vector(@means[@i,] -	@grandmean');#ll(@s,@ll,@q,@addcon,@offset)
	@q[@i] <- -.5*solve(@s)
	@l[@i] <- vector(solve(@s,@offset),labels:@varlabs)
	@addcon[@i] <- 	-.5*@l[@i] %c% @offset - .5*log(det(@s))
}

delete(@groups,@p,@g,@i,@N,@varlabs,@facname,@means,@n)
structure(Q:delete(@q,return:T),L:delete(@l,return:T),\
	addcon:delete(@addcon,return:T),grandmean:delete(@grandmean,return:T))
%discrimquad%

probsquad macro dollars
) Macro to compute posterior probabilities based on a quadratic
) discriminant functions with coefficients computed by discrimquad()
)
) Usage:
)  probs <- probsquad(x,Str [, prior:P])
)
)   x            REAL m by p matrix or length p vector with no missing
)                values
)   Str          structure(Q:q,L:l,addcon:c,grandmean:ybar) as returned by
)                macro discrimquad().
)    q           structure of g square matrices
)    l           structure of g vectors of length p
)    c           1 by g row vector
)    ybar        length p vector
)    P           length g positive REAL vector proportional to
)                prior probabilities; default is rep(1,g)
)   probs        REAL m by g matrix of posterior probabilities for
)                rows of x or length g vector of posterior
)                probabilities for vector x
)   If x is a vector and p = length(grandmean) > 1, x must have length p
)   Otherwise x must be a matrix with p columns and probabilities are
)   for each row of x.
)
)) 991119, written by, C. Bingham, kb@stat.umn.edu
)) 011121 modified method of computing probabilities to avoid overflow
))        problem in computing exp(...)
#$S(x,structure(Q,L,addcon,grandmean) [,prior:p])
@x <- argvalue($1,"data","real matrix nonmissing")
@stuff <- argvalue($2,"quadratic discriminant info","structure")
if (!isstruc(@stuff[1]) || !isstruc(@stuff[2]) ||\
	!isvector(@stuff[3]',real:T) ||\
	!isvector(@stuff[4],real:T)){
	error("Argument 2 to $S has wrong form")
}
@p <- nrows(@stuff[1][1]) # dimension
@g <- ncomps(@stuff[1]) # number of groups
@P <- keyvalue($K,"prior","positive vector",default:rep(1,@g))

if (!isvector(@x) && ncols(@x) != @p ||\
	isvector(@x) && @p > 1 && length(@x) != @p){
	error("Wrong number of variables in data")
}

if(length(@P) != @g){
	error("Wrong length prior probability vector")
}

@n <- if(!isvector(@x) || @p == 1){
	nrows(@x)
}else{
	@x <- vector(@x)'
	1
}

@x <- @x' - @stuff$grandmean #transpose x and center it
 # The following is arcane, but it takes advantage of MacAnova's ability
 # to do matrix arithmetic with structures of matrices
 # @stuff$Q and @stuff$L are both structures with g components
 # @stuff$addcon is a 1 by g matrix (row vector)
@probs <- (matrix(vector(sum(@x * (@stuff$Q %*% @x + @stuff$L))),@n) +\
	 @stuff$addcon)'
if (haslabels(@stuff$Q)){
	@probs <- matrix(@probs,labels:structure(getlabels(@stuff$Q),"@"))
}
delete(@stuff,@x,@p,@g,@n)
@probs <- @P*exp(@probs - max(@probs))
@probs <- (@probs/sum(@probs))'
delete(@P)
delete(@probs,return:T)
%probsquad%

===> hotellval <===
hotellval MACRO DOLLARS
)Macro to compute a one-sample Hotelling's T^2 statistic
)Usage:
) t_sq <- hotellval(x)
)  x             n by p REAL matrix with no MISSING elements
)  t_sq          REAL scalar Hotellings T^2 for testing H_0:mu = 0
)
) result <- hotellval(x, pval:T)
)  result        structure(hotelling:t_sq, pvalue:P), where P =
)                P value computed from F-distribution
)) 030114. written by C. Bingham (kb@stat.umn.edu)
)) 030926 P-value calculation corrected
# $S(x [, pval:T])
@x <- matrix(argvalue($1,"x","real matrix nonmissing"))
@pval <- keyvalue($K,"pval*","TF",default:F)
@n <- nrows(@x)
@p <- ncols(@x)
if (@n <= @p){
	error("nrows(x) <= ncols(x)")
}
@xbar <- sum(@x)/@n
@x <-- @xbar
@result <- @n * (@n - 1) * @xbar %*% solve(@x %c% @x,@xbar')
if (@pval){
	@fe <- @n - 1
	@fe1 <- @fe - @p + 1
	@P <- 1 - cumF(@fe1*@result/(@fe*@p),@p,@fe1)
	@result <- strconcat(hotelling:@result, pvalue:delete(@P,return:T))
	delete(@fe, @fe1)
}
delete(@x,@n,@p,@xbar,@pval)
delete(@result,return:T)
%hotellval%

===> hotell2val <===
hotell2val MACRO DOLLARS
)Macro to compute a one-sample Hotelling's T^2 statistic
)Usage:
) t_sq <- hotell2val(x1, x2)
)  x1            n1 by p REAL matrix with no MISSING elements
)  x2            n2 by p REAL matrix with no MISSING elements
)  t_sq          REAL scalar Hotellings T^2 for testing
)                H_0: mu1 = mu2
) result <- hotellval(x1, x2, pval:T)
)  result        structure(hotelling:t_sq, pvalue:P), where P =
)                P value computed from F-distribution
)) 030114 Written by C. Bingham (kb@stat.umn.edu)
)) 030926 P-value calculation corrected
# $S(x1, x2 [, pval:T])
@x1 <- argvalue($1,"x1","real matrix nonmissing")
@x2 <- argvalue($2,"x2","real matrix nonmissing")
@pval <- keyvalue($K,"pval*","TF",default:F)
@n1 <- nrows(@x1)
@n2 <- nrows(@x2)
@p <- ncols(@x1)
if (ncols(@x2) != @p){
	error("ncols(x1) != ncols(x2)")
}
@fe <- @n1 + @n2 - 2
if (@fe <= @p - 1){
	error("n1 + n2 < ncols(x1) + 1")
}
@xbar1 <- describe(@x1,mean:T)
@xbar2 <- describe(@x2,mean:T)
@d <- @xbar1 - @xbar2
@x1 <-- @xbar1'
@x2 <-- @xbar2'
@vhat <- (1/@n1 + 1/@n2)*(@x1 %c% @x1 + @x2 %c% @x2)/@fe
delete(@x1,@x2,@xbar1,@xbar2,@n1,@n2)
@result <- @d %c% solve(@vhat, @d)
if (@pval){
	@P <- 1 - cumF((@fe - @p + 1)*@result/(@fe*@p),@p,@fe-@p+1)
	@result <- strconcat(hotelling:@result,pvalue:delete(@P,return:T))
}
delete(@pval,@p,@vhat,@d)
delete(@result,return:T)
%hotell2val%


===> dastepsetup <===
dastepsetup        MACRO  DOLLARS
) Setup stepwise discriminant
) Usage:
)   dastepsetup(Model [,allin:T or in:ins] [,silent:T])
)     Model            glm model, usually of the form "y = groups"
)     ins              LOGICAL vector with length = ncols(y) or
)                      vector of unique positive integers <= ncols(y)
)
)   dastepsetup([allin:T or in:ins] [,silent:T]), with no model,
)    is equivalent to 
)   dastepsetup(STRMODEL [,allin:T or in:ins] [,silent:T])
)
) dastepsetup() initializes an "invisible" structure _DASTEPSTATE and,
) except with silent:T, prints F-to-enter with P-values
)
) With allin:T, stepsetup() starts with all dependent variables "in".
)
) When ins is an integer vector, stepsetup() starts with the indicated
) variables "in".  in:run(ncols(y)) has same effect as allin:T
)
) When ins is a LOGICAL vector, in:ins has the same effect as 
) in:run(ncols(y))[ins], that is T indicates a variable is "in"
) in:rep(T,ncols(y)) has the same effect as allin:T
)
)) 000530 Written by C. Bingham, kb@umn.edu
)) 031226 pushmodel() and popmodel() invoked to preserve current
))        model information and variables
)) 041206 Value of 'in' can be vector of unique integers between 1 and p
#$S(model [,silent:T] [,in:logVector  or  allin:T])
@model <- if ($v > 0){
	$01
}else{
	NULL
}
@model <- if (!isnull(@model)){
	argvalue(@model,"model","string")
}elseif(isscalar(STRMODEL,char:T)){
	STRMODEL
}else{
	error("No model provided and STRMODEL does not exist")
}
if(match("*=*",@model, 0, exact:F) == 0){
	error(paste("Illegal model \"",@model,"\"",sep:""))
}

@allin <- keyvalue($K,"allin","TF",default:F)
@instart <- keyvalue($K,"in","vector nonmissing")
if (@allin && !isnull(@instart)){
	error("keyword 'in' illegal with 'allin:T'")
}

@nocon <- match("*-1",@model,0, exact:F) != 0 ||\
	match("*- 1",@model,0, exact:F) != 0

@pushed <- isfunction(pushmodel)
if (@pushed) {
	@pushed <- pushmodel(canpush:T)
	if (@pushed){
		pushmodel()
	}
}
manova(@model,silent:T) #compute side effect variables
@nvars <- dim(SS)[2]

if (!isnull(@instart)){
	if (islogic(@instart)){
		if(length(@instart) != @nvars){
			error("Length of LOGICAL 'instart' != number of variables")
		}
	} else {
		if (!isvector(@instart,positive:T,integer:T) || \
			length(unique(@instart)) != length(@instart) || \
			max(@instart) > @nvars) {
			error("Value of 'in' not vector of positive integers <= ncols(y)")
		}
		@instart <- match(run(@nvars),@instart,0) != 0
	}
	@allin <- (sum(@instart) == @nvars)
}else{
	@instart <- rep(@allin,@nvars)
}

@nterms <- length(DF) - 1

@fe <- DF[@nterms+1]
@fh <- sum(DF[run(2,@nterms)])

@e <- matrix(SS[@nterms+1,,])
@hpluse <- matrix(sum(SS[-1,,]))
if (@nocon){
	@fh <-+ DF[1]
	@hpluse <-+ matrix(SS[1,,])
}

if(haslabels(SS)){
	@varnames <- getlabels(SS,2)
	if(!isname(@varnames[1])){
		delete(@varnames)
	}
}
if (delete(@pushed,return:T)){
	popmodel()
}
if (!isvector(@varnames)){
	@varnames <- getlabels(vector(rep(0,@nvars),labels:"X"))
}

@Fs <- matrix(rep(0, 2*@nvars), 2, labels:structure(vector("F", "P"),\
	@varnames))

if (!@allin){
	if (sum(@instart) > 0){
		@hpluse <- swp(@hpluse, run(@nvars)[@instart])
		@e <- swp(@e, run(@nvars)[@instart])
		@fe <-- sum(@instart)
	}
}else{
	@hpluse <- solve(@hpluse)
	@e <- solve(@e)
	@fe <-- @nvars
}

@history <- if(sum(@instart) > 0){
	run(@nvars)[@instart]
}else{
	NULL
}

setlabels(@hpluse, structure(@varnames,@varnames))
setlabels(@e, structure(@varnames,@varnames))
@fe <- vector(@fe,@fe+1,labels:vector("enter", "remove"))

_DASTEPSTATE <- structure(model:delete(@model,return:T),\
	hpluse:delete(@hpluse,return:T),\
	e:delete(@e,return:T),\
	in:vector(delete(@instart,return:T),labels:@varnames),\
	F:delete(@Fs,return:T), fh:delete(@fh,return:T), fe:delete(@fe,return:T),\
	history:delete(@history,return:T) )

delete(@nocon,@nvars,@varnames,@allin)

if(!ismacro(dastepstatus)){
	getmacros(dastepstatus,quiet:T,printname:F)
}
dastepstatus($K)
%dastepsetup%

===> dasteplook <===
dasteplook  MACRO DOLLARS
) Macro to acess elements of _DASTEPSTATE
) Usage:
)  answers <- dasteplook(name_1, name_2, ...)
)   name_k             one of 'model', 'hpluse', 'e', 'in', 'F',
)                      'fh', 'fe', 'history' and 'varnames'
)   answers            scalar, vector or matrix (with one name)
)                      or k-component structure (with k names)
)  answers <- dastelook(all)
)   answers            structure identical to _DASTEPSTATE; does
)                      not contain variable names
)
) Argument 'varnames' requests a vector of variable names, extracted as
) the labels of _DASTEPSTATE$in.
)
) Except for 'varnames', the names recognized are names of _DASTEPSTATE
) components
)
)) Version 000826
)) 031224 exact match of names is no longer required; use enough to specify
))        the name uniquely.  Also, the variable names are returned when
))        'varnames' is an argument
# $(name1, name2, ...) names model, hpluse, e, in, F, fh, fe, history, varnames
if ($k > 0){
	error("no keywords allowed as arguments")
}
@names <- $A
for(@i,run($v)){
	# strip off double quotes, if any
	if (match("\"*\"",@names[@i], 0,exact:F) != 0){
		@names[@i] <- <<@names[@i]>>
	}
}
@all <- match("all",@names,0) != 0

if(delete(@all,return:T)){
	if ($v > 1){
		error("No other arguments are allowed with 'all'")
	}
	delete(@names,@i)
	@result <- _DASTEPSTATE
} else {
	@compnames <- compnames(_DASTEPSTATE)
	@ncomps <- length(@names)
	@result <- split(rep(0,@ncomps)')
	for(@i,1,@ncomps){
		if (match("var*",@names[@i],0,exact:F) > 0) {
			@names[@i] <- "varnames"
			@result[@i] <- getlabels(_DASTEPSTATE$in)
		} else {
			if (@names[@i] == "h" || @names[@i] == "f") {
				error(paste("\"",@names[@i],"\" doesn't fully specify name",sep:""))
			}
				
			@k <- match(paste(@names[@i],"*",sep:""),@compnames,0,exact:F)
			if(@k == 0){
				error(paste("\"",@names[@i],"\" not a legal name",sep:""))
			}
			@names[@i] <- @compnames[@k]
			@result[@i] <- _DASTEPSTATE[@k]
			delete(@k)
		}
	}
	@result <- if (delete(@ncomps,return:T) == 1){
		@result[1]
	} else {
		strconcat(@result,compnames:@names)
	}
	delete(@compnames,@names,@i)
}
delete(@result,return:T)
%dasteplook%

===> dastepstatus <===
dastepstatus       MACRO DOLLARS
) Macro to update the F-statistics and P-values in _DASTEPSTATE
) Usage:
)  dastepstatus()
)   Updates the F-statistics and P-values in _DASTEPSTATE and prints
)   a summary of the updated information
)  dastepstatus(silent:T)
)   Updates _DASTEPSTATE but prints nothing
)
) In both usages, dastepstatus() returns as value an updated copy of
) _DASTEPSTATE as an invisible variable
)) Version 000530
# $S([silent:T])
@silent <- keyvalue($K,"silent","TF",default:F)
@hpluse <- _DASTEPSTATE$hpluse
@e <- _DASTEPSTATE$e
@nvars <- nrows(@e)
@Fs <- _DASTEPSTATE$F
@in <- _DASTEPSTATE$in
@out <- !@in
@u <- diag(@hpluse)
@v <- diag(@e)
@fh <- _DASTEPSTATE$fh
@fe <- _DASTEPSTATE$fe
if (sum(@out) > 0){
	@Fs[1,@out] <- (@fe[1]/@fh)*(@u[@out]/@v[@out] - 1)
	@Fs[2,@out] <- 1 - cumF(@Fs[1,@out],@fh,@fe[1])
}
if (sum(@in) > 0){
	@Fs[1,@in] <- ((@fe[2])/@fh)*(@v[@in]/@u[@in] - 1)
	@Fs[2,@in] <- 1 - cumF(@Fs[1,@in],@fh,@fe[2])
}
@ii <- match("F",compnames(_DASTEPSTATE))
_DASTEPSTATE[@ii] <- delete(@Fs,return:T)
delete(@ii)

if (!delete(@silent,return:T)){
	@nin <- sum(@in)
	print(paste("Model: \"",_DASTEPSTATE$model,"\"",sep:""))
	if(@nin > 0) {
		print(paste("F(",@fh, ",",@fe[2],") to delete",sep:""))
		print(_DASTEPSTATE$F[, @in], header:F)
	}else{
		print("No variables are \"in\"")
	}
	print("")
	if(@nin < @nvars) {
		print(paste("F(",@fh,",",@fe[1],") to enter",sep:""))
		print(_DASTEPSTATE$F[, !@in], header:F)
	}else{
		print("All variables are \"in\"")
	}
	print("")
	delete(@nin)
}
delete(@fe,@fh,@in,@nvars)
_DASTEPSTATE # invisible so it won't print but can be assigned
%dastepstatus%

===> daentervar <===
daentervar     MACRO DOLLARS
) Macro to take one or more steps forward, entering one or more
) response variables into the current stepwise discriminant
) analysis
) Usage:
)  daentervar(var1 [,var2 ...] [,silent:T])
)    var1, var2 ...     quoted or unquoted names or numbers of variables
)                       to be entered, or CHARACTER or integer scalars
)                       specifying the variables			
)
)  Without silent:T, a summary of the state of the stepwise process
)  is printed
)) Version 000530
)) 031224 Modified to use davarid() to decode arguments
)) 041011 Load davarid() silently if not present0
# $S(var1 [,var2 ...] [,silent:T]), varj a variable name or number
if ($v == 0 || $k > 1){
	error("usage: $S(var1 [,var2 ...] [,silent:T])", macroname:F)
}
if (!isstruc(_DASTEPSTATE)){
	error("You must run dastepsetup(Model) before using $S()",macroname:F)
}
@varnames <- getlabels(_DASTEPSTATE$in)
@nvars <- length(@varnames)
@vars <- $A
if(!ismacro(dastepstatus)){
	getmacros(dastepstatus,quiet:T,printname:F)
}
if(!ismacro(davarid)){
	getmacros(davarid,quiet:T,printname:F)
}
for (@j, 1, length(@vars)){
	@var <- @vars[@j]
	if(match("*:*",@var,0,exact:F) != 0){
		next
	}

	# not a keyword phrase
	@newvar <- davarid(@var)
	if (@newvar[1] == 0) {
		error(@newvar[2])
	}
	
	if (_DASTEPSTATE$in[@newvar[1]]){
		error(paste("Variable \"",@newvar[2],"\" already 'in'",sep:""))
	}
	@newvar <- @newvar[1]
	delete(@var)

	@ii <- match("hpluse",compnames(_DASTEPSTATE))
	_DASTEPSTATE[@ii] <- swp(_DASTEPSTATE$hpluse, @newvar)
	@ii <- match("e",compnames(_DASTEPSTATE))
	_DASTEPSTATE[@ii] <- swp(_DASTEPSTATE$e, @newvar)
	@in <- _DASTEPSTATE$in
	@in[@newvar] <- T
	@ii <- match("in",compnames(_DASTEPSTATE))
	_DASTEPSTATE[@ii] <- @in
	@ii <- match("fe",compnames(_DASTEPSTATE))
	_DASTEPSTATE[@ii] <- _DASTEPSTATE$fe - 1
	@ii <- match("history",compnames(_DASTEPSTATE))
	_DASTEPSTATE[@ii] <- vector(_DASTEPSTATE$history,@newvar)

	delete(@ii,@newvar,@in)
	dastepstatus($K)
} #for (@j, 1, length(@vars))
delete(@vars, @j, @varnames)
_DASTEPSTATE
%daentervar%

davarid    MACRO  DOLLARS
) Utility macro used by daentervar() and daremovevar() to identify
) a variable (column) in a data matrix.  It extracts information from
) variable _DASTEPSTATE which must be defined
) Usage:
)  Id <- davarid(id)
)    id         CHARACTER scalar or quoted string identifying
)               a column of the response matrix in stepwise
)               variable selection.
)    Id         structure(varnumber:k, varname:Name) where k is
)               an integer and Name is a CHARACTER scalar, the
)               k-th label of IN = _DASTEPSTATE$in
) Forms for id:
)   "name"      "name" matches the label j of IN
)               Id$varnumber = j, Id$varname = "name"
)   "3"         Id$varname = 3, Id$varname = label 3 of IN
)   "J"         J an integer scalarvariable ; Id$varnumber = J,
)               Id$varname = label J of IN
)   "x3"        Id$varname = 3, Id$varname = label 3 of IN; if
)               "x3" matches a label of IN, the label number takes
)               precedence.  In place of "x", any name not starting
)               with '@' can be used.
)) 031227 written by C. Bingham, kb@umn.edu
#$S(id), id a CHARACTER scalar
@var <- argvalue($1,"arg 1", "string")
if (match("\"*\"", @var, 0, exact:F) != 0){
	@var <- <<@var>> #strip off quotes if present
}

if (!isstruc(_DASTEPSTATE)){
	error("You must run dastepsetup(Model) before using $S()",macroname:F)
}
@nvars <- nrows(_DASTEPSTATE$in)
@varnames <- getlabels(_DASTEPSTATE$in)

@outvar <- match(@var, @varnames, 0)
if (@outvar == 0) { # not a direct match
	if (isnumber(@var)) {
		@outvar <- <<@var>>
		if (!isscalar(@outvar,integer:T)) {
			@outvar <- -1000
		} elseif (@outvar <= 0) {
			@outvar <-- 1
		}
	} else {
		if (isname(@var)) {
			if (match("@*",@var,0,exact:F) == 0 && \
				max(match(vector("*0","*1","*2","*3","*4","*5","*6","*7","*8","*9"),\
						  @var,0,exact:F)) > 0) {
				# arg may have form "v3",  "x1" or something similar
				@outvar <- vecread(string:@var,silent:T,byfields:F)
				if (!isscalar(@outvar,integer:T)){
					@outvar <- 0
				}
			} elseif (isscalar(<<@var>>,integer:T)) {
				@outvar <- <<@var>>
				if (@outvar <= 0) {
					@outvar <-- 1
				}
			}
		}
		
		if (alltrue(@outvar == 0, isscalar(<<@var>>, character:T))) {
			@var <- <<@var>>
			@outvar <- match(@var, @varnames, 0)
			if (alltrue(@outvar == 0,\
				max(match(vector("*0","*1","*2","*3","*4","*5","*6","*7",\
								 "*8","*9"),\
					  @var,0,exact:F)) > 0)) {
				@outvar <- vecread(string:@var,silent:T,byfields:F)
				if (!isscalar(@outvar,integer:T)){
					@outvar <- 0
				}
			}
		}
	}
}
@var <- if (@outvar == 0) {
	paste($1,"does not specify a variable in the model")
} elseif (@outvar == -1000) {
	@outvar <- 0
	paste($1,"does not specify integer")
} elseif (@outvar < 0) {
	@outvar <- 0
	"Variable number <= 0"
} elseif (@outvar > @nvars) {
	@outvar <- 0
	"Variable number > number of variables"
} else {
	@varnames[@outvar]
}
delete(@nvars,@varnames)
structure(varnumber:delete(@outvar,return:T),varname:delete(@var,return:T))
%davarid%

===> daremovevar <===
daremovevar    MACRO  DOLLARS
) Macro to take one or more steps backward, removing a variable or
) several variables from the current stepwise disciminant analysis
) Usage:
)  daremovevar(var1 [,var2 ...] [,silent:T])
)   var1, var2 ...  names or numbers of variables to be removed
)                   or scalar variables containing names or numbers
)   silent:T        Suppress output
)
) Adapted from a macro by Gary Oehlert
)) Version 000530
)) 031219 modification so that a varj can be a integer or character
)) scalar variable
# $S(var1 [,var2 ...] [,silent:T]), varj a variable name or number
if ($v == 0 || $k > 1){
	error("usage is $S(var1 [,var2 ...] [,silent:T])")
}
if (!isstruc(_DASTEPSTATE)){
	error("You must run dastepsetup(Model) before using $S()",macroname:F)
}
@nvars <- nrows(_DASTEPSTATE$e)
@varnames <- getlabels(_DASTEPSTATE$in)
@vars <- $A
if(!ismacro(dastepstatus)){
	getmacros(dastepstatus,quiet:T,printname:F)
}
if(!ismacro(davarid)){
	getmacros(davarid,quiet:T,printname:F)
}
for (@j, 1, length(@vars)){
	@var <- @vars[@j]
	if(match("*:*",@var,0,exact:F) != 0){
		next
	}
	# not a keyword phrase
	@outvar <- davarid(@var)
	if (@outvar[1] == 0) {
		error(@outvar[2])
	}
	delete(@var)
	
	if (!_DASTEPSTATE$in[@outvar[1]]){
		error(paste("Variable \"", @outvar[2],"\" not currently 'in'",sep:""))
	}
	@outvar <- @outvar[1]

	@ii <- match("hpluse",compnames(_DASTEPSTATE))
	_DASTEPSTATE[@ii] <- swp(_DASTEPSTATE$hpluse, @outvar)
	@ii <- match("e",compnames(_DASTEPSTATE))
	_DASTEPSTATE[@ii] <- swp(_DASTEPSTATE$e, @outvar)
	@in <- _DASTEPSTATE$in
	@in[@outvar] <- F
	@ii <- match("in",compnames(_DASTEPSTATE))
	_DASTEPSTATE[@ii] <- @in
	@ii <- match("fe",compnames(_DASTEPSTATE))
	_DASTEPSTATE[@ii] <- _DASTEPSTATE$fe + 1
	@ii <- match("history",compnames(_DASTEPSTATE))
	_DASTEPSTATE[@ii] <- vector(_DASTEPSTATE$history,-@outvar)
	dastepstatus($K)
}
delete(@ii,@outvar, @varnames, @vars, @j, @nvars)
_DASTEPSTATE
%daremovevar%

===> dastepselect <===
dastepselect  macro  dollars
) Macro to perform automatic forward or backward stepwise variable select
) for discriminant analysis
) Usage:
)  ins <- dastepselect(Model, alpha [, allin:T] [,updown:T] [,bonferroni:F]\
)               [,verbose:T or silent:T])
)    Model        quoted string or CHARACTER scalar specifying a
)                 manova() model, usually of the form "y = a", where y
)                 is a REAL matrix and a is a factor
)    alpha        positive REAL scalar 1 < alpha < 50 or 0 < alpha < .5,
)                 significance level (in percent when alpha > 1)
)    allin:T      Start with all variables in and step backwards;
)                 otherwise, start with all variables out and step
)                 forwards
)    updown:T     Before entering a variable, check to see whether one
)                 can be removed.
)    bonferroni:F Do not Bonferronize P-values
)    quiet:T      Suppress printing of action taken at each step and
)                 summary and make value returned "visible"
)    silent:T     Same as 'quiet:T' except warning messages are also
)                 suppressed;  incompatible with 'quiet:F'
)    ins          Integer vector of variable numbers selected or
)                 NULL if no variable selected; "invisible" unless
)                 quiet:T
)
) A decision to enter or remove a variable is based on a P-value for an
) F-to-enter (forward step) or F-to-remove (backward).
)
) Unless bonferroni:F is an argument, P-values are Bonferronized by
) by multiplication by an integer K.
)
) On a forward step, K = number of "outs".
) On a backward step, K = (number of "outs") + 1.
)
) This choice ensures a variable that has just entered will not immediately
) be a candidate for removal, and a variable that has just been removed
) will not immediately be a candidate to enter.
)
)) 031226 Written by C. Bingham, kb@umn.edu
)) 041011 Warning message when no variables selected unless suppressed
))        by new keyword phrase silent:T; 'quiet' replaces 'verbose'
)) 041206 Added keyword 'in' to specify a starting position
)) 041215 Corrected bug.  It now uses F-statistic to choose the candidate
))        for entering or removal rather than P-values.  The problem was
))        having multiple P's = 0, only one of which went with the largest F
#$S(Model, alpha [, allin:T or in:logvec] [,updown:T] [,verbose:T or silent:T] [,bonferroni:F])
@model <- $01
@model <- if (!isnull(@model)){
	argvalue(@model,"model","string")
}elseif(isscalar(STRMODEL,char:T)){
	STRMODEL
}else{
	error("No model provided and STRMODEL does not exist")
}
if(match("*=*",@model, 0, exact:F) == 0){
	error(paste("Illegal model \"",@model,"\"",sep:""))
}

@alpha <- argvalue($2,"alpha","positive number")
if (@alpha >= 50) {
	error("alpha >= 50%")
}
if (@alpha > 1) { # treat alpha between 1 and 50 as percent
	@alpha <-/ 100
} elseif (@alpha >= .5) {
	error("alpha >= .5")
}
@allin <- keyvalue($K, "allin", "TF",default:F)
@updown <- keyvalue($K, "updown","TF",default:F)
@silent <- keyvalue($K, "silent","TF",default:F)
@quiet <- keyvalue($K, "quiet","TF",default:@silent)
@instart <- keyvalue($K,"in","logical vector nonmissing")
if (@allin && !isnull(@instart)){
	error("keyword 'in' illegal with 'allin:T'")
}

if (!@quiet && @silent){
	error("'silent:T' can't be used with 'quiet:F'")
}
@bonfer <- keyvalue($K,"bonfer*","TF",default:T)

if (!ismacro(dastepsetup)) {
	getmacros(dastepsetup,silent:T)
}
if (!ismacro(dasteplook)) {
	getmacros(dasteplook, silent:T)
}
if (@allin || @updown) {
	if (!ismacro(daremovevar)) {
		getmacros(daremovevar, silent:T)
	}
}
if (!@allin || @updown) {
	if (!ismacro(daentervar)) {
		getmacros(daentervar, silent:T)
	}
}

if (!isnull(@instart)){
	dastepsetup(@model,in:@instart,silent:T)
} else {
	dastepsetup(@model,allin:@allin,silent:T)
}
@in <- dasteplook(in)
@names <- getlabels(@in)
@what <- if (@bonfer) {
	"Bonferronized P-value-to-"
} else {
	"P-value-to-"
}

@direction <- if (@allin) {
	-1
} else {
	1
}

while (@direction != 0) {
	@out <- !@in
	if (@direction  < 0) { # back stepping
		if (sum(@in) == 0) {
			# all out, so we are done
			@direction <- 0
		} else {
			@f_info <- dasteplook(F)
			@i <- grade(@f_info[1,@in]',up:T)[1]
			# see if we should delete a variable
			@pval <- @f_info[2,@in][@i]
			if (@bonfer) {
				@pval <-* sum(@out) + 1
			}
			@direction <- if (@pval <= @alpha) {
				if (!@updown || @direction < -1) {
					# can't delete, so we are done
					0
				} else {
					# Try forward next time
					2 # indicates failure to delete
				}
			} else {
				# delete variable and continue backwards
				@var <- @names[@in][@i]
				if (!@quiet) {
					print(paste("Removing ",@var,". ", @what, "remove = ",\
								@pval,sep:""))
				}
				daremovevar(@var, silent:T)
				-1 # indicates success at deleting
			}
		}
	} else { # forward stepping
		if (sum(@out) == 0) {
			# all are in so stop
			@direction <- 0
		} else {
			@f_info <- dasteplook(F)
			@i <- grade(@f_info[1,@out]',down:T)[1]
			@pval <- @f_info[2,@out][@i]
			if (@bonfer) {
				@pval <-* sum(@out)
			}
			@direction <- if (@pval > @alpha) {
				if (!@updown || @direction > 1) {
					# can't enter so we are done
					0
				} else {
					# try forward next time
					-2 # indicates failure to enter
				}
			} else {
				@var <- @names[@out][@i]
				if (!@quiet) {
					print(paste("Entering ",@var,". ",@what,"enter = ",@pval,sep:""))
				}
				daentervar(@var, silent:T)
				if (!@updown) {
					1
				} else {
					# Next time try to delete
					-1
				}
			}
		}
	}
	@in <- dasteplook(in)
}
delete(@f_info,@i,@pval,@var,silent:T)	
@out <- !@in
if (!@quiet) {
	if ((@allin || @updown) && sum(@in) > 0) {
		@pval <- max(vector(dasteplook(F)[2,@in]))
		if (@bonfer) {
			@pval <-* sum(@out) + 1
		}
		print(paste("Largest ",@what,"remove = ",@pval," < ",@alpha,sep:""))
		delete(@pval)
	}
	if ((!@allin || @updown) && sum(@out) > 0) {
		@pval <- min(vector(dasteplook(F)[2,@out]))
		if (@bonfer) {
			@pval <-* sum(@out)
		}
		print(paste("Smallest ",@what,"enter = ",@pval," > ",@alpha,sep:""))
		delete(@pval)
	}
}

@names <- if (sum(@in) > 0) {
	@names[@in]
} else {
	NULL
}

if (!@quiet) {
	if (!isnull(@names)) {
		print(paste("Variables selected:",paste(@names, sep:", ")))
	} else {
		print("No variables selected")
	}
}

@result <- if (!isnull(@names)) {
	vector(run(length(@in))[@in],labels:@names)
} else {
	if (!@silent){
		print("WARNING: no variables selected; result is NULL")
	}
	NULL
}

delete(@in,@out,@names,@bonfer,@alpha,@what,@silent)
delete(@result, return:T, invisible:!delete(@quiet,return:T))
%dastepselect%

===> dascreen <===
dascreen       macro   dollars
) Macro to examine all possible subsets for best linear
) discrimination
) Usage:
)  result <- dascreen(model [,mbest:m] [, method:Method] [,penalty:pen]\
)                       [, forced:varlist] [,maxvar:maxv] [,silent:T])
)   model     CHARACTER scalar of the form "x=a", where a is a factor
)             x is a real matrix, and nrows(x) = nrows(a)
)   m         positive integer <= 2^p-1, default = 5
)   pen       REAL scalar >= 0, used as multiplier of number of
)             parameters in computing AIC; default is pen = 2
)   varlist   vector of distinct integers > 0, a list of variables
)             to be included in every subset considered
)   maxv      positive integer, default:15, specifying the maximum
)             number of variables allowed, not counting those
)             specified by 'forced:varlist'.  If it is exceeded a
)             warning is printed, suggesting a value for maxv.
)   Method    "aic" or "penaic", default is "aic"
)   silent:T  Suppress any warning messages
)   result    structure(subsets:S, criterion:Crit)
)    S        length p1 vector or a p1 by m matrix, with leading
)             elements of each column the variable numbers retained
)             and p1 = size of largest set
)    Crit     m vector of values of the criterion for each subset
)
) This actually works with any model with only one error term
) If CONSTANT is the first term it is not considered part of the model
) in computing the criterion
)
) dascreen() works by brute force, examining all 2^p1 - 1 non empty
) subsets of variables on the LHS of the model, where p1 = p - number of
) variables forced in.  The amount of time goes up by about a factor of
) 2 for every additional response variable
)
) The default criterion is AIC =
)   -N Log(det(H + E)/det(E) + pen*nparams, nparams = q*fh + q*(q+1)/2,
)   where q = size of subset and pen = 2 or value of keyword 'penalty'
)
) See Yasunori Fujikoshi, Selection of variables for discriminant
) analysis in a high-dimensional case, Sankhya Series A, 64, (2002)
) 256-267
)
))031127 written by C. Bingham, kb@umn.edu
))031227 added keywords maxvar, penalty, forced and silent
))041011 When mbest > 1, subset matrix is mbest by p1
))041012 Changed default for mbest to 5, made sets correspond to
))       columns of S rather than rows; added labels to S and Crit
))041220 Corrected computation of number of parameters following
))       Fujikoshi 2002
))       Criterion for empty set is now computed when no variables
))       are forced in
# $S(model [,mbest:m] [,method:Method] [,penalty:pen] \
#    [,forced:varlist] [,silent:T])
@model <- argvalue($1, "model", "string")

@keys <- if ($k > 0){
	structure($K,notakey:NULL)
} else {
	structure(notakey:NULL)
}
@mbest <- keyvalue(@keys,"mbest","positive count",default:5)
@penalty <- keyvalue(@keys,"pen*","nonneg number",default:2)
@maxvar <- keyvalue(@keys,"max*","positive count",default:15)
@method <- keyvalue(@keys,"meth*","string",default:"aic")
@forced <- keyvalue(@keys,"forc*","postive integer vector")
@silent <- keyvalue(@keys,"silent","TF",default:F)
delete(@keys)
@Methods <- vector("aic","penaic")
@jmeth <- match(paste(@method,"*",sep:""),@Methods,0,exact:F)

if (@jmeth == 0){
	error(paste("Unrecognized screening method\"",@method,"\"",sep:""))
}
@method <- delete(@Methods,return:T)[@jmeth]

if (match("*=*",@model,0,exact:F) == 0) {
	error(paste("\"",@model,"\" is not a legal model",sep:""))
}

@nforced <- if (isnull(@forced)){
	0
} else {
	length(@forced)
}
if (alltrue(@nforced > 0, length(unique(@forced)) != @nforced)){
	error("Duplicate values in value of 'forced'")
}
if (alltrue(isfunction(pushmodel),pushmodel(canpush:T))) {
	pushmodel()
}

manova(@model,silent:T)

@p <- ncols(RESIDUALS)
@N <- nrows(RESIDUALS)
if (alltrue(@nforced > 0, max(@forced) > @p)){
	error(paste("Maximum forced value > dimension =",@p))
}

if (@maxvar + @nforced < @p) {
	@result <- NULL
	if (!@silent){
		print(paste(sep:"","WARNING: # of non-forced variables > ", @maxvar,\
			"; use 'maxvar:",@p - @nforced," if you've got lots of time"))
	}
} else {
	@nterms <- length(DF)
	@jterms <- run(@nterms-1)
	if (TERMNAMES[1] == "CONSTANT") {
		if (@nterms == 2){
			error(paste("No non-constant terms in \"",@model,"\"",sep:""))
		}
		@jterms <- @jterms[-1]
	}
	@fh <- sum(DF[@jterms])
	@fe <- DF[@nterms]
	if (@fh == 0) {
		error(paste("No non-constant terms in \"",@model,"\"",sep:""))
	}
	if (@fe == 0) {
		error("No degrees of freedom for error")
	}
	@e <- matrix(SS[@nterms,,])
	@hpluse <- matrix(sum(SS[@jterms,,])) + @e
	delete(@jterms,@nterms)

	if (@nforced == @p) {
		if (!@silent){
			print("WARNING: all variables forced in; no selection done")
		}
		#if (@jmeth <= 2){
			@crit <- -@N*log(det(@hpluse)/det(@e)) +\
				@penalty*@p*(@p+1 + 2*@fh)/2
		#}
		@result <- structure(subsets:run(@p),criterion:@crit)
	} else {
		@p1 <- @p - @nforced
		@nsets <- 2^@p1
		if (@mbest > @nsets){
			if (!@silent){
				print(paste("WARNING: value for 'mbest' too big; reduced to",\
					@nsets))
			}
			@mbest <- @nsets
		}
		@eligible <- if (@nforced == 0){
			run(@p)
		} else {
			run(@p)[-@forced]
		}
		@crit <- rep(0, @nsets)
		@nparamsig <- @p*(@p+1)/2
		if (@nforced > 0){
			# the criterion for forced-in set
			@loglik <-\
			   -@N*log(det(@hpluse[@forced,@forced])/det(@e[@forced,@forced]))
			@nparam <- @nparamsig + @nforced*@fh + @p
		} else {
			# the criterion for the empty set
			@loglik <- 0
			@nparam <- @nparamsig + @p
		}
		@crit[@nsets] <- @loglik + @penalty*@nparam
		
		@powersof2 <- 2^run(0,@p1-1)
		for(@i,1,@nsets-1){
			@set <- vector(@eligible[(@i %& @powersof2) != 0],@forced)
			@q <- length(@set) # size of subset

			#if (@jmeth <= 2){
				@loglik <- \
				  -@N*log(det(@hpluse[@set,@set])/det(@e[@set,@set]))
				@nparam <- @nparamsig + @q*@fh + @p
				@crit[@i] <- @loglik + @penalty*@nparam
			#}
		}
		delete(@nparam,@nparamsig,@loglik)

		@jbest <- grade(@crit)[run(@mbest)]
		@crit <- @crit[@jbest]

		@maxp <- if (@nforced > 0 && @mbest == 1 && @jbest[1] == @nsets) {
			@nforced
		} else {
			max(nbits(@jbest)) + @nforced
		}
		@subsets <- matrix(rep(0,@mbest*@maxp),@maxp)
		for(@i,1,@mbest) {
			@set <- @jbest[@i]
			if (@set == @nsets) {
				if (@nforced > 0){
					@subsets[run(@nforced), @i] <- sort(@forced)
				}
			} else {
				@subsets[run(nbits(@set) + @nforced), @i] <- \
					sort(vector(@eligible[(@set %& @powersof2) != 0], @forced))
			}
		}
		if (@mbest == 1) {
			@subsets <- vector(@subsets)
		} else {
			setlabels(@subsets,structure("@","Set "))
			setlabels(@crit,"Set ")
		}
		delete(@i,@mbest,@set,@jbest,@p,@p1,@maxp,@powersof2, @nsets)
		@result <-structure(subsets:delete(@subsets,return:T),\
						  criterion:delete(@crit,return:T))
	}
	delete(@fh,@fe,@e,@hpluse)
}
if (alltrue(isfunction(popmodel),popmodel(canpop:T))) {
	popmodel()
}
delete(@result,return:T)
%dascreen%

===> mlcrit <===
mlcrit       MACRO DOLLARS
) Macro to compute ML factor analysis criterion function and its gradient
) The parameters x are log(psi_sub_i), i = 1,...,p
) This is to be used with optimization macros dfp(), bfs() and broyden()
) Usage
)  result <- mlcrit(logpsi, k, structure(s, nfact [, del]))
)   logpsi       REAL vector of log uniquenesses
)   k            integer scalar to control what is calculated
)                k = 0, return criterion (scalar) (MISSING when s not P.D.)
)                k = 1, return gradient (vector) (rep(MISSING,nrows(s))
)                  when s not P.D.)
)                k = -2, return structure(loadings, eigenvals)
)   s            REAL symmetric covariance or correlation matrix
)   nfact > 0    integer number of factors
)   del >= 0     Small REAL scalar [1e-5].  Gradient is computed as
)     gradient[i] = (critfn(logpsi + del*(run(p)==i)) - critfn(logpsi))/del
)                del = 0 results in analytic gradient being used
)   When s not positive definite, MISSING (k = 0) or
)   rep(MISSING,nrows(s)) (k = 1) is returned
) No argument checking is done
)) 011120 First draft
)) 011124 del = 0 signals use of analytic gradient
)) Version 011124, C. Bingham, kb@stat.umn.edu
@x <- $1
@n <- length(@x)
_PSIML <- exp(@x)

@param <- $3
@s <- @param[1]
@J <- run(@param[2])
@result <- if (($2) < 0) {
	if (($2) == -2) {
	  # return structure(loadings, eigenvals)
		@stuff <- releigen(@s, dmat(_PSIML))
		@gamma <- @stuff$values[@J]-1
		@crit <- NULL
		@loadings <- if(min(@gamma) < 0){
			NULL
		} else {
			_PSIML * @stuff$vectors[,@J] * sqrt(@gamma')
		}
		structure(loadings:delete(@loadings,return:T),\
				  eigenvals:delete(@stuff,return:T)$values)
	} else {
		NULL
	}
} else {
	@stuff <- releigen(@s, dmat(_PSIML))
	@gamma <- @stuff$values[-@J]
	if (min(@gamma) <= 0) {
		delete(@stuff)
		if (($2) == 0) {
			?
		} else {
			rep(?, @n)
		}
	} else {
		@crit <- sum(@gamma - 1 - log(@gamma))
		if (($2) == 0){
			@crit
		} else {
			@del <- if (ncomps(@param) > 2) {
				@param[3]
			} else {
				1e-5
			}
			if (@del > 0) {
				@gradient <- @x
				for (@i,1,@n){
					@tmp <- @x
					@tmp[@i] <- @tmp[@i] + @del
					@gamma <- releigenvals(@s, dmat(exp(@tmp)))[-@J]
					@gradient[@i] <- if (min(@gamma) <= 0) {
						?
					} else {
						(sum(@gamma - 1 - log(@gamma)) - @crit)/@del
					}
				}
				delete(@tmp,@i)
			} else {
				@omega <- sqrt(_PSIML)*@stuff$vectors[,-@J]
				@gradient <- vector(sum((1 - @gamma)*(@omega'^2)))
			}
			
			delete(@del,@stuff,@crit)
			delete(@gradient,return:T)
		}
	}
}

delete(@x,@s,@gamma,@crit,@J,@param,@n,silent:T)
delete(@result,return:T)
%mlcrit%

===> glscrit <===
glscrit       MACRO DOLLARS
) Macro to compute GLS factor analysis criterion function and its gradient
) The parameters x are log(psi_sub_i), i = 1,...,p
) This is to be used with optimization macros dfp(), bfs() and broyden()
) Usage
)  result <- glscrit(logpsi, k, structure(s, nfact [, del]))
)   logpsi       REAL vector of log uniquenesses
)   k            integer scalar to control what is calculated
)                k = 0, return criterion (scalar) (MISSING when s not P.D.)
)                k = 1, return gradient (vector) (rep(MISSING,nrows(s))
)                  when s not P.D.)
)                k = -2, return structure(loadings, eigenvals)
)   s            REAL symmetric covariance or correlation matrix
)   nfact > 0    integer number of factors
)   del >= 0     Small REAL scalar [1e-5].  Gradient is computed as
)     gradient[i] = (critfn(logpsi + del*(run(p)==i)) - critfn(logpsi))/del
)                del = 0 results in analytic gradient being used
)   When s not positive definite, MISSING (k = 0) or
)   rep(MISSING,nrows(s)) (k = 1) is returned
) No argument checking is done
)) Version 011120, C. Bingham, kb@stat.umn.edu
)) 011124 del = 0 signals use of analytic gradient
@x <- $1
@n <- length(@x)
_PSIGLS <- exp(@x)

@param <- $3
@s <- @param[1]
@J <- run(@param[2])
@result <- if (($2) < 0) {
	if (($2) == -2) {
		# return structure(loadings, eigenvalues)
		@stuff <- releigen(@s, dmat(_PSIGLS))
		@gamma <- @stuff$values[@J]-1
		@crit <- NULL
		@loadings <- if(min(@gamma) < 0){
			NULL
		} else {
			_PSIGLS * @stuff$vectors[,@J] * sqrt(@gamma')
		}
		structure(loadings:delete(@loadings,return:T),\
				  eigenvals:delete(@stuff,return:T)$values)
	} else {
		NULL
	}
} else {
	@stuff <- releigen(@s, dmat(_PSIGLS))
	@gamma <- @stuff$values[-@J]
	if (min(@gamma) <= 0) {
		delete(@stuff)
		if (($2) == 0) {
			?
		} else {
			rep(?, @n)
		}
	} else {
		@crit <- .5*sum((1 - 1/@gamma)^2)
		@result <- if (($2) == 0){
			@crit
		} else {
			@del <- if (ncomps(@param) > 2) {
				@param[3]
			} else {
				1e-5
			}
			if (@del > 0) {
				@gradient <- @x
				for (@i,1,@n){
					@tmp <- @x
					@tmp[@i] <- @tmp[@i] + @del
					@gamma <- releigenvals(@s, dmat(exp(@tmp)))[-@J]
					@gradient[@i] <- if (min(@gamma) <= 0) {
						?
					} else {
						@gradient[@i] <- (.5*sum((1 - 1/@gamma)^2) - @crit)/@del
					}
				}
				delete(@tmp,@i)
			} else {
				@omega <- sqrt(_PSIGLS)*@stuff$vectors[,-@J]
				@gradient <- vector(sum((1 - @gamma)*(@omega'^2)/@gamma^2))
			}
			delete(@del,@stuff)
			delete(@gradient,return:T)
		}
	}
}
delete(@x,@s,@gamma,@crit,@J,@n,silent:T)
delete(@result,return:T)
%glscrit%

===> ulscrit <===
ulscrit       MACRO DOLLARS
) Macro to compute ULS factor analysis criterion function and its gradient
) The parameters x are log(psi_sub_i), i = 1,...,p
) This is to be used with optimization macros dfp(), bfs() and broyden()
) Usage
)  result <- ulscrit(logpsi, k, structure(s, nfact [, del]))
)   logpsi       REAL vector of log uniquenesses
)   k            integer scalar to control what is calculated
)                k = 0, return criterion (scalar) (MISSING when s not P.D.)
)                k = 1, return gradient (vector) (rep(MISSING,nrows(s))
)                  when s not P.D.)
)                k = -2, return structure(loadings, eigenvals)
)   s            REAL symmetric covariance or correlation matrix
)   nfact > 0    integer number of factors
)   del >= 0     Small REAL scalar [1e-5].  Gradient is computed as
)     gradient[i] = (critfn(logpsi + del*(run(p)==i)) - critfn(logpsi))/del
)                Currently del = 0 gets the default.  At some point del = 0
)                will signal the use of analytic derivatives
)   When s not positive definite, MISSING (k = 0) or
)   rep(MISSING,nrows(s)) (k = 1) is returned
) No argument checking is done
)) Version 011120, C. Bingham, kb@stat.umn.edu
@x <- $1
@n <- length(@x)
_PSIULS <- exp(@x)

@param <- $3
@s <- @param[1]
@J <- run(@param[2])
@result <- if (($2) < 0) {
	if (($2) == -2) {
		# return structure(loadings, eigenvals)
		@crit <- NULL
		@stuff <- eigen(@s - dmat(_PSIULS))
		@delta <- @stuff$values[@J]
		@loadings <- if(min(@delta) < 0){
			NULL
		} else {
			@stuff$vectors[,@J] * sqrt(@delta')
		}
		structure(loadings:delete(@loadings,return:T),\
				  eigenvals:delete(@stuff,return:T)$values)
	} else {
		NULL
	}
} else {
	@stuff <- eigen(@s - dmat(_PSIULS))
	@delta <- @stuff$values[-@J]
	@crit <- .5*sum(@delta^2)
	if (($2) == 0){
		@crit
	} else {
		@del <- if (ncomps(@param) > 2) {
			@param[3]
		} else {
			1e-5
		}
		if (@del > 0) {
			delete(@stuff)
			@gradient <- @x
			for (@i,1,@n){
				@tmp <- @x
				@tmp[@i] <- @tmp[@i] + @del
				@delta <- eigenvals(@s - dmat(exp(@tmp)))[-@J]
				@gradient[@i] <- (.5*sum(@delta^2) - @crit)/@del
			}
			delete(@tmp,@i,@del)
		} else {
			@omega <- delete(@stuff,return:T)$vectors[,-@J]
			@gradient <-\
			   -_PSIULS*vector(sum(@delta*(delete(@omega,return:T)'^2)))
		}
		
		delete(@gradient,return:T)
	}
}

delete(@x,@s,@delta,@crit,@J,@n,silent:T)
delete(@result,return:T)
%ulscrit%

===> facanal <===
facanal   MACRO   DOLLARS
) Macro to use optimization macro dfp(), bfs() or broyden()
) to find ULS, GLS or ML estimates of uniquenesses and loadings.
) It uses macro dfp(), bfs() or broyden() to minimize the criterion as
) a function of log(psi), psi = uniquenesses.
)
) It requires macros ulscrit(), glscrit() or mlcrit() to compute
) the objective function, gradient vector and loading matrix.
)
) Usage:
)  facanal(s, m [, method:Meth] [,rotate:Rotmeth] \
)           [,start:psi0  or  presteps:nsteps]\
)           [,maxit:nmax] [,minit:nmin] [,crit:vector(nsig1,nsig2,dg)]\
)           [,minimizer:M] [,gold:ngold] [,quiet:T] [,silent:T] \
)           [,kaiser:F] [,recordwhen:d1] [,printwhen:d2])
)   s          REAL p x p positive definite covariance or correlation
)              matrix
)   m          0 < integer < (2*p + 1 - sqrt(8*@p+1))/2
)   Meth       one of "uls", "gls", "ml" or "mle" ["mle"]
)   Rotmeth    one of "varimax", "quartimax", "equimax" ["none"]
)   psi0       Positive length p vector of starting values for psi
)   nsteps     Integer > 0, number of stepuls(), stepgls() or stepml()
)              steps before using minimizer [1]; 'presteps' illegal with
)              'start:psi0'
)   nmax       integer > 0, maximum number of iterations permitted [30]
)   nmin       integer >= 0, minimum number of iterations [0]
)   nsig1      target number of significant digits in log uniquenesses [5]
)   nsig2      target number of significant digits in criterion [8]
)   dg         target upper limit for ||gradient|| [-1]
)   M          one of "dfp", "bfs", or "broyden" ["bfs"]
)   ngold      integer >= 0, # of steps in golden mean linear search [1]
)   quiet:T    suppresses printing of results
)   silent:T   suppresses printing of results and warning messages
)   kaiser:F   don't use Kaiser normalization
)   d1         integer >= 0. When d1 > 0, partial results are printed at
)              iterations d1, 2*d1, ...  [0]; not affected by silent:T
)   d2         integer >= 0. When d2 > 0, partial results are recorded at
)              iterations d2, 2*d2, ...  [0]
)
)   nsig1, nsig2 and dg specify convergence critera; any <= 0 are ignored
)   ngold is ignored if M is "broyden"
) The result is an "invisible" structure
) structure(psihat:psi,loadings:l,criterion:crit,eigenvals:vals,\
)           gradient:g, method:Meth, rotation:Rotmeth,iter:n, status:k)
)
)  psi         REAL vector of estimated uniquenesses
)  l           REAL p by m matrix of loadings
)  crit        Minimized criterion
)  vals        eigenvalues s relative to psi
)  g           REAL gradient vector at optimum
)  n > 0       integer = number of iterations
)  k >= 0      integer specifying which of the convergence criterion
)              signalled convergence; k == 0 means not converged
)
) When d2 > 0, side effect structure DFPRECORD, BFSRECORD or
) BROYDNRECORD with components funvals, xvals and gradients is created.
) On each call, ulscrit(), glscrit() or mlcrit() saves a copy of psi
) in invisible variable _PSIULS, _PSIGLS or _PSIML.
)
)) Version of 011117, written by C. Bingham, kb@stat.umn.edu
)) 011125 added keyword 'golden' with default 1
))        added keyword 'rotation' with default "none"
))        changed default for delta to 0 so as to use analytic gradient
))        modified printed output
))        result now includes method and rotation
)) 030115 minor bug fixed so silent:T is recognized
)) 031217 Fixed bug so it would recognize when macro minimizer was loaded
)) 031222 Fixed minor bug so it gets labels right when dim(s) is 1,p,p
)) 041111 Changed default starting value to result of one principal
))        factor step from 1/diag(solve(s)); added default labels to
))        PSI and LOADINGS; cosmetic changes to output
)) 041115 Changed default starting value to the result of one step away
)         from 1/diag(solve(s)) using stepuls(), stepgls() or stepml();
))        added keyword phrase 'presteps:nsteps' to allow more steps
# $S(s or r, m [, start:psi0  or presteps:n] [, maxit:m, silent:T, quiet:T])
@s <- matrix(argvalue($1,"variance or correlation","square real nonmissing"))
@m <- argvalue($2,"number of factors","positive count")
@keys <- if ($k > 0) {
	structure($K,NULL)
} else {
	structure(NoTaKey:NULL)
}

@method <- keyvalue(@keys,"method","string",default:"ml")
@maxit <- keyvalue(@keys, "maxit*","positive count",default:30)
@minit <- keyvalue(@keys, "minit*","count",default:0)
@delta <- keyvalue(@keys,"delta","nonnegative number",default:0)
@silent <- keyvalue(@keys, "silent", "TF", default:F)
@quiet <- keyvalue(@keys, "quiet", "TF", default:@silent)
@record <- keyvalue(@keys, "recordwhen", "count", default:0)
@print <- keyvalue(@keys, "printwhen", "count", default:0)
@minimizer <- keyvalue(@keys,"minimizer","string",default:"bfs")
@crit <- keyvalue(@keys,"crit","real nonmissing vector",default:vector(5,8,-1))
@gold <- keyvalue(@keys,"gold*","count",default:1)
@rotmethod <- keyvalue(@keys,"rot*","string",default:"none")
@kaiser <- keyvalue(@keys,"kais*","TF",default:T)
@presteps <- keyvalue(@keys,"prestep*","positive count")
@x0 <- keyvalue($K,"start","positive vector")
delete(@keys)

@imeth <- min(match(@method,vector("uls","gls","ml","mle"),0),3)
if (@imeth == 0) {
	error(paste("unrecognized method \"",@method,"\"",sep:""))
}
@kmethod <- match(vector("var*","quart*","equi*","none"),\
			  @rotmethod,0,exact:F)
if (sum(@kmethod) != 1) {
	error(paste(@rotmethod,\
				"is not one of \"varimax\",\"quartimax\" or \"equimax\""))
}
@kmethod <- run(4)[@kmethod != 0]
@rotmethod <- \
	vector("Varimax","Quartimax","Equimax","none")[delete(@kmethod,return:T)]
if (@minimizer == "broyden") {
	@gold <- 0
}

@p <- nrows(@s)
if (isnull(@x0)) {
	if (isnull(@presteps)){
		@presteps <- 1
	}
	@stepper <- vector("stepuls","stepgls","stepml")[@imeth]
	if (!ismacro(<<@stepper>>)) {
		if (@imeth == 1){
			getmacros(stepuls,silent:T)
		} elseif (@imeth == 2) {
			getmacros(stepgls,silent:T)
		} else {
			getmacros(stepml,silent:T)
		}
	}
} else {
	if (!isnull(@presteps)) {
		error("You can't use both 'presteps:nsteps' and 'start:psi0'")
	}
	if (nrows(@x0) != @p) {
		error("Wrong number of starting uniquenesses")
	}
}

@kind <- if(sum(diag(@s) == 1) == @p){
	@kind1 <- "R"
	"Correlation"
}else{
	@kind1 <- "S"
	"Covariance"
}
if (min(eigenvals(@s)) <= 0){
	error(paste(@kind,"matrix is not positive definite"))
}
if (F && !ismacro(<<@minimizer>>)){
	if (@minimizer == "dfp") {
		getmacros(dfp,quiet:T,printname:F)
	} elseif (@minimizer == "bfs") {
		getmacros(bfs,quiet:T,printname:F)
	} elseif (@minimizer == "broyden") {
		getmacros(broyden,quiet:T,printname:F)
	} else {
		error(paste("unrecognized 'minimizer' \"",@minimizer,"\"",sep:""))
	}
} elseif (!ismacro(minimizer)){
		getmacros(minimizer,quiet:T,printname:F)
}

@critfn <- vector("ulscrit","glscrit","mlcrit")[@imeth]
if (!ismacro(<<@critfn>>)) {
	if (@imeth == 1){
		getmacros(ulscrit,quiet:T,printname:F)
	} elseif (@imeth == 2){
		getmacros(glscrit,quiet:T,printname:F)
	} else {
		getmacros(mlcrit,quiet:T,printname:F)
	}
}

if (isnull(@x0)){ # no value; compute starting value
	# start with 1 or more steps of stepuls(), stepgls() or stepml()
	@x0 <- @psi0 <- 1/diag(solve(@s))
	for (@i,1,@presteps) {
		@x0 <- <<@stepper>>(@s, @x0, @m)
	}
	@x0 <- @x0$psi
	if (min(@x0) <= 0) {
		@x0 <- @psi0
	}
	delete(@psi0)
}

@x0 <- log(@x0)
 #@x0 <- log(keyvalue($K,"start","positive vector",\
 #		   default:1/diag(solve(@s)))) # starting value


@param <- structure(@s, @m, @delta)
if (F) {
	@result <- <<@minimizer>>(@x0,<<@critfn>>,@param,\
		maxit:@maxit, minit:@minit, crit:@crit,gold:@gold,\
		recordwhen:@record, printwhen:@print)
} else {
	@result <- minimizer(@x0,<<@critfn>>,@param,method:@minimizer,\
		maxit:@maxit, minit:@minit, crit:@crit,gold:@gold,\
		recordwhen:@record, printwhen:@print,silent:@silent)
}

 # The result is structure(x:xmin, f:minVal, gradient:gradient,\
 # h:invhessian, iterations:niter)

@iter <- @result$iterations
@status <- @result$status
if (@status < 0) {
	print(paste("WARNING: criterion or gradient not defined on iteration", \
			@iter))
}

@x <- @result$x
PSI <- vector(exp(@x))
CRITERION <- @result$f
GRADIENT <- vector(@result$gradient)

@result <- if (@imeth == 1) {
	ulscrit(@x, -2, @param)
} elseif (@imeth == 2){
	glscrit(@x, -2, @param)
} else {
	mlcrit(@x, -2, @param)
}
LOADINGS <- if (@rotmethod == "none") {
	@result$loadings
} else {
	@H <- if (!@kaiser) {
		1
	} else {
		vector(sqrt(sum(@result$loadings'^2)))
	}
	@H*rotation(@result$loadings/@H, method:@rotmethod, reorder:T)
}

EIGENVALUES <- delete(@result,return:T)$eigenvals
delete(@x,@param,@kind)

@labels <- if (haslabels(@s)) {
	getlabels(@s,1)
} else {
	"V"
}
@labels <- structure(@labels,vector("F1","F")[1 + (@m > 1)])
setlabels(PSI, @labels[1])
setlabels(LOADINGS, delete(@labels,return:T))

if (!@silent) {
	if (2*@m >= 2*@p + 1 - sqrt(8*@p+1)){
		print(paste("WARNING: With m = ",@m," and p = ",@p,\
					", m > (2*p + 1 - sqrt(8*p+1))/2",sep:""))
	}
	if (isnull(LOADINGS)){
		print(paste("WARNING: eigenvalue",@m,"of Psi^(-1) %*%",@kind1,"< 1"))
	}
	if (@status <= 0){
		print(paste("WARNING: no convergence in",@iter,"iterations"))
	}
	if (!@quiet) {
		if(@status > 0){
			print(paste("Convergence in",@iter,"iterations by criterion",@status))
		}
		@loadname <- if (@rotmethod != "none") {
			paste(@rotmethod,"rotated estimated loadings")
		} else {
			"Unrotated estimated loadings"
		}

		print(PSI,name:"Estimated uniquenesses")
		print(LOADINGS, name:delete(@loadname,return:T))
		print(CRITERION, name:paste("Minimized",@method,"criterion"))
	}
}
@silent <- @silent || @quiet
delete(@quiet,@p,@m,@delta,@kind1)
@result <- structure(psihat:PSI,loadings:LOADINGS,criterion:CRITERION,\
	eigenvals:EIGENVALUES, gradient:GRADIENT,\
	method:delete(@method,return:T),rotation:delete(@rotmethod,return:T),\
	iter:delete(@iter,return:T),status:delete(@status,return:T))
delete(@result,return:T,invisible:!delete(@silent,return:T))
%facanal%

===> ulsfactor <===
ulsfactor   MACRO   DOLLARS
) Macro to use nonlinear least squares macro levmar() to find
) ULS estimates of uniquenesses and loadings.  It uses macro ulsresids()
) to compute residuals used by levmar().
)
) Usage:
)	result <- ulsfactor(s, m [, quiet:T,silent:T,start:psi0,uselogs:T],\
)       maxit:nmax, minit:nmin, print:T,crit:crvec])
)   s          positive definite covariance or correlation matrix
)   m          positive integer < (2*p + 1 - sqrt(8*@p+1))/2
)   nmax       nonnegative integer, maximum number of iterations permitted
)              default = 30
)   nmin       nonnegative integer, minimum of iterations required,
)              default=0 [0]
)   psi0       positive REAL vector of length p, starting values
)   crvec      vector(numsig, nsigsq, delta), 3 criteria for convergence;
)              criterion < means it is not used
)   print:T    something is printed on each interation
)   uselogs:T  parametrize in terms of log(psi)
)   result     "invisible" structure(psihat:psi,loadings:l,criterion:crit,\
)              eigenvals:vals,status:iconv,iter:n)
)    psi       vector of estimated uniquenesses
)    l         p by m matrix of loadings with orthogonal columns
)    crit      sum(vector((s-sigmahat)^2)), sigmahat = l%*%l' + dmat(psi)
)    vals      vector of eigenvalues of s - psi
)    iconv     convergence information
)        = 0 or 4 not converged
)        = 1, converged with relative change in coeffs < relcon = 10^-numsig
)        = 2, converged with relative change in rss <= relssq = 10^-nsigsq
)        = 3, converged with ||gradient|| < delta
)        = 4, interation stopped; could not reduce rss.
)    n         number of iterations required
)
) There is no guarantee that min(psihat) > = 0.  In that case you can try
)   using uselogs:T to parametrize in terms of log(psi).  This ensures
)    min(psi) > 0 if it converges.
)   However, if the minimum is actually outside the permissible
)    range, you will probably get a 'argument to solve() singular'
)    error message.
)
) Macro ulsresids saves a copy of psi on every call in invisible variable
) _PSIULS.  Type 'print(_PSIULS) to see this value.
)) Version of 991114, written by C. Bingham, kb@stat.umn.edu
)) It uses macro levmar to minimize the sum of squared residuals.
# $S(s or r, m [,start:psi0,maxit:m,print:T,quiet:T])
@s <- argvalue($1,"variance or correlation","square real nonmissing")

@quiet <- keyvalue($K, "quiet","TF",default:F)
@silent <- keyvalue($K, "silent","TF",default:F)
@maxit <- keyvalue($K, "maxit*","positive integer",default:30)

__USELOGSULS <- keyvalue($K, "uselogs","TF",default:F)

@m <- argvalue($2,"number of factors","positive integer")
@p <- nrows(@s)
@kind <- if(sum(diag(@s) == 1) == @p){
	"correlation"
}else{
	"covariance"
}
if (min(eigenvals(@s)) <= 0){
	error(paste(@kind,"matrix is not positive definite"))
}
if (2*@m >= 2*@p + 1 - sqrt(8*@p+1)){
	error(paste("With m = ",@m," and p = ",@p,\
		", m > (2*p + 1 - sqrt(8*p+1))/2",sep:""))
}

if (!ismacro(levmar)){
	getmacros(levmar,quiet:T,printname:F)
}
if (!ismacro(ulsresids)){
	getmacros(ulsresids,quiet:T,printname:F)
}

@psi0 <- keyvalue($K,"start","positive vector",\
	default:1/diag(solve(@s))) # starting value
if (length(@psi0) != @p){
	error("Wrong number of starting uniqueness")
}
if (__USELOGSULS){
	@psi0 <- log(@psi0)
}
@result <- levmar(@psi0,0,vector(@s),@m,resid:ulsresids,$K)

@psihat <- if (__USELOGSULS){
	exp(@result$coefs)
}else{
	@result$coefs
}
@iconv <- @result$iconv
@iter <- @result$iter

if (@iconv == 0){
	print(paste("WARNING: failed to converge in",@maxit,"iterations"))
}
@eigs <- eigen(@s - dmat(@psihat))
@vals <- @eigs$values
@J <- run(@m)
@lhat <- if (@vals[@m] < 0){
	NULL
}else{
	@eigs$vectors[,@J] * sqrt(@vals[@J]')
}
delete(@eigs,@J)

if (!@silent){
	if (isnull(@lhat)){
		print(paste("WARNING: eigenvalue",@m,"of (R or S) - Psi is negative"))

	}
	if (min(@psihat) < 0){
		print("WARNING: psi not acceptable since min(psihat) < 0")
	}
	if (@iconv == 4){
		print("WARNING: could not reduce RSS on iteration",@iter)
	}elseif (@iconv <= 0){
		print(paste("WARNING: no convergence in",@iter,"iterations"))
	}
}
if (!@quiet && !@silent){
	if(@iconv > 0 && @iconv < 4){
		print(paste("Convergence in",@iter,"iterations by criterion",@iconv))
	}	
	print(psihat:@psihat,loadings:@lhat,criterion:@result$rss,eigenvals:@vals)
}

delete(@quiet,@silent,@p,@m,@psi0,__USELOGSULS)
@result <- structure(psihat:delete(@psihat,return:T),\
	loadings:delete(@lhat,return:T),criterion:delete(@result$rss,return:T),\
	eigenvals:delete(@vals,return:T),status:delete(@iconv,return:T),\
	iter:delete(@iter,return:T))
delete(@result,return:T,invisible:T)
%ulsfactor%

===> ulsresids <===
ulsresids    MACRO    DOLLARS OUTLINE
) Macro to be used with levmar to do ULS factor extraction
)
) ulsresids() usage:
)   ulsresids(psi,0,vector(r),m)
)     psi        vector of uniquenesses
)     r          correlation or covariance matrix
)     m          positive integer number of factors
)
) levmar usage:
)   levmar(if(__USELOGSULS){log(psi)}else{psi},0,\
)     vector(r),m,resid:ulsresids)
)) 0000828
# $S(psi,0,vector(r),m)
_PSIULS <- if(__USELOGSULS){
	exp($1)
}else{
	$1
}
@psi <- dmat(_PSIULS)
@r <- $3
@p <- sqrt(length(@r))
@J <- run($4) #run(m)
@eigs <- eigen(matrix(@r,@p)  - @psi)
@l <- @eigs$vectors[,@J]
@res <- delete(@r,return:T) - vector((@l * @eigs$values[@J]') %C% @l + @psi)
delete(@psi,@p,@J,@eigs,@l)
delete(@res, return:T)
%ulsresids%

===> glsfactor <===
glsfactor   MACRO   DOLLARS
) Macro to use nonlinear least squares macro levmar to find
) GLS estimates of uniquenesses and loadings.  It uses macro glsresids
) to compute residuals used by levmar.
)
) Usage:
)	result <- glsfactor(s, m [, quiet:T,silent:T,start:psi0,uselogs:T],\
)       maxit:nmax, minit:nmin, print:T,crit:crvec])
)   s          positive definite covariance or correlation matrix
)   m          positive integer < (2*p + 1 - sqrt(8*@p+1))/2
)   nmax       nonnegative integer, maximum number of iterations permitted
)              default = 30
)   nmin       nonnegative integer, minimum of iterations required,
)              default=0 [0]
)   psi0       positive REAL vector of length p, starting values
)   crvec      vector(numsig, nsigsq, delta), 3 criteria for convergence;
)              criterion < means it is not used
)   print:T    something is printed on each interation
)   uselogs:T  parametrize in terms of log(psi)
)   result     "invisible" structure(psihat:psi,loadings:l,criterion:crit,\
)              eigenvals:vals,status:iconv,iter:n)
)    psi       vector of estimated uniquenesses
)    l         p by m matrix of loadings with orthogonal columns satisfying
)              l' %*% dmat(1/psi) %*% l is diagonal
)    crit      sum(vector((s-sigmahat)^2)), sigmahat = l%*%l' + dmat(psi)
)    vals      vector of eigenvalues of s - psi
)    iconv     convergence information
)        = 0 or 4 not converged
)        = 1, converged with relative change in coeffs < relcon = 10^-numsig
)        = 2, converged with relative change in rss <= relssq = 10^-nsigsq
)        = 3, converged with ||gradient|| < delta
)        = 4, interation stopped; could not reduce rss.
)
) There is no guarantee that min(psihat) > = 0.  In that case you can try
)   using uselogs:T to parametrize in terms of log(psi).  This ensures
)    min(psi) > 0 if it converges.
)   However, if the minimum is actually outside the permissible
)    range, you will probably get a 'argument to solve() singular'
)    error message.
)
)  l is the p by m matrix of loadings satisfying
)  crit is sum of squared differences of s from sigmathat
)  vals are eigenvalues of s - psi
)
) Macro glsresids saves a copy of psi on every call in invisible variable
) _PSIGLS.  Type 'print(_PSIGLS) to see this value.
)) Version of 991115, written by C. Bingham, kb@stat.umn.edu
)) It uses macro levmar to minimize the sum of squared residuals.
# $S(s or r, m [,start:psi0,maxit:m,print:T,quiet:T])
@s <- argvalue($1,"variance or correlation","square real nonmissing")

@quiet <- keyvalue($K, "quiet","TF",default:F)
@silent <- keyvalue($K, "silent","TF",default:F)
@maxit <- keyvalue($K, "maxit*","positive integer",default:30)

__USELOGSGLS <- keyvalue($K, "uselog*","TF",default:F)

@m <- argvalue($2,"number of factors","positive integer")
@p <- nrows(@s)
@kind <- if(sum(diag(@s) == 1) == @p){
	@what <- "r"
	"correlation"
}else{
	@what <- "s"
	"covariance"
}
@kind <- if(sum(diag(@s) == 1) == @p){
	"correlation"
}else{
	"covariance"
}
if (min(eigenvals(@s)) <= 0){
	error(paste(@kind,"matrix is not positive definite"))
}
if (2*@m >= 2*@p + 1 - sqrt(8*@p+1)){
	error(paste("With m = ",@m," and p = ",@p,\
		", m > (2*p + 1 - sqrt(8*p+1))/2",sep:""))
}

if (!ismacro(levmar)){
	getmacros(levmar,quiet:T,printname:F)
}
if (!ismacro(glsresids)){
	getmacros(glsresids,quiet:T,printname:F)
}

@psi0 <- keyvalue($K,"start","positive vector",\
	default:1/diag(solve(@s))) # starting value
if (length(@psi0) != @p){
	error("Wrong number of starting uniqueness")
}
if (__USELOGSGLS){
	@psi0 <- log(@psi0)
}
@result <- levmar(@psi0,vector(@s),vector(@s),@m,resid:glsresids,$K)

@psihat <- if (__USELOGSGLS){
	exp(@result$coefs)
}else{
	@result$coefs
}
@iconv <- @result$iconv
@iter <- @result$iter

if (@iconv == 0){
	print(paste("WARNING: failed to converge in",@maxit,"iterations"))
}
@eigs <- releigen(solve(@s), dmat(1/@psihat))
@J <- run(@m)
@vals <- 1/reverse(@eigs$values) - 1
@lhat <- if (@vals[@m] < 0){
	NULL
}else{
	@eigs$vectors[,@p + 1 - @J] * sqrt(@vals[@J]')
}
delete(@eigs,@J)

if (!@silent){
	if (isnull(@lhat)){
		print(paste("WARNING: eigenvalue",@m,"is negative"))
	}
	if (min(@psihat) < 0){
		print("WARNING: psi not acceptable since min(psihat) < 0")
	} elseif (min(eigenvals(@s - dmat(@psihat))) < 0) {
		print(paste("WARNING:",@what,\
			"- dmat(psihat) not non-negative definite"),macroname:T)
	}
	if (@iconv == 4){
		print("WARNING: could not reduce RSS on iteration",@iter)
	}elseif (@iconv <= 0){
		print(paste("WARNING: no convergence in",@iter,"iterations"))
	}
}
if (!@quiet && !@silent){
	if(@iconv > 0 && @iconv < 4){
		print(paste("Convergence in",@iter,"iterations by criterion",@iconv))
	}	
	print(psihat:@psihat,loadings:@lhat,criterion:@result$rss,vals:@vals)
}

delete(@quiet,@silent,@p,@m,@psi0,__USELOGSGLS)
@result <- structure(psihat:delete(@psihat,return:T),\
	loadings:delete(@lhat,return:T),criterion:delete(@result$rss,return:T),\
	vals:delete(@vals,return:T),status:delete(@iconv,return:T),\
	iter:delete(@iter,return:T))
delete(@result,return:T,invisible:T)
%glsfactor%

===> glsresids <===
glsresids   MACRO    DOLLARS
) Macro to be used with levmar to do GLS factor extraction
) This is very much a work in progress
) Usage:
)    stuff <- levmar(if(__USELOGSGLS){log(psi)}else{psi},\
)    vector(r),vector(r),m,resid:glsresids), m = number of factors
) glsresids itself returns vector(Ip - solve(r,sigmahat))
) It saves its most recent argument as _PSIGLS
)
) You need to set __USELOGSGLS to T or F before calling levmar
) and use log(psi) as an argument to levmar when __USELOGSGLS is True
) This is done automatically by macro glsfactor
)
) The result from levmar is the structure returned by levmar.
)   stuff$coefs is psihat (or log(psihat))
)   stuff$rss is the sum(vector(Ip - solve(r,sigmahat))^2)
)   stuff$iconv is convergence status -- 0 means not converged
)   stuff$iter is the number of iterations
) Version 991115
# $S(psi,vector(r),vector(r),m)
 # residualstr((Ip - solve(s,sigmahat))^2)
_PSIGLS <- if(__USELOGSGLS){
	exp($1)
}else{
	$1
}
@r <- $2 #x
@p <- sqrt(length(@r))
@r <- matrix(@r,@p)
@m <- $4
@eigs <- releigen(solve(@r), dmat(1/_PSIGLS))
@J <- run(@p, @p - @m + 1)
@tmpl <- @eigs$vectors[,@J]
@vals <- 1/@eigs$values[@J] - 1
@sighat <- (@tmpl * @vals') %C% @tmpl + dmat(_PSIGLS)
@resids <- sqrt(abs(dmat(@p,1) - solve(@r, @sighat)))
@resids <- @resids * @resids'
delete(@m,@J,@eigs,@tmpl,@vals,@p,@r,@sighat)
vector(delete(@resids, return:T))
%glsresids%

===> goodfit <===
goodfit             MACRO DOLLARS
) Usage:
)  goodfit(S, lhat, psihat)
)  goodfit(R, lhat, psihat)
)   S             p by p sample variance matrix
)   R             p by p sample correlation matrix
)   lhat          p by m matrix of factor loadings
)   psihat        length p vector of factor uniquenesses
)
) Prints the maximum absolute deviation, the sum of squared
) deviations and the sum of absolution deviations from S or R of the
) fitted covariance or correlation matrix S_hat or R_hat computed
) as lhat %*% lhat' + dmat(psihat)
)
))Version 991114
# usage: $S(s,lhat,psihat), s a covar or correl matrix
@lhat <- argvalue($2,"lhat","real matrix nonmissing")
@dev <- vector(argvalue($1,"s","real square") -\
	(@lhat %C% @lhat)  -\
	dmat(argvalue($3[[1]],"psi","nonmissing real vector")))

@tmp <- vector(max(abs(@dev)),sum(@dev^2), sum(abs(@dev)),\
	labels:vector("max(|dev|)","sum(dev^2)","sum(|dev|)"))
delete(@dev)
delete(@tmp,return:T)
%goodfit%

===> stepuls <===
stepuls            MACRO DOLLARS
) Macro to do 1 step of crude iteration for ULS factor analysis
) Usage:
)    result <- stepuls(s, psi, m [,nsteps:n] [,print:T])
)    result <- stepuls(s, structure(psi,...), m [,nsteps:n][,print:T])
)  s          p by p variance or correlation matrix
)  psi        length p vector of current values for uniquenesses
)  m          positive integer = number of factors
)  n          positive integer = number of steps to take, default 1
)  print:T    new values of uniqueness and loadings printed on
)             each step
)  result     structure(psi:newPsi, loadings:oldL, crit:criterion)
)   newPsi    length p vector of updated uniquenesses
)   oldL      p by m matrices of loadings computed from psi (not
)             from newPsi)
)   criterion ULS criterion computed from newPsi
)
) stepuls() also creates side effect variabls PSI, LOADINGS, and
) CRITERION matching the components of result or values at last
) successful step when it can't complete a step
))041111 added keyword 'nsteps'
))Version 041111
# usage: $S(r,psi,m [,print:T]), s=covar or corr, psi=current value
# m = order = number of factors
if ($v != 3 || $k > 2){
	error("usage: $S(s, psi, m [,nsteps:n] [,print:T])",macroname:F)
}
@s <- matrix(argvalue($1,"s","real square"))
@psi <- argvalue($2[[1]],"psi","positive vector")
@m <- argvalue($3,"number of factors","positive integer scalar")
@Jm <- run(@m)
@print <- keyvalue($K,"print","TF",default:F)
@nsteps <- keyvalue($K,"nstep*","positive count",default:1)
@sii <- diag(@s)
@labels <- if (haslabels(@s)) {
	getlabels(@s,1)
} else {
	"V"
}
@labels <- structure(@labels, vector("F1","F")[1 + (@m > 1)])
for (@istep, 1, @nsteps) {
	@eigs <- eigen(@s - dmat(@psi))
	@vals <- @eigs$values[@Jm]
	if(min(@vals)<= 0){
		error(paste("s - psi not positive definite on step",@istep))
	}
	@l <- sqrt(@vals)' * @eigs$vectors[,@Jm]
	@psi <- @sii - vector(sum(@l'^2))
	if(min(@psi) <= 0){
		error(paste("min(psi) <= 0 on step", @istep))
	}
	PSI <- @psi
	setlabels(PSI,@labels[1])
	LOADINGS <- @l
	setlabels(LOADINGS,@labels)
	CRITERION <- sum(@eigs$values[-@Jm]^2)
	if(@print){ # print goodness of fit & psi
		print(psi:PSI ,criterion:CRITERION)
	}
}
delete(@sii,@s,@nsteps,@eigs,@vals,@psi,@l,@Jm,@m,@istep,@print, @labels)
structure(psi:PSI, loadings:LOADINGS, crit:CRITERION)
%stepuls%

===> stepml <===
stepml             MACRO DOLLARS
) Macro to do 1 or more steps of crude iteration for MLE factor analysis
) Usage:
)    result <- stepml(s, psi, m [,nsteps:n] [,print:T])
)    result <- stepml(s, structure(psi,...), m [,nsteps:n][,print:T])
)  s          p by p variance or correlation matrix
)  psi        length p vector of current values for uniquenesses
)  m          positive integer = number of factors
)  n          positive integer = number of steps (default 1)
)  print:T    new values of uniqueness and loadings printed
)  result     structure(psi:newPsi, loadings:oldL, crit:criterion)
)   newPsi    length p vector of updated uniquenesses
)   oldL      p by m matrices of loadings computed from psi (not
)             from newPsi)
)   criterion ML criterion computed from newPsi
)
) stepml() also creates side effect variabls PSI, LOADINGS, and
) CRITERION matching the components of result.
))041111 added keyword nsteps; use eigen() instead of releigen()
))Version 041111
# usage: $S(s,psi,m [,print:T]), s=covar or corr, psi=current value
# m=order, print=T or F
if ($v != 3 || $k > 2){
	error("usage: $S(s, psi, m [,nsteps:n] [,print:T])",macroname:F)
}
@s <- matrix(argvalue($1,"argument 1","real square"))
@m <- argvalue($3,"number of factors","positive integer scalar")
@Jm <- run(@m)
@psi <- argvalue($2[[1]],"psi","positive vector")
@print <- keyvalue($K,"print","TF",default:F)
@nsteps <- keyvalue($K,"nstep*","positive count",default:1)
@sii <- diag(@s)
@labels <- if (haslabels(@s)) {
	getlabels(@s,1)
} else {
	"V"
}
@labels <- structure(@labels, vector("F1","F")[1 + (@m > 1)])

for(@istep, 1, @nsteps) {
	@eigs <- eigen(@s/sqrt(@psi*@psi'))
	@vals <- @eigs$values[@Jm] - 1

	if(min(@vals) <= 0){
		error(paste("s - psi not positive definite on step",@istep))
	}
	@l <- sqrt(@psi) * @eigs$vectors[,@Jm] * sqrt(@vals)'

	@psi <- @sii - vector(sum(@l'^2))
	if(min(@psi) <= 0){
		error(paste("non-positive psi on step", @istep))
	}
	CRITERION <- sum((@eigs$values - log(@eigs$values)-1)[-@Jm])
	PSI <- @psi
	setlabels(PSI,@labels[1])
	LOADINGS <- @l
	setlabels(LOADINGS,@labels)
	if(@print){ # print goodness of fit & psi
		print(psi:PSI ,criterion:CRITERION)
	}
}

delete(@s,@m,@Jm,@psi,@print,@nsteps,@sii,@istep,@eigs,@vals,@l, @labels)
structure(psi:PSI, loadings:LOADINGS, crit:CRITERION)
%stepml%

===> stepgls <===
stepgls            MACRO DOLLARS
) Macro to do 1 step of crude iteration for GLS factor analysis
) Usage:
)    result <- stepgls(s, psi, m [,nsteps:n] [,print:T])
)    result <- stepgls(s, structure(psi,...), m [,nsteps:n][,print:T])
)  s          p by p variance or correlation matrix
)  psi        length p vector of current values for uniquenesses
)  m          positive integer = number of factors
)  n          positive integer = number of steps to take, default 1
)  print:T    new values of uniqueness and loadings printed
)  result     structure(psi:newPsi, loadings:oldL, crit:criterion)
)   newPsi    length p vector of updated uniquenesses
)   oldL      p by m matrices of loadings computed from psi (not
)             from newPsi)
)   criterion GLS criterion computed from newPsi
)
) stepgls() also creates side effect variabls PSI, LOADINGS, and
) CRITERION matching the components of result.
))041111 added keyword 'nsteps' and use eigen() instead of releigen()
))Version 041111
# usage: $S(r, psi, m [,print:T]), s=covar or corr, psi=current value
# m = order = number of factors
if ($v != 3 || $k > 2){
	error("usage: $S(s, psi, m [,nsteps:n] [,print:T])",macroname:F)
}
@s <- matrix(argvalue($1,"argument 1","real square"))
@psi <- argvalue($2[[1]],"psi","positive vector")
@m <- argvalue($3,"number of factors","positive integer scalar")
@Jm <- run(@m)
@print <- keyvalue($K,"print","TF",default:F)
@nsteps <- keyvalue($K,"nstep*","positive count",default:1)

@labels <- if (haslabels(@s)) {
	getlabels(@s,1)
} else {
	"V"
}
@labels <- structure(@labels, vector("F1","F")[1 + (@m > 1)])

for (@istep, 1, @nsteps) {
	@eigs <- eigen(@s/sqrt(@psi*@psi'))
	@vals <- @eigs$values[@Jm] - 1
	if(min(@vals) <= 0){
		error(paste("s - psi not positive definite on step", @istep))
	}
	@l <- sqrt(@psi) * @eigs$vectors[,@Jm] * sqrt(@vals)'

	@sinv <- solve(@s)
	@psi <- vector(solve(@sinv^2,\
		diag(@sinv) - vector(sum((@l %c% @sinv)^2))))
	if(min(@psi) <= 0){
		error(paste("non-positive psi; cannot continue on step",@istep))
	}
	CRITERION <- sum((1-1/@eigs$values[-@Jm])^2)
	PSI <- @psi
	setlabels(PSI, @labels[1])
	LOADINGS <- @l
	setlabels(LOADINGS, @labels)
	if(@print){ # print goodness of fit & psi
		print(psi:PSI ,criterion:CRITERION)
	}
}
delete(@eigs, @sinv, @s, @psi, @l, @nsteps, @print, @istep, @Jm, @m, @labels)
structure(psi:PSI, loadings:LOADINGS, crit:CRITERION)
%stepgls%


===> compf <===
compf              MACRO DOLLARS
) Usage (only when integer vector INS is defined):
)  f <- compf(h,e,fh,fe)
)   h        p by p MANOVA hypothesis matrix
)   e        p by p MANOVA error matrix
)   fh       integer scalar hypothesis degrees of freedom
)   fe       integer scalar error degrees of freedom
)   f        REAL vector of 'F-to-enter' statistics for each
)            "out" variable in a forward stepwise discriminant
)            analysis, when the variables listed in integer vector
)            INS are "in".
) When INS[1] = 0, no variables are "in"
)
) compf() is obsolete.  Use dastepsetup(), dastepstatus(), daentervar()
) daremovevar() and dasteplook(), or dastepselect() to do stepwise
) discriminant analysis
)) Version 000529, use argvalue()
# usage: f <- $S(h,e,fh,fe) ; INS must be defined
@h <- matrix(argvalue($1,"hypothesis matrix", "real square nonmissing"))
@e <- matrix(argvalue($2,"error matrix", "real square nonmissing"))
@fh <- argvalue($3, "hypothesis DF", "positive count")
@fe <- argvalue($4, "error DF", "positive count")
@p <- nrows(@h)
if (nrows(@e) != @p){
	error("hypothesis and error matrices are different sizes")
}
if (alltrue(!isvector(INS,real:T), !isnull(INS))){
	print("WARNING: INS not defined; initialized to 0")
	INS <- 0
} elseif (!isnull(INS)) {
	if(sum(INS != ceiling(INS)) != 0 || min(INS) < 0 || max(INS) > @p ||\
			!isscalar(INS) && min(INS) == 0){
		error("INS not 0 or vector of positive integers <= nrows(h)")
	}
	if(length(unique(INS)) != length(INS)){
		error("duplicate elements in INS")
	}
	if(@p < length(INS)) {
		error("length(INS) > p = nrows(h); not possible")
	}
}
if (alltrue(!isnull(INS), @p == length(INS))){
	print("WARNING: all variables IN", macroname:T)
	@out <- structure(f:NULL,df:vector(@fh,@fe - @p),ins:INS,OUTS:NULL)
} else {
	if(anytrue(isnull(INS), INS[1]==0)){ # initialize
		@h <- diag(@h)
		@e <- diag(@e)
		OUTS <- run(@p)
	} else {
		OUTS <- run(@p)[-INS]
		@h <- diag(swp(@h + @e,INS))[OUTS]
		@e <- diag(swp(@e,INS))[OUTS]
		@h <-- @e
		@fe <-- length(INS)
	}
	@out <- structure(f:(@h/@fh)/(@e/@fe),df:vector(@fh,@fe),ins:INS,outs:OUTS)
}
delete(@p,@fh,@fe,@h,@e)
delete(@out,return:T)
%compf%

===> forstep <===
forstep            MACRO DOLLARS
) forstep(i,h,e,fh,fe)
) Macro to add i-th variable to list of "in" variables in multigroup
) stepwise discriminant analysis.
) Integer vector INS must be defined and is a list of numbers of
) variables already "in".  INS = 0 means no variables are "in"
)
) Returns structure(f, df, ins, outs) where f is the vector of
) F-to-enter values, df = vector(fh,fe) contains their degrees
) of freedom, ins and outs are vectors of variable numbers that
) are "in" and "out" respectively.  INS = 0 means no variables are "in"
)
) forstep() is obsolete.  Use dastepsetup(), dastepstatus(),
) daentervar(), daremovevar() and dasteplook(), or dastepselect()
) to do stepwise discriminant analysis
)) Version 000826 use argvalue()
# usage: $S(i,h,e,fh,fe) ; INS must be defined
if ($v != 5){
	error("usage is $S(ivar,h,e,fh,fe)")
}
@h <- matrix(argvalue($2,"hypothesis matrix","nonmissingreal square"))
@e <- matrix(argvalue($3,"error matrix","nonmissingreal square"))
if(nrows(@h) != nrows(@e)){
	error("hypothesis and error matrices are different sizes")
}
@fh <- argvalue($4,"hypothesis DF","positive count")
@fe <- argvalue($5,"error DF","positive count")
@i <- argvalue($1,"entering variable number","positive count")

if (!isvector(INS, real:T) && !isnull(INS)){
	print("WARNING: initializing INS to 0")
	INS <- 0
}

if (alltrue(!isnull(INS), match(@i, INS, 0) != 0)){
	error("variable already IN")
}
INS <- if(alltrue(!isnull(INS), INS[1] != 0)){
	vector(INS,@i)
}else{
	@i
}
OUTS <- if(isscalar(OUTS) == 1){
	0
}else{
	OUTS[OUTS != @i]
}
if (!ismacro(compf)){
	getmacros(compf,quiet:T,printname:F)
}
@out <- compf(@h, @e, @fh, @fe)
delete(@h, @e, @fh, @fe, @i)
delete(@out, return:T)
%forstep%

===> backstep <===
backstep           MACRO DOLLARS
) Macro for doing a back step in stepwise discriminant analysis.
) Usage:
)   backstep(h,e,fh,fe)
)   INS must be defined and contain indices of in variables
)   h and e are hypothesis and error matrices with fh and fe degrees of
)   freedom
) backstep deletes the variable with smallest F-to-delete, and returns
) structure(f:F_to_delete,df:vector(fh,fe),ins:INS,outs:OUTS)
) INS is updated and OUTS is created
)) Version 000826 Use argvalue()
# usage: $S(h,e,fh,fe)
if ($v != 4){
	error("usage is $S(h,e,fe,fh)")
}
@h <- matrix(argvalue($1,"hypothesis matrix","nonmissing square real"))
@e <- matrix(argvalue($2,"error matrix","nonmissing square real"))
@fh <- argvalue($3,"hypothesis DF","positive count")
@fe <- argvalue($4,"error DF","positive count")
@p <- nrows(@h)
if (!isvector(INS,real:T)){
	print(paste("WARNING: INS not defined; set to run(",@p,")",sep:""))
	INS <- run(@p)
}
if(anytrue(isnull(INS), INS[1] == 0)){
	error("no variables left IN")
}
   # get submatrices of h & e for variables in INS
@h <- @h[INS,INS]
@e <- @e[INS,INS]
@fe <--  length(INS) - 1
@u <- 1/diag(solve(@h + @e))
@v <- 1/diag(solve(@e))
@f <- ((@u - @v)/@fh)/(@v/@fe)
delete(@u,@v,@h,@e)
@deleted <- grade(@f)[1]
INS <- if(isscalar(INS)){
	0
}else{
	INS[-@deleted]
}
delete(@deleted)
OUTS <- if(INS[1] == 0){
	run(@p)
}else{
	run(@p)[-INS]
}
delete(@p)
structure(f:delete(@f,return:T),df:vector(delete(@fh,return:T),\
	delete(@fe,return:T)),ins:INS,outs:OUTS)
%backstep%


===> mulvarhelp <===
mulvarhelp   MACRO
) Macro to get help on macros in mulvar.mac
# usage $S(topic1 [, topic2 ...] [help keywords])
if(!ismacro(_gethelp)){
	getmacros(_gethelp,quiet:T,printname:F)
}
___HELPFILE_ <- "mulvar.mac"
___INDEXTOP_ <- "mulvar_index"
___MACRO_ <- "$S"
_gethelp($0)
%mulvarhelp%

_E_N_D_O_F_M_A_C_R_O_S_ Marker for end of macros and data

Help file for mulvar.mac for MacAnova
(C) 2004 by Gary W. Oehlert and Christopher Bingham
Updated 041024 CB

!!!! Starting marker for message of the day
!!!! Ending marker for message of the day

???? Starting marker for list of up to 32 comma/newline separated keys
Classification
Covariance
Descriptive statistics
Diagnostics
Distances
Discrimination
Factor Analysis
Hotelling Tsq
Hypothesis tests
Iteration
MANOVA
Plotting
Random numbers
Stepwise
Subset Selection
Transformations
???? Ending marker for keys

====backstep()#classification,discrimination,stepwise,subset selection
%%%%
backstep(H,E,fh,fe), H and E symmetric REAL matrices of the same size
  with no MISSING values
%%%%
NOTE: This macro is OBSOLETE and is retained only for backward
compatibility because it was in file MacAnova.mac in earlier versions of
MacAnova.  For doing stepwise variable selection in discriminant
analysis you should use newer macros dastepsetup(), daentervar(),
daremovevar(), dastepstatus() and dasteplook().

@@@@when_used#When backstep() is used
You can use macro backstep() to perform a variable elimination step in
backwards stepwise variable selection in linear discriminant analysis.

Macro backstep() is intended to be used after you have used manova() to
compute hypothesis and error matrices H and E, with fh and fe degrees of
freedom respectively.

@@@@status_information#Status information
Status information about the variables currently "in" and "out" is
maintained in integer vectors INS and OUTS containing numbers of
variables currently included and currently excluded.  When no variables
are "in", INS = NULL (INS = 0 means the same thing); when all variables
are "in", OUTS = NULL.  INS must be initialized, usually to run(p),
before backstep() can be used.

@@@@usage#Usage
backstep(H,E,fh,fe) computes F-to-delete for all variables currently
included in vector INS.  It then updates INS by removing the index of
the variable with the smallest F-to-delete, setting INS to 0 if no
variable is left "in".  OUTS is computed as run(p) when there are no
variables "in" and as run(p)[-INS] otherwise.

@@@@value_returned#Value returned
The value returned is structure(f:F_to_delete, df:vector(fh,fe-k+1),
ins:INS,outs:OUTS), where F_to_delete is the vector of F-to-delete
statistics and k = length(INS) before deletion.  INS and OUTS are
copies of the updated INS and OUTS vectors.  The F-to-delete statistics
have nominal degrees of freedom fh an fe - k + 1.

The variable that is deleted is the one with the smallest F-to-delete
statistic.  If this is large enough, you may want to reenter the
variable using forstep().  See forstep().

@@@@example#Example
You can initialize things by
  Cmd> manova("y = groups",silent:T) # response matrix y, factor groups

  Cmd> H <- matrix(SS[2,,]); E <- matrix(SS[3,,])

  Cmd> fh <- DF[2]; fe <- DF[3]

  Cmd> INS <- run(nrows(H)) # all variables "in"
@@@@______

Macro forstep() is available for doing a forward step (variable
inclusion).  One difference between backstep() and forstep() is that
backstep() determines the variable to eliminate, and then updates INS
and OUTS; you must tell forstep() which variable to include.  See
forstep() for details.  See also compf() which computes F-to-enter for
variables not in INS.  compf() does not compute F-to-delete.

Both forstep() and compf() are OBSOLETE and are retained only for
backward compatibility.

@@@@see_also#Cross references
See also manova(), daentervar(), daremovevar(), dastepsetup(),
dastepstatus(), dastepselect(), and dasteplook().
@@@@______

====chiqqplot()#plotting,diagnostics
%%%%
chiqqplot(y [,sqrt:T] [,graphics keywords phrases]), REAL matrix y
chiqqplot(dsq, df [,sqrt:T] [, graphic keyword phrases]), REAL vector
  or matrix dsq of squared distances, REAL scalar df > 0
%%%%
@@@@usage#Usage
You use chiqqplot() to make a chisquare or sqrt(chisquare) quantile-
quantile (Q-Q) plot of ordered generalized distances against chisquare
or sqrt(chisquare) quantiles.

chiqqplot(y [, graphics keyword phrases]), where y is a REAL matrix
with no MISSING elements, plots the ordered values of generalized
distances against x[i] = invchi((i - 1/2)/n, df), where df = p =
ncols(y).  Common graphics keywords are 'xlab', 'ylab', 'title' and
'symbols'.  'symbols' may be used as with chplot(), except that
'symbols:0' labels the points by row number.

The generalized distance of the data vector y[i,]' for a case from ybar,
the sample mean vector = tabs(y,mean:T) is
       dsq[i] = (y[i,]' - ybar)' %*% solve(s) %*% (y[i,]' - ybar)
with s the sample variance-covariance matrix with divisor n-1.

@@@@assessment_of_normality#Assessment of normality
When the rows of y are a random sample from a multivariate normal
distribution, d[i] is distributed approximately as chi-squared on p
degrees of freedom and the plot should approximate a straight line with
slope 1.

Very commonly, the "data" matrix y = RESIDUALS, the matrix of residuals
computed by manova().  The Q-Q plot then helps assess the multivariate
normality of the errors.

@@@@sqrt_keyword#Keyword sqrt
chiqqplot(y, sqrt:T ...) does the same, except the ordered values of
sqrt(dsq[i]) are plotted against sqrt(invchi((i - 1/2)/n, df)).
Again the plot should be linear when y is multivariate normal.  This is
often an easier plot to assess linearity than a chi-squared QQ-plot.

@@@@distance_argument#Distance argument
chiqqplot(dsq, df, [,sqrt:T] ...) does the same except the ordered
columns of REAL vector or matrix dsq or sqrt(dsq) are plotted against
quantiles of chisquare or sqrt(chisquare) with df degrees of freedom.

When dsq is a matrix, the columns are separately ordered and plotted
with different symbols.  With symbol:0, the symbols are the row numbers
when ncols(dsq) = 1 and are the column numbers otherwise.

The usage chiqqplot(y [,sqrt:T] ...) is equivalent to
chiqqplot(distcomp(y), ncols(y) [,sqrt:T] ...).

@@@@see_also#Cross reference
See also distcomp() and tabs().

@@@@______

====compf()#
%%%%
f <- compf(h,e,fh,fe), REAL symmetric matrices h and e, positive
  integers fh and fe; integer vector INS must be defined
%%%%
NOTE: This macro is OBSOLETE and is retained only for backward
compatibility because it was in file MacAnova.mac in earlier versions of
MacAnova.  For doing stepwise variable selection in discriminant
analysis you should use newer macros dastepsetup(), daentervar(),
daremovevar(), dastepselect(), dastepstatus() and dasteplook().

@@@@when_used#When compf() is used
Macro compf() computes Fs-to enter at any stage in stepwise variable
selection in linear discriminant analysis.

You can use compf() after manova() has computed hypothesis and error
matrices H and E, with fh and fe degrees of freedom respectively.

@@@@setup#Setup
You must define Integer vector INS, containing the variable numbers
that are "in".  INS = 0 means no variables are in and the Fs-to-enter
are simply the separate ANOVA Fs.  When INS != 0, the Fs-to-enter are
the analysis of covariance Fs for each "out" variable, with the "in"
variables being used as covariates.

@@@@usage#Usage
compf(H,E,fh,fe) returns structure(f:f_to_enter,df:vector(fh,fe-k),
ins:INS, outs:OUTS), where OUTS is run(p) when INS = 0, is NULL when INS
contains all integers 1, ..., p and is run(p)[-INS] otherwise, where p =
ncols(H).  k is the number of variables "in".  The F-to-enter statistics
have nominal degrees of freedom fh and fe - k.

@@@@example#Example
Here is an example of starting forward stepwise variable selection.

  Cmd> manova("y = groups",silent:T)#, response matrix y, factor groups

  Cmd> H <- matrix(SS[2,,]); E <- matrix(SS[3,,])

  Cmd> fh <- DF[2]; fe <- DF[3]

  Cmd> INS <- 0 # no variables in

  Cmd> results <- compf(H,E,fh,fe)

  Cmd> j <- grade(results$f,down:T)[1] # index of largest F

  Cmd> # now continue with forstep(j,H,E,fh,fe)

@@@@see_also#Cross references
See also manova(), forstep(), backstep() daentervar(), daremovevar(),
dastepsetup(), dastepselect(), dastepstatus() and dasteplook().

@@@@______

====covar()#covariance,descriptive statistics
%%%%
covar(x), REAL matrix x with no MISSING values, returns
  structure(n:sampleSize,mean:xbar,covariance:s)
%%%%
@@@@usage#Usage
covar(x), where x is a REAL matrix with no MISSING values computes the
sample mean and variance-covariance matrix of x.  It returns
         structure(n:N, mean:xbar, covariance:s)
where N = nrows(x) is the sample size, xbar is a row vector of the
sample means of the columns, and s is the sample variance-covariance
matrix (with divisor N - 1).

Macro covar() is obsolete but is retained for backward compatibility.

@@@@comparison_with_tabs#Comparison with tabs()
Essentially the same output can be obtained from tabs(y, count:T,
covar:T, mean:T) which returns structure(mean:means, covar:s,count:n).
The components of covar() output differ from those of tabs() in three
ways:
 (a) Component names ('covariance' instead of 'covar' and 'n'
     instead of 'count')
 (b) The value of 'mean' is a row vector (1 by ncols(x)) for covar() and
     a (column) vector for tabs()
 (c) The value of covar() component 'n' is the scalar N, while tabs()
     component 'count' is rep(N,ncols(x)), with a count for each column
     of x

For both covar() and tabs() with covar:T, it is an error if y has any
MISSING elements.

@@@@see_also#Cross references
See also tabs() and groupcovar().
@@@@______

====cumpillai()#manova,hypothesis tests
%%%%
P <- cumpillai(tr,fh,fe,p [,upper:T] [,nterms:m or useF:T])
P <- cumpillai(tr,fh,fe,p,pillai:T [,upper:T] [,nterms:m or useF:T])
allP <- cumpillai(tr,fh,fe,p [,upper:T],all:T)
allF <- cumpillai(tr,fh,fe,p [,upper:T],useF:T,all:T)
  tr REAL variables, fh, fe, p integer variables > 0, integer scalar
  0 < m <= 4; non-scalar tr, fh, fe or p must have same dimensions
%%%%
@@@@introduction#Introduction
cumpillai() computes approximate values for lower and upper tail
probabilities for Pillai's trace statistic for testing MANOVA hypotheses.

If H is a p by p MANOVA hypothesis matrix with fh degrees of freedom and
E is a p byh p MANOVA error with fe degrees of freedom, the test
statistic is

  trpillai = trace(solve(H+E,H)) = sum(vals/(1 + vals) = V/(fh + fe)

where V is Pillai's statistic and vals is the vector of eigenvalues of H
relative to E.

cumpillai() is an interface to cumtrace() which is invoked with
'pillai:T' as an argument.  See cumtrace() for details on the
computational method.

@@@@usage#Usage
P <- cumpillai(trpillai, fh, fe, p, pillai:T) computes lower tail
probabilities of trpillai = trace(solve(H+E) %*% H), where p by p
matrices H and E are MANOVA hypothesis and error matrices with fh and fe
degrees of freedom.

Any nonscalar arguments must all have the same dimensions and P will be
a REAL variable with the same dimensions.

When min(p,fh) > 1, P is computed from up to four terms of a series due
to Fujikoshi whose first term is cumchi(m3*trpillai,fh*p), m3 = fh+fe,
the standard large sample chi-squared approximation.  Terms 2, 3 and 4
are O(1/m3), O(1/m3^2) and O(1/m3^3) adjustments, respectively.

When min(p,fh) = 1, exact probabilities based on the F-distribution are
used.

P <- cumpillai(trpillai, fh, fe, p, nterms:m) does the same except
exactly m terms of Fujikoshi's series are used, where 1 <= m <= 4, even
when min(fh,p) = 1.

P <- cumpillai(trpillai, fh, fe, p, useF:T) does the same except
probabilities are computed from an approximation based on the
F-distribution used in SPSS and SAS.  This is exact when min(p,fh) = 1,
but otherwise appears to be less accurate than Fujikoshi's
approximation.

Note: The default 'useF:F' differs from cumwilks() for which 'useF:T' is
the default.

@@@@upper_keyword#Keyword phrase 'upper:T'
P <- cumpillai(trpillai, fh, fe, p, upper:T [,nterms:m or useF:T])
computes upper tail probabilities.  This is usually what is wanted since
large values of tr are evidence against the null hypothesis.

@@@@all_keyword#Keyword phrase 'all:T'
See cumtrace() for information about the effect of 'all:T'.

@@@@reference#Reference
See Fujikoshi, Yasunori (1973) Asymptotic formulas for the distributions
of three statistics for multivariate linear hypothesis,  Ann. Inst.
Statist. Math. 25 423-437

@@@@see_also#Cross References
See also manova(), cumwilks(), cumtrace(), seqF().

====cumtrace()#manova,hypothesis tests
%%%%
P <- cumtrace(tr,fh,fe,p [,upper:T] [,nterms:m or useF:T] [,pillai:T])
P <- cumtrace(tr,fh,fe,p,pillai:T [,upper:T] [,nterms:m or useF:T] \
   [,pillai:T])
allP <- cumtrace(tr,fh,fe,p [,pillai:T] [,upper:T],all:T [,pillai:T])
allF <- cumtrace(tr,fh,fe,p [,pillai:T] [,upper:T],useF:T,all:T \
   [,pillai:T])
  tr REAL variables, fh, fe, p integer variables > 0, integer scalar
  0 < m <= 4; non-scalar tr, fh, fe or p must have same dimensions
%%%%
@@@@introduction#Introduction
cumtrace() computes approximate lower or upper tail probabilities for
Hotelling's and Pillai's trace statistics for testing MANOVA hypotheses.

The test statistics are

  trhotelling = trace(solve(E,H)) = T_0^2/fe = sum(vals)
  trpillai = trace(solve(H+E,H)) = V/(fh + fe) = sum(vals/(1 + vals)

where T_0^2 is Hotelling's generalized T^2, V is Pillai's statistic and H
and E are p by p MANOVA hypothesis and error matrices with fh and fe
degrees of freedom.  vals is the vector of eigenvalues of H relative to E.

By default, cumtrace() uses a series due to Fujikoshi whose first term is
the usual chi-squared approximation for the distribution of m2*T_0^2/fe, m2
= fe - p - 1 or m3*V/(fh+fe), m3 = fh + fe.  When min(p,fh) = 1, cumtrace()
computes exact probabilities based on the F-distribution.  

@@@@usage#Usage
P <- cumtrace(trhotelling, fh, fe, p) computes lower tail probabilities
of trhotelling = trace(solve(E) %*% H) where p by p matrices H and E are
MANOVA hypothesis and error matrices with fh and fe degrees of freedom.

P <- cumtrace(trpillai, fh, fe, p, pillai:T) computes lower tail
probabilities of trpillai = trace(solve(H+E) %*% H).  This is identical to
P <- cumpillai(trpillai, fh, fe, p).

Any nonscalar arguments must all have the same dimensions and P will be
a REAL variable with the same dimensions.

When min(p,fh) > 1, P is computed from up to four terms of a series due
to Fujikoshi whose first term is cumchi(mj*tr,fh*p), the standard large
sample chi-squared approximation.  Without pillai:T, tr = trhotelling
and mj = m2 = fe - p - 1.  With pillai:T, tr = trpillai and mj = m3 = fh
+ fe.  Terms 2, 3 and 4 are O(1/mj), O(1/mj^2), and O(1/mj^3)
adjustments, respectively.

When min(p,fh) = 1, exact probabilities based on the F-distribution are
returned.

P <- cumtrace(tr, fh, fe, p, nterms:m [,pillai:T]) does the same except
exactly m terms of Fujikoshi's series are used, where 1 <= m <= 4, even
when min(fh,p) = 1.

P <- cumtrace(tr, fh, fe, p, useF:T [,pillai:T]) does the same except
probabilities are computed from approximations based on the F-distribution
that are used in SPSS and SAS.  These are exact when min(p,fh) = 1, but
otherwise appear to be less accurate than Fujikoshi's series.

Note: The default for 'useF' differs from cumwilks() which uses a
F-based approximation by default.

@@@@upper_keyword#Keyword phrase 'upper:T'
P <- cumtrace(tr, fh, fe, p, upper:T [,nterms:m or useF:T] [,pillai:T])
computes upper tail probabilities.  This is usually what is wanted since
large values of trhotelling or trpillai are evidence against the null
hypothesis.

@@@@all_keyword#Keyword phrase 'all:T'
allP <- cumtrace(tr, fh, fe, p, all:T [,upper:T] [,pillai:T]) sets allP
to structure(P:cumPStr, terms:termStr).  cumPStr = structure(P_1:p1,
P_2:p2, P_3:p3, P_4:p4), where p1, p2, p3, and p4 are the 1, 2, 3 and 4
term aproximations to the cumulative or tail probabilities, and termStr
= structure(Term_1:term1, Term_2:term2, Term_3:term3, Term_4:term4),
where term1, ..., term4 are the actual terms in the series.

allF <- cumtrace(tr, fh, fe, p, all:T, useF:T [,upper:T] [,pillai:T])
sets allF to structure(P:cumP, f:fstat, df1:dfnum, df2:dfdenom), where
cumP contains cumulative or tail probabilities, f contains F-statistics
and dfnum and dfdenom contain numerator and denominator degrees of
freedom.

@@@@examples#Examples
  Cmd> list(Y,groups) # response matrix and factor
  groups          REAL   50    FACTOR with 4 levels
  Y               REAL   50    4

  Cmd> manova("Y = groups",silent:T)

  Cmd> h <- matrix(SS[2,,]);fh <- DF[2]

  Cmd> e <- matrix(SS[3,,]);fe <- DF[3]

  Cmd> p <- ncols(Y); vector(fh,fe,p)
  (1)           3          46           4

  Cmd> trhotelling <- trace(solve(e,h)); trhotelling
  (1)        0.77

  Cmd> trpillai <- trace(solve(h+e,h)); trpillai
  (1)     0.48115

  Cmd> cumtrace(trhotelling,fh,fe,p,upper:T) # default
  (1)   0.0047263

  Cmd> cumtrace(trpillai,fh,fe,p,pillai:T,upper:T)
  (1)    0.013864

  Cmd> cumtrace(trhotelling,fh,fe,p,upper:T,all:T)
  component: P
    component: P_1
  (1)   0.0016117        [usual chi-squared approximaition]
    component: P_2
  (1)    0.003814
    component: P_3
  (1)   0.0047398
    component: P_4
  (1)   0.0047263        [value returned by default usage above]
  component: terms
    component: Term_1
  (1)   0.0016117
    component: Term_2
  (1)   0.0022023
    component: Term_3
  (1)  0.00092578
    component: Term_4
  (1) -1.3487e-05

  Cmd> cumtrace(trpillai,fh,fe,p,upper:T,pillai:T,nterms:1)
  (1)    0.023213        [usual chi-squared approximation]

  Cmd> cumtrace(trhotelling,fh,fe,p,upper:T,useF:T) # F-approximation
  (1)   0.0031321        [F-approximation probability]

  Cmd> cumtrace(trpillai,fh,fe,p,upper:T,useF:T,all:T)
  component: P
  (1)    0.081057        [F-approximation probability]
  component: f
  (1)      1.6707
  component: df1
  (1)          12
  component: df2
  (1)         125

From a simulation using 100,000 samples, the P-value for trhotelling is
estimated as 0.00451 and for trpillai is estimated to be 0.01431.  These
are quite close to the values computed by the default usage of
cumtrace() and less close to the values computed with useF:T


@@@@reference#Reference
See Fujikoshi, Yasunori (1973) Asymptotic formulas for the distributions
of three statistics for multivariate linear hypothesis,
Ann. Inst. Statist. Math. 25 423-437

@@@@see_also#Cross References
See also manova(), cumwilks(), cumpillai(), seqF().

====cumwilks()#manova,hypothesis tests
%%%%
cumP <- cumwilks(lambdastar, fh, fe, p [,upper:T] [,useF:F]\
  [,nterms:m])
allF <- cumwilks(lambdastar, fh, fe, p, all:T, [,upper:T])
allP <- cumwilks(lambdastar, fh, fe, p, all:T, useF:F [,upper:T])
  lambdastar REAL variables, fh, fe, p integer variables > 0, integer
  scalar m > 0; all non-scalar lambdastar, fh, fe or p must have same
  dimensions
%%%%
@@@@introduction#Introduction
cumwilks() computes exact or approximate values for the cumulative
distribution of the MANOVA test statistic lambdastar = lambda^(2/N) =
det(E)/det(H+E) = 1/prod(1+vals).  lambda is the LR statistic and vals
is the vector of eigenvalues of a hypothesis matrix H relative to an
error matrix E.

The default usage uses a formula due to Rao which involves F and is
exact when min(p,fh) <= 2.  Alternatively, you can use 1, 2, 3 or 4
terms of a series due to Rao and Box, the first term of which is the
standard large sample chi-squared approximation to the distribution of
-(fe-(p-fh+1)/2)*log(lambdastar).

@@@@usage#Usage
cumP <- cumwilks(lambdastar, fh, fe, p) computes lower tail
probabilities of lambdastar = det(E)/det(H+E) = 1/prod(1 + vals).  H and
E are p by p MANOVA hypothesis and error matrices with fh and fe degrees
of freedom, and vals are eigenvalues of H relative to E.  Any nonscalar
arguments must all have the same dimensions and cumP will be a REAL
variable with the same dimensions.

cumP is computed from Rao's formula approximating the distribution of
(1 - lambdastar^(1/t))/lambdastar^(1/t) by that of a multiple of F,
where t is computed from fh, fe and p.  This is exact when s = min(p,fh)
<= 2.

cumP <- cumwilks(lambdastar, fh, fe, p, useF:F) does the same except,
when min(fh,p) > 2, up to 4 terms of a series due to Rao and Box are
used.  Term 1 is the usual chi-squared(p*fh) approximation for the
distribution of -m1*log(lambdastar) with m1 = fe - (p - fh + 1)/2.
Terms 2, 3 and 4 are O(1/m1^2), O(1/m1^4) and O(1/m1^4) adjustments to
the chi-squared probability.  The exact probability is computed when
min(fh,p) <= 2.

cumP <- cumwilks(lambdastar, fh, fe, p, nterms:m) uses m terms of the
Rao-Box series, where 1 <= m <= 4, even when min(fh,p) <= 2.  'useF:F'
is implied by 'nterms:m'.

@@@@upper_keyword
tailP <- cumwilks(lambdastar,fh,fe,p,upper:T ...) computes upper tail
probabilities.  This is usually not what is wanted since small values of
lambdastar are evidence against the null hyothesis.

@@@@all_keyword
allF <- cumwilks(lambdastar, fh, fe, p, all:T [,upper:T]) sets allF to
structure(P:cumP,f:fstat,df1:dfnum,df2:dfdenom), where cumP contains
lower or upper tail probabilities, fstat contains the F-statistics and
dfnum and dfdenom contain the numerator and denominator degrees of
freedom.

allP <- cumwilks(lambdastar, fh, fe, p, all:T [,upper:T]) sets allP to
structure(P:cumPStr, terms:termStr), where cumPStr =
structure(P_1:cumP1, P_2:cumP2, P_3:cumP3, P_4:cumP4) and termStr =
structure(Term_1:t1, Term_2:t2, Term_3:t3, Term_4:t4).  cumP1,..., cumP4
are 1, ...,4 term approximations and t1,...,t4 are the terms themselves
so that P_1 = t1, P_2 = t1 + t2, P_3 = t1 + t2 + t3 and P_4 = t1 + t2 +
t3 + t4.

@@@@examples#Examples
  Cmd> list(Y,groups) # response matrix and factor
  groups          REAL   50    FACTOR with 4 levels
  Y               REAL   50    4

  Cmd> manova("Y=groups",silent:T)

  Cmd> fh <- DF[2]; fe <- DF[3], p <- ncols(Y)

  Cmd> vector(fh,fe,p)
  (1)           3          46           4

  Cmd> h <- matrix(SS[2,,]); e <- matrix(SS[3,,])

  Cmd> lambdastar <- det(e)/det(e + h);lambdastar
  (1)      0.54836

  Cmd> cumwilks(lambdastar,fh,fe,p) # Rao's F-approximation
  (1)    0.0077151

  Cmd> cumwilks(lambdastar,fh,fe,p,all:T) # same with details
  component: P
  (1)   0.0077151
  component: f
  (1)      2.4232
  component: df1
  (1)          12
  component: df2
  (1)      114.06

  Cmd> cumwilks(lambdastar,fh,fe,p,useF:F) # Rao/Box series
  (1)    0.0077152

  Cmd> cumwilks(lambdastar,fh,fe,p,nterms:1) # use 1st term (chisq)
  (1)    0.0076322

  Cmd> cumwilks(lambdastar,fh,fe,p,nterms:2) # use two terms
  (1)    0.0077147

  Cmd> cumwilks(lambdastar,fh,fe,p,useF:F,all:T) # find all terms
  component: P
    component: P_1
  (1)   0.0076322           [chi-squared(12) approximation]
    component: P_2
  (1)   0.0077147           [2 term value]
    component: P_3
  (1)   0.0077152           [3 term value]
    component: P_4
  (1)   0.0077152           [4 term value]
  component: terms
    component: Term_1
  (1)   0.0076322           [O(1) term]
    component: Term_2
  (1)  8.2563e-05           [O(1/m1^2) term]
    component: Term_3
  (1)   4.444e-07           [O(1/m1^4) term]
    component: Term_4
  (1)  2.1917e-09           [O(1/m1^6) term]


@@@@reference#Reference
See GEP Box, A general distribution theory for a class of likelihood
criteria, Biometrika 36 (1949) 317-346 and CR Rao, Tests of significance
in multivariate analysis, Biometrika 35 (1948) 58-79.

@@@@see_also#Cross References
See also manova(), cumtrace(), cumpillai(), seqF(), releigen().

====daentervar()#classification,discrimination,stepwise,subset selection
%%%%
daentervar(var1 [,var2 ...] [,silent:T]), var1, var2 ... names
  or numbers of variables to be entered
%%%%
@@@@usage#Usage
daentervar(j), where j is a positive integer, enters dependent variable
j as part of a stepwise dependent variable selection, usually as one
stage in stepwise discriminant analysis.  This is what is sometimes
called a "forward" step.  It is an error if variable j is already an
"in" variable.

daentervar() updates variable _DASTEPSTATE which encapsulates the
current state of the variable selection process and prints a report with
the new values of F-to-enter or F-to-remove, and their P-values.  See
topic '_DASTEPSTATE' and dastepsetup() for more details.  It returns a
copy of the updated _DASTEPSTATE as an invisible variable that can be
assigned but is not normally printed.

You normally choose the variable to enter as the variable with the
largest value of F-to-enter as printed by macro dastepsetup(),
dastepstatus() or a preceding use of daentervar()

@@@@entering_named_variable#Entering named variable
daentervar(varname), where varname is the quoted or unquoted name of
a variable in the model, does the same.  Thus if the original data
matrix had column labels "SepLen", "SepWid", "PetLen" and "PetWid",
either daentervar(SepWid) or daentervar("SepWid") would be equivalent
to daentervar(2).  If varname ends in a number,

@@@@entering_by_numbered_name
daentervar(x3), say, is an alternate way to specify variable 3.  Any
variable name like x3 or var17 that doesn't start with '@' and ends with
a positive integer will work.  If the numbered name is the name
of a variable in the model, that variable will be entered, even if the
number is wrong.

@@@@entering_several_variables#Entering several variables
daentervar(j1, j2, ...) and daentervar(varname1, varname2, ...)
successively enter several variables, printing a report after each is
entered.  The value returned is _DASTEPSTATE after all the variables
have been entered.

@@@@silent_keyword#Keyword silent
daentervar(j1 [,j2, ...], silent:T) and daentervar(varname1 [,varname2,
...], silent:T) do the same, except no report is printed.  This would
normally be followed by dastepstatus() to print a report after all the
variables have been entered.

@@@@see_also#Cross references
See also dastepstatus(), daremovevar()  and dasteplook().
@@@@______

====daremovevar()#classification,discrimination,stepwise,subset selection
%%%%
daremovevar(var1 [,var2 ...] [,silent:T]), var1, var2 ... names or
  numbers of variables to be removed

%%%%
@@@@usage#Usage
daremovevar(j), where j is a positive integer, removes dependent
variable j as part of a stepwise dependent variable selection, usually
as one stage in stepwise discriminant analysis.  This is what is
sometimes called a "backward" step.  It is an error if variable j is
not an "in" variable (is already an "out" variable).

daremovevar() updates variable _DASTEPSTATE which encapsulates the
current state of the variable selection process and prints a report with
the new values of F-to-remove or F-to-remove, and their P-values.  See
topic '_DASTEPSTATE' and dastepsetup() for more details.  It returns a
copy of the updated _DASTEPSTATE as an invisible variable that can be
assigned but is not normally printed.

You normally choose which variable to remove as the variable with the
smallest value of F-to-remove as printed by macro dastepsetup(),
dastepstatus() or a preceding use of daremovevar().

@@@@removing_named_variable#Removing named variable
daremovevar(varname), where varname is the quoted or unquoted name of a
variable in the model, does the same.  Thus if the original data matrix
had column labels "SepLen", "SepWid", "PetLen" and "PetWid", either
daremovevar(SepWid) or daremovevar("SepWid") would be equivalent to
daremovevar(2).

@@@@removing_by_numbered_name
daremovevar(x3), say, is an alternate way to specify variable 3.  Any
variable name like x3 or var17 that doesn't start with '@' and ends with
a positive integer will work.  If the numbered name is the name of a
variable in the model, that variable will be removed, even if the number
is wrong.

@@@@removing_several_variables#Removing several variables
daremovevar(j1, j2, ...) and daremovevar(varname1, varname2, ...)
successively remove several variables, printing a report after each is
removed.  The value returned is _DASTEPSTATE after all the variables
have been removed

@@@@silent_keyword#Keyword silent
daremovevar(j1 [,j2, ...], silent:T) and daremovevar(varname1
[,varname2,...], silent:T) do the same, except no report is printed.
This would normally be followed by dastepstatus() to print a report
after all the variables have been removed.

@@@@see_also#Cross references
See also dastepstatus(), daentervar()  and dasteplook().
@@@@______

====dascreen()#classification,discrimination,subset selection
%%%%
result <- dascreen(Model [,mbest:m] [, method:Method] [,penalty:pen]\
  [,forced:varlist] [,maxvar:maxv] [,silent:T]), CHARACTER scalars Model
  and Method, positive integers m and maxv, REAL scalar pen >= 0, vector
  of distinct positive integer varlist
%%%%
@@@@introduction#Introduction
dascreen() examines all possible subsets of the columns of a data matrix
that might be used as variables in linear discrimination or
classification.  It finds the best subset or subsets according to a
variant of the AIC criterion.  Other criteria may be added in future
versions.  You can force every subset examined to contain specified
variables.

dascreen() works by brute force, examining all eligible non empty
subsets of variables on the LHS of the model.  The amount of time goes
up by about a factor of 2 for every additional response variable.  On a
reasonably fast computer, dascreen() can easily handle up to 15
variables (the default limit) in addition to any variables forced in by
keyword 'forced' (see below).  You can use keyword 'maxvar' to increase
this limit.

@@@@usage#Usage
result <- dascreen(Model) sets result to structure(subsets:Sets,
criterion:Crit).  Sets is a p1 by m matrix of integers, with Sets[,j]
containing the variable numbers (possible padded with 0's) of the j-th
"best" subset, and p1 is the size of the largest set in Sets.
Crit is a REAL vector with Crit[i] = value of the ordering criterion for
the subset encoded in Set[,j].  The default value for m = max(2^p-1, 5)
where p is the number of response variables in the model.

Model is a CHARACTER scalar or quoted string of the form "y = a", where
y is a N by p data matrix with no MISSING values and a is a factor of
length N specifying the groups to be discriminated from each other.  In
fact, dascreen() will work with any Model including at least one
non-constant term and only one error term.

result <- dascreen(Model, mbest:m) does the same, except Sets is a
matrix with m' columns numbers specifying the m' "best" subsets and Crit
is a vector of length m', containing the criterion values for each
subset, with m' = min(m,p^2-1). 'mbest:5' corresponds to the default
usage.  With mbest:1, component Sets will be a vector rather than a
matrix.

The default criterion is AIC = -N Log(det(H + E)/det(E) + pen*nparams,
where nparams = q*fh + q*(q+1)/2, q = size of subset, fh = hypothesis
degrees of freedom and pen = 2 or value of keyword 'penalty' (see
below).  The CONSTANT term, if any, is not included in the parameter
count.

@@@@forced_keyword#Keyword 'forced'
When you have reason to believe that a particular set of variables
should be used in discrimination, you can use keyword phrase
'forced:varlist' as an argument to force them to be included in every
subset considered.  varlist must be a vector of distinct positive
integers <= p, and will be included in every column of Sets.

@@@@penalty_keyword#Keyword 'penalty'
You can modify the formula for AIC by including keyword phrase
'penalty:pen' as an argument, where pen >= 0 is a scalar.  pen will be
used instead of 2 as multiplier of the number of parameters in computing
AIC.  Larger values of pen tend to reduce the size of subsets selected;
smaller values tend to reduce subset size.

@@@@silent_keyword#Keyword 'silent'
When 'silent:T' is an argument, any warning messages are suppressed.

@@@@maxvar_keyword#Keyword 'maxvar'
If the number of variables is more than the default limit (15 + number
of variables forced in), you can increase the limit by including
'maxvar:M' as an argument.  If the number of variables is large enough
this can result in prohibitively long computing time.

@@@@example#Example
  Cmd> list(y,layer)
  layer           REAL   75    1     FACTOR with 3 levels (labels)
  y               REAL   75    12    (labels)

  Cmd> dascreen("y=layer")
  component: subsets
            Set 1       Set 2       Set 3       Set 4       Set 5
  (1)           1           1           1           1           1
  (2)           2           6           2           2           6
  (3)           6          10           6           3           7
  (4)          10           0           8           6          10
  (5)           0           0          10          10           0
  component: criterion
       Set 1       Set 2       Set 3       Set 4       Set 5
      -152.6     -148.23     -147.68     -146.12     -145.99

  Cmd>  dascreen("y=layer",mbest:3,forced:vector(1,2))
  component: subsets
            Set 1       Set 2       Set 3
  (1)           1           1           1
  (2)           2           2           2
  (3)           6           6           3
  (4)          10           8           6
  (5)           0          10          10
  component: criterion
         Set 1       Set 2       Set 3
        -152.6     -147.68     -146.12

  Cmd>  dascreen("y=layer",penalty:3,forced:vector(1,2))
  component: subsets
            Set 1       Set 2       Set 3       Set 4       Set 5
  (1)           1           1           1           1           1
  (2)           2           2           2           2           2
  (3)           6          10           9           8           6
  (4)          10           0          10          10           8
  (5)           0           0           0           0          10
  component: criterion
         Set 1       Set 2       Set 3       Set 4       Set 5
        -134.6     -127.81     -126.78     -124.67     -122.68

@@@@see_also#Cross references
See also dastepselect(), screen().

====dasteplook()#classification,discrimination,stepwise,subset selection
%%%%
dasteplook(name1, name2, ...), names selected from  'model', 'hpluse',
  'e', 'in', 'F', 'fh', 'fe', 'history' and 'varnames'
dasteplook(all)
%%%%
@@@@usage#Usage
dasteplook(name), where name is the name of a component of variable
_DASTEPSTATE, that is one of 'model', 'hpluse', 'e', 'in', 'F', 'fh',
'fe', and 'history', returns that component of _DASTEPSTATE.  See topic
'_DASTEPSTATE' for details on these components

dasteplook(varnames) returns the names of the variables as a CHARACTER
vector.

name can either be quoted, as in dasteplook("hpluse"), or unquoted, as
in dasteplook(hpluse).

@@@@viewing_several_components#Viewing several components
dasteplook(name1, name2, ...), where the names are quoted or unquoted
_DASTEPSTATE component names or 'varnames' returns a structure
consisting of the requested items.

dasteplook(all) or dasteplook("all") return all of _DASTEPSTATE as
value.

@@@@purpose#Purpose of dasteplook()
dasteplook() is designed for use in a macro such as dastepselect() which
controls the use of macros dastepsetup(), daentervar(), daremovevar() to
carry out an entire stepwise dependent variable selection.

@@@@see_also#Cross references
See also topics dastepsetup(), daentervar(), daremovevar() and
dastepstatus().
@@@@______

====dastepselect()#classification,discrimination,stepwise,subset selection
%%%%
Ins <- dastepselect(Model, alpha [,updown:T] [,allin:T] [,bonferroni:F]\
  [,quiet:T  or  silent:T]), Model a CHARACTER scalar specifying a model
  for manova(), alpha > 0 a REAL scalar between 0 and .5 or between 1 and
  50
%%%%
@@@@introduction#Introduction
dastep() performs complete forward, backward or mixed direction stepwise
variable selection for linear discriminant analysis.  It provides a
"black box" interface to dastepsetup(), daentervar() and daremovevar().

@@@@usage#Usage
Ins <- dastepselect(Model, alpha) uses a stepwise procedure to try to
find subset of variables which can be used in place of a complete set of
variables in discriminant analysis.

After each successful step (variable entered or removed), the variable
name, the action taken, and a P-value are printed.  At completion, the
selected subset is printed and variable numbers are returned as an
"invisible" variable.

Model is a CHARACTER scalar or quoted string specifying a GLM model
suitable for manova().  Its usual form is "y = a", where y is an n by p
REAL matrix and a is a factor with length(a) = n defining the groups to
be discriminated among.  The columns of the dependent variable (y) are
the complete set of variables.

alpha is a positive scalar between 0 and .5 or between 1 and 50
specifying a significance level.

Ins becomes a vector of distinct positive integers (column numbers
selected) or NULL (when no columns are selected).

By default the subset is selected by forward selection with a stopping
rule based on Bonferronized P-values of analysis of covariance
F-statistics in analyses of covariance of each "out" variable (variable
not yet entered) with the "in" variables (variables already entered) as
covariates.  When there are K "out" variables, the Bonferronized P-value
is K*P, where P is the ordinary P-value.  Forward selection stops when
K*P > alpha for all "out" variables, or when all variables are "in".

@@@@backward_selection#Backward selection
Ins <- dastepselect(Model, alpha, allin:T) does the same, except
backward selection is used to eleminate one variable at a time until the
largest Bonferronized P-value > alpha or all variables are "out".  When
there are K "out" variables", the Bonferronized P-value for an "in"
variable is (K+1)*P where P is the P-value for the F-statistic in an
analysis of covariance of the variable with the "out" variables as
covariates.

@@@@updown_selection#Up and down selection
Ins <- dastepselect(Model, alpha, updown:T [,allin:T]) does the same
except both forward and backward steps can occur.  After every forward
step, a check is made to see if a backward step involving another
variable would eliminate that variable.  And after every unsuccessful
attempt at a backward step, a forward step is attempted to try to
include another variable.

The process stops when both the forward and backwards stopping
conditions are satisfied.

@@@@bonferroni_keyword#Keyword phrase 'bonferroni:F'
Ins <- dastepselect(Model, alpha, bonferroni:F ... ) does the same,
except P-values are not Bonferronized.  That is forward steps succeed
when P <= alpha and backward steps succeed wgeb P > alpha.  This almost
always results in larger subsets being selected.

@@@@quiet_and_silent_keywords#Keyword phrases 'quiet:T' and 'silent:T'
When 'quiet:T' or 'silent:T' is an additional argument, printing of the
steps taken and the final subset is suppressed.  The selected variable
numbers are returned in a "visible" vector.  'silent:T' also suppresses
any warning messages.

@@@@example#Example
  Cmd> list(y,layer)
  layer           REAL   75    1     FACTOR with 3 levels (labels)
  y               REAL   75    12    (labels)

  Cmd> print(paste(getlabels(y,2))) # column labels
  La Eu Co Fe Ta Sc Hf Yb Cr Th Lu Au

  Cmd> ins <- dastepselect("y=layer",.05) # forward selection
  Entering Th. Bonferronized P-value-to-enter = 1.1819e-10
  Entering La. Bonferronized P-value-to-enter = 0
  Entering Sc. Bonferronized P-value-to-enter = 0.0024604
  Entering Eu. Bonferronized P-value-to-enter = 0.0048223
  Smallest Bonferronized P-value-to-enter = 0.13078 > 0.05
  Variables selected: La, Eu, Sc, Th

  Cmd> ins
             La           Eu           Sc           Th
              1            2            6           10

  Cmd> dastepselect("y=layer",.05,allin:T,quiet:T) # backward selection
             La           Eu           Sc           Th
              1            2            6           10

Without 'quiet:T', this last would have reported the removal of 8
variables.

@@@@see_also#Cross references
See also dastepsetup(), daentervar(), daremovevar(), dastepstatus(),
'variables:"invisible_variables"'

====dastepsetup()#classification,discrimination,stepwise,subset selection
%%%%
dastepsetup([Model] [,allin:T or in:ins] [,silent:T]), Model a
  CHARACTER scalar glm model, usually of the form "y = groups",
  ins either a LOGICAL vector or vector of unique integers > 0
%%%%
@@@@purpose#Purpose of dastepsetup()
You use macro dastepsetup() at the start of a forward or backward
stepwise selection of dependent variables in a discriminant analysis or
more generally in a multivariate linear model.

What it actually does is create and initialize invisible variable
_DASTEPSTATE which encapsulates information on which dependent
variables are "in" and which are "out" at any stage of the variable
selection process. See topic '_DASTEPSTATE'.

@@@@linear_discrimination#Linear discrimination
The most common use is in stepwise linear discriminant analysis where
you are trying to select a subset of reponse variables that effectively
discriminate among two or more groups.  It can also be used in any
linear model when you are trying to select a subset of reponse
variables that are responsible for any violation of the overall null
hypothesis H0: all model coefficients except constant term are 0.

@@@@usage#Usage
dastepsetup(Model), where Model is a CHARACTER scalar specifying a GLM
model, initializes _DASTEPSTATE so that no variables are "in" and all
are "out".  This is appropriate at the start of forward stepwise
dependent variable selection.  In linear discrimination analysis, Model
has the form "y = groups", where groups is a factor defining the groups
to be discriminated.

@@@@printed_output#Printed output
A report of the current status is printed.  This includes all the
F-to-enter statistics and their P-values.

@@@@value#Value returned
A copy of _DASTEPSTATE is returned as an "invisible" variable which can
be assigned but is not automatically printed.

@@@@silent_keyword#Keyword silent
dastepsetup(Model, silent:T) does the same, except the printed report
is suppressed.

@@@@default_model#Default model
dastepsetup([,silent:T]) does the same, except variable STRMODEL,
usually the most recent GLM model used, is taken as Model.

@@@@allin_keyword#Keyword allin
dastepsetup([Model], allin:T [,silent:T]) does the same, except that
all response variables are "in" and no variables are "out".  Component
'history' of _DASTEPSTATE is initialized to run(p), where p is the
number of variables.

@@@@in_keyword#Keyword in
dastepsetup([Model], in:Ins [,silent:T]) does the same, except only the
variables specified by vector Ins are "in".  Ins can be either LOGICAL
or REAL.

When Ins = vector(j1,j2, ...) is REAL, its elements must be unique
integers between 1 and p = number of variables.

When Ins is LOGICAL, length(Ins) must be p and variables j1, j2, ...
will be "in", where Ins[j1], Ins[j2] ... are the only elements with value
True (T).  Thus ins:Ins has the same effect as ins:run(p)[Ins].

Component 'history' is initialized to vector(j1,j2,...).

@@@@what_you_do_next#What you do next
After dastepsetup(), your next step is to use daentervar() to enter a
new variable or daremovevar() to remove a variable.  The choice of which
variable to enter or remove is usually made on the basis of the
F-to-enter and/or F-to-remove statistics in the printed report.

@@@@see_also#Cross references
See also topics daentervar(), daremovevar(), dastepstatus() and
dasteplook().
@@@@______

====_DASTEPSTATE%#classification,discrimination,stepwise,subset selection
%%%%
Invisible variable _DASTEPSTATE encapsulates the current state of a
process of stepwise variable selection in linear discriminant analysis.

You use macro dastepsetup() to initialize _DASTEPSTATE, macros
daentervar() and daremovevar()e to update it, macro dastepstatus() to
report F-to-enter and F-to-remove statistics in _DASTEPSTATE and macro
dasteplook() to extract information from _DASTEPSTATE.
%%%%
@@@@description#Description
_DASTEPSTATE is an invisible variable which encapsulates the current
state of a process of stepwise response variable selection for a
multivariate linear model.  The model is most commonly of the form "y =
groups", groups a factor, and the stepwise process corresponds to
stepwise linear discriminant analysis.  You normally don't need to be
concerned with _DASTEPSTATE itself, since all interaction with it can
be done using macros.

In the following, p = number of variables, E = p by p error matrix from
manova() and H = p by p hypothesis matrix.  For a model of the form "y
= groups", H = matrix(SS[2,,]) and E = matrix(SS[3,,]).  H excludes
contributions from a constant term.  See topic 'models'

At a particular stage there are from 0 to p "in" variables, variables
tentatively considered important in rejecting the null hypothesis H0
associated with H; the remainder are "out" and either have not yet been
considered or are tentatively considered unimportant in rejecting H0.

At each stage there is an F-to-enter statistics for each "out"
variable, if any.  This is the F-statistic in an analysis of covariance
of the variable with the "in" variables as covariates.  When no
variables are "in", this is just the usual ANOVA F-statistic for the
variable.

Similarly, at each stage, there is an F-to-remove statistic for each
"in" variable, if any.  This is the F-to-enter statistic for the
variable that would be computed if it were to be removed and turned
into an "out" variable.

@@@@contents#Contents of _DASTEPSTATE
_DASTEPSTATE has the form
 structure(model:glmModel, hpluse:(H+E)swept, e:Eswept, in:ins,\
     F:F_stats, fh:hypDF, fe:errDF, history:hist)

  glmModel      CHARACTER scalar specifying a GLM model, usually of the
                form "y = groups", groups a factor
  ins           LOGICAL vector of length p with ins[j] = T when variable
                j is "in"
  (H+E)swept    H+E when no variables are "in"; swp(H+E,run(p)[ins]),
                otherwise;  that is any "in" variables have been swept
  Eswept        E when no variables are "in"; swp(E,run(p)[ins])
                otherwise
  F_statistics  A REAL vector of F-statistics, with F[ins] being values
                of F-to-remove and F[!ins] being values of F-to-enter
  hypDF         Numerator d.f. of F-to-enter and F-to-remove = fh, the
                degrees of freedom associated with H
  errDF         Denominator d.f of F_statistics = fe - k - 1 for "in"
                and fe - k otherwise, where fe are the error d.f. and
                k = sum(ins) = number of "in variables"
  hist          An integer vector summarizing the path followed to
                reach the current state.  when hist[i] = j > 0, variable
                j was entered at step i; when hist[i] = -j < 0, variable
                j was removed at step i.

@@@@when_to_use#When to use dastepsetup()
You normally use macro dastepsetup() to initialize _DASTEPSTATE.
dastepsetup() uses manova() to find E, H, fh and fe.  When variables j1,
j2, ... jk are specified as being "in", hist is initialized to
vector(j1,...,jk); when no variables are "in" at the start, hist is
NULL.  dastepsetup() calls stepstatus to compute component 'F' and
optionally report on the initial status.

@@@@what_you_do_next#What you do next
You use macros daentervar() and daremovevar() to update _DASTEPSTATE by
changing an "out" variable to an "in" or vice versa.  These optionally
print a report of the new status.

You use macro dastepstatus() to compute component F and optionally print
a report of the information in _DASTEPSTATE.

You use macro dasteplook() to extract components of _DASTEPSTATE without
changing it.

Macros daentervar(), daremovevar(), dastepsetup() and dasetupstatus()
return the new value of _DASTEPSTATE as an "invisible" result which can
be assigned but is not printed.

@@@@see_also#Cross references
See also topics dastepsetup(), daentervar(), daremovevar(),
dasteplook() and dastepstatus().
@@@@______

====dastepstatus()#classification,discrimination,stepwise,subset selection
%%%%
dastepstatus([silent:T])
%%%%
@@@@usage#Usage
dastepstatus() computes and prints out the values of F-to-enter or
F-to-remove, and their P-values at the current stage of stepwise
dependent variable selection.  You can use this information to select
the next variable to be entered (the one with the largest F-to-enter)
or remove (the one with the smallest F-to remove).  Component
'F_statistics' of variable _DASTEPSTATE is updated.

@@@@value#Value returned
dastepstatus() returns as value an "invisible" copy of variable
_DASTEPSTATE.  This can be assigned but is not normally printed.
_DASTEPSTATE must have been previously initialized or updated by macros
dastepsetup(), daentervar() or daremovevar().

@@@@silent_keyword#Keyword silent
dastepstatus(silent:T) sets component 'F_statistics' of _DASTEPSTATE
and returns an invisible copy of _DASTEPSTATE but does not print the
F-statistics.
@@@@______

You ordinarily don't need to use dastepstatus() directly since macros
daentervar(), daremovevar() and dastepstatus() all use dastepstatus() to
report F-statistic values and their P-values.  An exception is when you
want to enter or remove several variables at once and want to see a
report only at the end.  Suppose, for example, that you want to add
variables 2 and 4 together, but don't want to see a report after
variable 2 is entered.

@@@@example#Example
  Cmd> daentervar(2,4,silent:T) # silently enter variables 2 and 4

  Cmd> dastepstatus() # print report

@@@@see_also#Cross references
See also topics dastepsetup(), daentervar(), daremovevar(), dasteplook()
and '_DASTEPSTATE'.
@@@@______

====davarid()#classification,discrimination,stepwise,subset selection
%%%%
Id <- davarid(name), name a CHARACTER scalar
%%%%
@@@@introduction#Introduction
davarid() is a utility used by daentervar() and daremovevar() to
identify a variable (column) in a data matrix.  As such, it is unlikely
to be useful to the ordinary user, but might be helpful in writing a
macro.  davarid() extracts information from variable _DASTEPSTATE which
must be defined previously.  See '_DASTEPSTATE' for details.

@@@@usage#Usage
Id <- davarid("name"), where name is the label of column j of
matrix _DASTEPSTATE$F, sets Id to structure(varnumber:j,varname:"name").
The argument must be a quoted string or CHARACTER scalar.

Id <- davarid("3"), say, sets Id to structure(varnumber:3,varname:"name")
where "name" is the label on column 3 of _DASTEPSTATE$F.

Id <- davarid("x3"), say, is the same as Id <- davarid("3"), unless "x3"
is a column label of _DASTEPSTATE$F in which case component varnumber is the
number of that column.

If var is an positive integer scalar with value 3, say, Id <-
davarid("var") is the same as Id <- davarid("3").

If var is a CHARACTER scalar with value "x3", say, Ic <- davarid("var")
is the same as Id <- davarid("x3").

@@@@example#Example
  Cmd> dastepsetup("iris = variety", silent:T)

  Cmd> dasteplook(varnames)
  (1) "SepLen"
  (2) "SepWid"
  (3) "PetLen"
  (4) "PetWid"

  Cmd> davarid("PetLen")
  component: varnumber
  (1)            3
  component: varname
  (1) "PetLen"

  Cmd> davarid("y3")
  component: varnumber
  (1)            3
  component: varname
  (1) "PetLen"

  Cmd> davarid("3")
  component: varnumber
  (1)            3
  component: varname
  (1) "PetLen"

@@@@see_also#Cross references
See also dastepsetup(), daentervar(), daremovevar(), dasteplook().

====discrim()#classification,discrimination
%%%%
discrim(groups, y), vector of positive integers groups, REAL matrix y
  with no MISSING values
%%%%
@@@@usage#Usage
discrim(groups, y), where groups is a factor or an integer vector, and
y is a REAL data matrix with no MISSING elements, computes the
coefficients of linear discriminant functions that can be used to
classify an observation into one of the populations specified by
argument groups.

The functions being estimated are optimal in the case when the
distribution in each population is multivariate normal and the
variance-covariance matrices are the same for all populations.

When there are g = max(groups) populations, and p = ncols(y) variables,
the value returned is structure(coefs:L, add:C) where L is a REAL p by
g matrix and C is a 1 by g row vector.

If y is a length p vector of data to be classified to one of the
populations, then f = L' %*% y + C' is the vector of discriminant
function values (scores) for the g populations.

If f[j] = max(f) is the largest element of f, then, assuming the g
populations are equally probable (each have prior probability 1/g), then
population j is the most probable population based on y.

@@@@prior_and_posterior_probabilities#Prior and posterior probabilities
If P is a length g vector such that P[j] = prior probability a randomly
selected case belongs to population j, then the estimated
posterior probability that y belongs to population k is
   P[k]*exp(f[k])/sum(P * exp(f)) =
            P[k]*exp(f[k] - f[1])/sum(P * exp(f - f[1]))
The second form is preferred since exp(f[k]) can be too large to
compute.

@@@@classifying_rows_of_matrix#Classifying rows of matrix
When Y is a m by p data matrix whose rows are to be classified, F = Y
%*% L + C is m by g matrix, with F[i,j] containing the value of the
discriminant function for population j evaluated with the data in row i
of Y.  A m by g matrix of posterior probabilities for each group and
case can be computed by
       P * exp(F  - F[,1])/((P * exp(F - F[,1])) %*% rep(1,g))

@@@@bias_in_estimating_posterior_probabilities#Bias in estimating posterior probabilities
It is well known that posterior probabilities computed for a case that
is in "training set", the data set from which a classification method
was estimated, are biased in an "optimistic" direction: The estimated
posterior probability for its actual population is biased upward.  For
this reason posterior probabilities should be estimated only for cases
that are not in the training set.  See macro jackknife() for a partial
remedy.

@@@@______

====discrimquad()#classification,discrimination
%%%%
discrimquad(groups, y), factor or vector of positive integers groups,
  REAL matrix y with nrows(y) = length(groups)
%%%%
@@@@usage#Usage
discrimquad(groups, y), where groups is a factor or an integer vector,
and y is a REAL data matrix, computes the coefficients of quadratic
discriminant functions that can be used to classify an observation into
one of the populations specified by argument groups.

It is an error if the smallest group has p or fewer members or if y has
any MISSING elements.

@@@@value_returned#Value returned
When there are g = max(groups) populations, and p = ncols(y) variables,
the value returned is structure(Q:q, L:l, addcon:c, grandmean:ybar),
where the components are as follows:
   q     structure(Q1,Q2,...Qg), each Qj a  REAL p by p matrix
   l     structure(L1,L2,...Lg), each L2 a REAL vector of length p
   c     vector(c1,c2,...cg), cj REAL scalars
   ybar  vector(ybar1,...ybarp), the vector of column means

@@@@quadratic_score#Quadratic score
When x is a vector of length p to be classified, the quadratic score
for group j is
    qs[j] = (x-ybar)' %*% q[j] %*% (x-ybar) + (x-ybar)' %*% l[j] + c[j]

The functions are optimal in the case when the distribution in each
population is multivariate normal with no assumption that the variance-
covariance matrices are the same for all populations.

@@@@prior_and_posterior_probabilities#Prior and posterior probabilities
When P = vector(P1,P2,...,Pg) is a vector of prior probabilities a
randomly selected case comes from the various populations, then the
posterior probabilities the elements of the vector
   P*exp(qs)/sum(P*exp(qs)) = P*exp(qs - qs[1])/sum(P*exp(qs - qs[1])

The latter form is usually preferable since it is possible for
exp(qs[1]) to be so large as to be uncomputable.  These probabilities
can be computed using macro probsquad().

@@@@bias_in_estimating_posterior_probabilities#Bias in estimating posterior probabilities
It is well known that posterior probabilities computed for a case that
is in "training set", the data set from which a classification method
was estimated, are biased in an "optimistic" direction: The estimated
posterior probability for its actual population is biased upward.  For
this reason posterior probabilities should be estimated only for cases
that are not in the training set.

@@@@see_also#Cross references
See also discrim() and probsquad().
@@@@______

====distcomp()#covariance,descriptive statistics
%%%%
distcomp(y [,center:cent] [,quadform:Q]), REAL matrix y with no MISSING
  values, REAL row or column vector cent, square positive definite REAL
  matrix Q
%%%%
@@@@usage#Usage
dsq <- distcomp(y) computes the squared generalized distances of each
row of a N by p REAL matrix y from the vector of column means.  dsq will
be a length N REAL vector of squared distances with elements

   dsq[i] = (y[i,]' - ybar)' %*% solve(S) %*% (y[i,]' - ybar)

where S and ybar are the p by p sample variance matrix and length p mean
(column) vector.

dsq <- distcomp(y, center:cent, quadform:Q) computes squared generalized
distances defined by a square p by p REAL matrix Q of each row of y from
the length p REAL vector cent.  Neither Q nor cent can have any missing
elements.  Q must be symmetric and positive definite.  In this usage

   dsq[i] = (y[i,]' - cent)' %*% solve(Q) %*% (y[i,]' - cent)

When 'center:cent' is omitted, cent = ybar is used and when 'quadform:Q'
is omitted, Q = S is used.

Note the use of the inverse of Q rather than Q itself.

@@@@normality_assessment#Assessing Multivariate Normality
When the rows of y are a random sample from a multivariate normal
distribution, the generalized distances with Q = S and cent = ybar are
approximately distributed as chi-squared on p degrees of freedom.  Thus
departure from a straight line of a chi-squared QQ plot of sort(dsq) or
a square root chi-squared QQ plot of sqrt(sort(dsq)) indicates that the
rows of y may not be multivariate normal.  You can compute the necessary
chi-square quantiles as invchi((run(n) - .5)/n, p),

This is automated by macro chiqqplot().

For other Q, the distribution may often be approximated by a multiple of
chi-squared and similarly be plotted against chi-squared quantiles.  The
degrees of freedom to use is problematical.  One useful approximation
for degrees of freedom is

  dsqdf <- 2*describe(dsq,mean:T)^2/describe(dsq,var:T)

When y = RESIDUALS, the matrix of residuals computed by manova(), this
helps assess the normality of the residuals.

@@@@see_also#Cross references
See also chiqqplot(), resvsrankits().
@@@@______

====facanal()#factor analysis
%%%%
facanal(s, m [, method:Meth] [,start:psi0] [,rotate:Rotmeth] \
  [,maxit:nmax] [,minit:nmin] [,crit:vector(nsig1,nsig2,dg)]\
  [,minimizer:M] [,gold:ngold] [,quiet:T] [,silent:T] \ [,recordwhen:d1]
  [,printwhen:d2]), REAL p.d. symmetric matrix s, integer m > 0, Meth
  one of "uls", "gls" or "mle" (default), psi0 REAL vector of starting
  uniquenesses > 0, Rotmeth one of "varimax", quartimax", "equimax",
  "none" (default), nmin >= 0, nmax > 0, ngold > 0, d1 >= 0, d2 >= 0,
  nsig1, nsig2 integers, dg REAL scalar, M one of "dfp", "bfs" (default)
  or "broyden",
%%%%
@@@@introduction#Introduction
facanal() uses an optimization macro (dfp(), bfs() or broyden()) to find
ULS, GLS or ML estimates of uniquenesses and loadings.  The optimizer
minimizes a criterion as a function of log(psi), psi = vector of
uniquenesses.

You can specify a rotation method and starting values.

facanal() is to be preferred to other macros, including ulsfactor() and
glsfactor(), for factor extraction.

By default, minimizer bfs() is used but you can specify another
minimizer.

@@@@usage#Usage
facanal(r, m, method:Meth) computes estimated uniquenesses and unrotated
estimated loadings for a orthogonal factor analysis based a p by p
correlation or covariance matix r which must be symmetric and positive
definite.

Integer m > 0 is the number of factors assumed.  It must satisfy m <
(2*p + 1 - sqrt(8*p+1))/2

Meth is must be one of "mle" or "ml" (ML or maximum likelihood method),
"uls" (ULS or unweighted least squares method) or "gls" (GLS or
generalized least squares method).  If you omit method:Meth, the default
is method:"mle".

In addition to the uniquenesses and loadings, facanal() prints the
minimized value of the criterion being minimized.  This can be used in a
goodness of fit test.

See below for the value returned by facanal().  It is "invisible" unless
'quiet:T' or 'silent:T' is an argument.

@@@@starting values#Starting values
facanal(r, m, method:Meth, start:psi0), does the same except that the
minimization iteration is started with uniquenesses psi0.  The default
starting values are psi0 = 1/diag(solve(r)).

@@@@rotation#Rotation
You control rotation of the estimated loadings using keywords 'rotate'
and 'kaiser'.

facanal(r, m, method:Meth, rotate:Rot), where Rot is one of "varimax",
"quartimax", "equimax" or "none" (the default), carries out the
indicated rotation of the factor loadings using Kaiser normalization.
It does not affect the values of the uniquenesses or criterion.

facanal(r, m, method:Meth, rotate:Rot,kaiser:F), does the same except
Kaiser normalization is not done.

@@@@optimization_control#Optimization control
There are several keywords you can use to control the actual
minimization of the criterion.  Default values are in [..]

  nmin:n1     integer n1 >= 0 = minimum number of interations performed
              [0]
  nmax:n2     integer n2 > 0 = maximum number of interations performed
              [30]
  crit:vector(nsig1,nsig2,dg)
     nsig1    target number of significant digits in log uniquenesses
              [5]
     nsig2    target number of significant digits in the criterion [8]
     dg       target upper limit for ||gradient|| [-1]
              Negative values of nsig1, nsig2 or dg are ignored.
  minimizer:M M = name of minimizer, one of "dfp" (Davidon-Fletcher-
              Powell), "bfs" (Broyden-Fletcher-Shanno) or "broyden"
              ["bfs"]
  ngold:n     integer n >= 0 = number of steps in golden mean linear
              search by minimizer [1]; ignored with minimizer:"broyden"

@@@@macros_needed#Macros needed
facanal() requires and loads automatically one of the macros ulscrit(),
glscrit() or mlcrit() to compute the objective function, gradient vector
and loading matrix.

@@@@printing_control#Printing control
There are three keywords that you can use to control what gets printed.

   quiet:T     suppresses printing of results;  the return value is
               visible
   silent:T    same as quiet:T except warning messages are also
               suppressed
   printwhen:d1  d1 >= 0 integer specifies (when d1 > 0) that partial
               results are printed at iterations d1, 2*d1, 3*d1.  These
               are the current values of x = log(uniquenesses), the F(x)
               = criterion and the gradient vector dF(x)/dx

@@@@sideeffect_variables#Side effect variables
In addition to returning them, facanal() saves the estimated
uniquenesses, rotated loadings, minimized criterion, eigenvalues and
gradient vector in side-effect variables PSI, LOADINGS, CRITERION,
EIGENVALUES and GRADIENT, respectively.

Each time they are entered, ulscrit(), glscrit() or mlcrit() save a copy
of their argument psi in invisible variable _PSIULS, _PSIGLS or _PSIML.

You can use keyword 'recordwhen' to specify that partial results are
saved in side-effect structure BFSRECORD, DFPRECORD or BROYDNRECORD on
some or all iterations:

   recordwhen:d2  d2 >= 0 integer specifies (when d2 > 0) that partial
               results are saved at iterations d2, 2*d2, 3*d2, ...
               Values of x = log(uniquenesses), F(x) = criterion and
               dF(x)/df go in components 'xvalx', 'funval' and
               'gradients' of the structure.

@@@@value_returned#Value returned
facanal() always returns as structure a value which can be assigned.
Unless 'silent:T' or 'quiet:T' is an argument, the return value is
"invisible" and won't be automatically printed.

The result has the form structure(psihat:psi,loadings:l,criterion:crit,
eigenvals:vals, gradient:g, method:Meth, rotation:Rot, iter:n, status:k)
where the values of the components are as follows:

   psi         REAL vector of estimated uniquenesses
   l           REAL p by m matrix of loadings
   crit        Minimized criterion
   vals        eigenvalues s relative to psi
   g           REAL gradient vector at optimum
   n > 0       integer = number of iterations
   k >= 0      integer specifying which (1, 2 or 3) of the convergence
               criterion signalled convergence; k = 0 means not
               converged

@@@@examples#Examples
  Cmd> facanal(cor(y),2,rotate:"varimax") # y a 50 by 5 matrix
  Convergence in 26 iterations by criterion 2
  estimated uniquenesses:
  (1)     0.39881  5.4107e-07     0.20203     0.32915     0.63907
  varimax rotated estimated loadings:
  (1,1)     0.71121     0.30883
  (2,1)     0.13876     0.99033
  (3,1)    -0.84817    -0.28032
  (4,1)    -0.80203    -0.16613
  (5,1)     0.59678   -0.069151
  minimized ml criterion:
  (1)    0.013439

  Cmd> result <- facanal(cor(y),2,rotate:"varimax",method:"uls",\
        silent:T)

  Cmd> compnames(result)
  (1) "psihat"
  (2) "loadings"
  (3) "criterion"
  (4) "eigenvals"
  (5) "gradient"
  (6) "method"
  (7) "rotation"
  (8) "iter"
  (9) "status"

  Cmd> result[vector(8,9)] # status = 0 => did not converge
  component: iter
  (1)          30
  component: status
  (1)           0

@@@@see_also#Cross references
See also ulsfactor(), glsfactor(), ulscrit(), glscrit(), mlcrit(),
minimizer().

====forstep()#factor analysis,iteration
%%%%
forstep(i,H,E,fh,fe), integer i > 0, fh > 0, fe > 0, REAL symmetric
  matrices H and E with no MISSING values
%%%%
NOTE: This macro is OBSOLETE and is retained only for backward
compatibility because it was in file MacAnova.mac in earlier versions of
MacAnova.  For doing stepwise variable selection in discriminant
analysis you should use newer macros dastepsetup(), daentervar(),
daremovevar(), dastepstatus() and dasteplook().

@@@@what_it_does#What forstep() does
Macro forstep() performs a variable inclusion step in forward stepwise
variable selection in linear discriminant analysis.

forstep() is intended to be used after you have used manova("y =
groups"), where y is a data matrix and groups is a factor, to compute
hypothesis and error matrices H = matrix(SS[2,,]) and E =
matrix(SS[3,,]), with fh = DF[2] and fe = DF[3] degrees of freedom
respectively.

@@@@status_information#Status information
Status information about the variables currently "in" and "out' is
maintained in integer vectors INS and OUTS containing numbers of
variables currently included and currently excluded.  When no variables
are "in", INS = 0; when all variables are "in", OUTS = NULL.  INS must
be initialized, usually to 0, before forstep() can be used.

@@@@usage#Usage
forstep(j,H,E,fh,fe), where j is the number of a variable not currently
"in", adds j to INS and removes it from OUTS, and then uses macro
compf() to compute F-to-enter for all variables included in the updated
INS.  The Fs-to-enter are the analysis of covariance Fs for each "out"
variable, with the "in" variables being used as covariates.  See topic
compf().  When no variables are "in", the Fs-to-enter are the ordinary
ANOVA F-statistics for each variable.

@@@@value_returned#Value returned
The value returned (which will normally be printed if not assigned) is
structure(f:F_to_enter, df:vector(fh,fe-k), ins:INS,outs:OUTS), where
F_to_enter is the vector of F-to-enter statistics, one for each
variable not in INS, INS and OUTs are copies of the status vectors INS
and OUTS.  k is the number of variables currently "in".

The F-to-enter statistics have nominal degrees of freedom fh and fe -
k.  The next variable to be entered, if any, is normally the variable
with the largest F-to-enter.  The decision to enter it is based on the
size of F-to-enter.

@@@@example#Example
You can somewhat automate the start of this process as follows:

  Cmd> manova("y = groups", silent:T) # response matrix y, factor groups

  Cmd> H <- matrix(SS[2,,]); E <- matrix(SS[3,,])

  Cmd> fh <- DF[2]; fe <-  DF[3]

  Cmd> INS <- 0; stuff <- compf(H,E,fh,fe)

  Cmd> stuff <- forstep(OUTS[grade(stuff$f,down:T)[1]],H,E,fh,fe);

The last step can be repeated to bring "in" variables.  Of course, in
practice, you want to examine the computed F-to-enter statistics to see
if another variable *should* be entered.

@@@@doing_back_step#Doing a backward step
You can do a backward step (variable deletion) using macro backstep().
One difference between backstep() and forstep() is that backstep()
determines the variable to eliminate, and then updates INS and OUTS; you
must tell forstep() which variable to include.  See backstep() for
details.  See also compf() which computes F-to-enter for variables not
in INS.
@@@@______

Both backstep() and compf() are OBSOLETE and are retained only for
backward compatibility.

@@@@see_also#Cross references
See also manova(), daentervar(), daremovevar(), dastepsetup(),
dastepstatus() and dasteplook().
@@@@______

====glscrit()#factor analysis
%%%%
result <- glscrit(logpsi, k, structure(r, nfact [, del])), REAL vector
  logpsi, integer scalar k, integer nfact > 0, small REAL scalar del
%%%%
@@@@usage#Usage
glscrit() is used by facanal() with method:"gls" to compute the
generalized least squares (GLS) criterion F_gls(x), x = log(psi) being
minimized, its gradient dF_gls(x)/dx and the current loadings and
eigenvalues.

trace((I - solve(r) %*% rhat) %*% (I - solve(r) %*% rhat))

result <- glscrit(logpsi, k, structure(r, nfact [,del]) computes a
scalar, vector or structure, depending on the value of integer k.

r is a REAL symmetric covariance or correlation matrix, nfact > 0 is the
number of factors, and del >= 0 is a small REAL scalar (default is 1e-5)
controlling how the gradient is computed.

When del = 0 (value used by facanal()), the gradient is computed
analytically.

When del > 0, the gradient vector has elements gradient[i] = (F_gls(x) +
del*(run(p)==i)) - F_gls(x))/del.

@@@@result#Result
The result computed is as follows:
   k           result
   0           scalar F_gls(x), x = logpsi
   1           vector gradient dF_gls(x)/dx
  -2           structure(loadings, eigenvals)

When r is not positive definite, MISSING (k = 0) or rep(MISSING,
nrows(s)) (k = 1) is returned.

No argument checking is done.

@@@@see_also#Cross references
See also facanal(), mlcrit(), ulscrit()

====glsfactor()#factor analysis
%%%%
glsfactor(s, m [, quiet:T,silent:T,start:psi0,uselogs:T], maxit:nmax,
  minit:nmin, print:T,crit:crvec]), REAL symmetric matrix s with no
  MISSING values, integers nmax > 0 and nmin > 0, crvec =
  vector(numsig, nsigsq, delta)
%%%%
@@@@usage#Usage
You use glsfactor() to find GLS (generalized least squares) estimates of
the vector psi of factor analysis uniquenesses and loading matrix L by
means of nonlinear least squares.  It uses macro glsresids() to compute
residuals used by non-linear least squares macro levmar().

glsfactor(r, m) estimates uniquenesses and loadings for a m-factor
analysis based on p by p symmetric correlation or variance-covariance
matrix r.  m > 0 must be an integer < (2*p + 1 - sqrt(8*p+1))/2.  The
default starting value for the uniquenesses is psi0 = 1/diag(solve(r)).

glsfactor() prints the estimated uniquenesses psihat and loading matrix,
the minimized criterion value, and vals, a vector of numbers computed
from certain eigenvalues; see below.

@@@@value_returned#Value returned
The result is an "invisible" (assignable but not automatically printed)
variable of the form
   structure(psihat:psi,loadings:L,criterion:crit,eigenvals:vals,\
      status:iconv,iter:n)
The values of these components are as follows

   psi     length p vector of estimated uniquenesses
   L       p by m matrix of estimated loadings satisfying
           L' %*% dmat(1/psi) %*% L is diagonal
   crit    sum of squared "residuals" of r from rhat
   vals    1/v - 1, where v is the vector of eigenvalues of
           dmat(psi) %*% solve(r) in increasing order
   iconv   integer >= 0 indicating convergence status, with 0 and 4
           indicating non-convergence
   n       number of iterations

@@@@keywords#Keyword phrase arguments
There are several keyword phrases that can be used as arguments to
control the behavior of glsfactor().
  start:psi0     psi0 a REAL vector of length p with min(psi0) > 0
                 to be used as starting values for the iteration
  uselogs:T      log(psi) is used in the iteration instead of psi;
                 this ensures psi does not become negative and may
                 speed convergence
  maxit:nmax     No more than nmax iterations are to be used; the
                 default is 30
  minit:nmin     At least nmin iterations will be performed; the
                 default is 1
  crit:crvec     crvec = vector(numsig, nsigsq, delta), 3 criteria
                 for convergence; a negative criterion is ignored;
                 see macro levmar() for details
  quiet:T        uniquenesses, loadings and vals are not printed;
                 only certain warning messages are printed
  silent:T       nothing is printed except error messages
  print:T        macro levmar() will print a status report on every
                 iteration

@@@@convergence_status_indicator#Convergence status indicator
Values of iconv indicate the following situations
 iconv  Meaning
   0    Not converged
   1    Converged with relative change in psi or log(psi) < 10^-numsig
   2    Converged with relative change in rss <= 10^-nsigsq
   3    Converged with ||gradient|| < delta
   4    Interation stopped; could not reduce rss.

@@@@caveats#Caveats
Without uselogs:T, there is no guarantee that min(psihat) > = 0.

Even with uselogs:T, there is no guarantee that the minimum is inside
the permissible range (r - dmat(psihat) positive semi-definite);
however, this often leads to an 'argument to solve() singular' error
message.

@@@@invisible_variable#Invisible variable _PSIGLS
Everytime it is called by levmar(), macro glsresids() saves a copy of
the current value of psi in invisible variable _PSIGLS.  Type
'print(_PSIGLS)' to see this value.
@@@@______

glsfactor() is very much a work in progress, that may in fact be
abandoned in favor of other methods using optimization macros in macro
file Math.mac distributed with MacAnova.

@@@@see_also#Cross references
See also ulsfactor() and glsresids().
@@@@______

====glsresids()#factor analysis
%%%%
glsresids() is used by macro glsfactor() to do GLS factor extraction by
a nonlinear least squares method.
%%%%
@@@@usage#Usage
glsresid(theta,0,vector(r),m) returns a length p^2 "residual vector"
computed from but not identical with dmat(p,1) - solve(r) %*% rhat,
where r is a REAL p by p correlation or variance-covariance matrix.  m >
0 is an integer specifying the number of factors to extract.

When __USELOGSGLS, a LOGICAL scalar defined by glsfactor(), has value
False, argument theta is interpreted as psi, the length(p) REAL vector
of uniquenesses.  When __USELOGSGLS is True, theta is interpreted as
log(psi).  glsresid saves a copy of psi in variable _PSIGLS.

rhat = dmat(psi) + V, where V is the best rank m approximation to r -
dmat(psi) in a generalized least squares sense.

glsresids() is used by macro glsfactor() which does generalized least
squares (GLS) factor extraction.
@@@@______

====goodfit()#factor analysis
%%%%
goodfit(s,lhat,psihat), REAL symmetric matrix lhat REAL vector or
  matrix, psihat REAL vector, all with no MISSING values
%%%%
@@@@usage#Usage
goodfit(r, L, psi), where r is a p by p symmetric correlation or
variance-covariance matrix, L is a REAL p by m matrix of purported
factor loadings, and psi is a length p REAL vector of factor analysis
uniquenesses computes three measures of lack of fit of r to the factor
analytic approximation rhat = L %*% L' + dmat(psi).

If dev = vector(r - rhat) is the length p^2 vector of deviations of r
from rhat, then goodfit() returns
       vector(max(abs(dev)), sum(dev^2), sum(abs(dev)))

When L and psi have been computed by some factor analysis estimation
method, and one or more of the elements of the result is large it may
indicate lack of fit of r to the factor analytic model.

These are intended to describe how bad the fit is, not to test it.
@@@@______

====groupcovar()#covariance,descriptive statistics
%%%%
groupcovar(groups, y), factor or vector of positive integers groups,
  REAL matrix y with no MISSING values
%%%%
@@@@usage#Usage
groupcovar(G, y), where G is a length N factor or vector of positive
integers and y is a REAL N by p data matrix with no MISSING values,
computes group means and the pooled sample covariance matrix for the
groups defined by the elements of G.  It is an error if no group has at
least 2 members.

@@@@value_returned#Value returned
The result is structure(n:n, means:ybar, covariance:V).
  n = vector(n1, ...,ng) is the vector of group sample sizes, where g =
      max(G).
  ybar is a g by p REAL matrix with ybar[i,] = sample mean for
      group i.  When group i is empty, ybar[i,] is all MISSING.
  V is the pooled variance-covariance matrix.  When all the groups are
      non-empty, V = ((n1-1)*S1 + (n2-1)*S2 + ... + (ng-1)*Sg)/(N - g)
      where S1, ..., Sg are the sample variance-covariance matrices for
      the g groups.

@@@@see_also#Cross reference
See also tabs().
@@@@______

====hotellval()#manova,hotelling tsq,hypothesis tests
%%%%
hotellval(x [, pval:T]), REAL matrix x with no MISSING elements
%%%%
@@@@usage#Usage
hotellval(x), where x is a n by p REAL matrix with no MISSING elements,
returns a REAL scalar containing Hotelling's T^2 statistic for testing
the null hypothesis H0:  mu = 0.  mu' is the expectation of each row of
x when the rows are assumed to be random sample from a multi-variate
population.

It is an error if n <= p.

hotellval(x, pval:T) does the same, except a P-value is also computed
under the assumption that the rows of x are independent multivariate
normal.  The result is structure(hotelling:tsq, pvalue:pval), where tsq
and pval are REAL scalars.

@@@@examples
Examples:
  Cmd> hotellval(x - mu0', pval:T)

where mu0 is a REAL vector of length p, computes a test statistic and
associated P-value for testing H0: mu = mu0.

  Cmd> hotellval(x1 - x2, pval:T)

where x1 and x2 are both n by p, computes Hotelling's paired P^2
statistic for testing H0: E(x1 - x2) = 0, together with a P-value.

@@@@see_also#Cross references
See also hotell2val(), tval() and t2val().
@@@@______

====hotell2val()#manova,hotelling tsq,hypothesis tests
%%%%
hotell2val(x1, x2 [, pval:T]), REAL matrices x1 and x2 with no MISSING
  elements and ncols(x1) = ncols(x2)
%%%%
@@@@usage#Usage
hotell2val(x1, x2), where x1 and x2 are REAL matrices with n1 and n2
rows and p columns, returns Hotelling's two-sample T^2 statistic for
testing the null hypothesis that the rows of x1 have the same mean as
the rows of x2 when they are considered as independent random samples
from two multivariate populations.  The statistic assumes that the
variance-covariance matrix is the same in the two populations.

It is an error if n1 + n2 <= p.

hotell2val(x1, x2, pval:T) does the same, except it also computes a
P-value under the assumption that the rows of x1 and of x2 constitute
independent random samples from multivariate normal distributions with
the same variance-covariance matrix.  The result is
structure(hotelling:tsq, pvalue:pval).

@@@@see_also#Cross references
See also hotellval(), tval() and t2val().
@@@@______

====mlcrit()#factor analysis
%%%%
result <- mlcrit(logpsi, k, structure(r, nfact [, del])), REAL vector
  logpsi, integer scalar k, integer nfact > 0, small REAL scalar del
%%%%
@@@@usage#Usage
mlcrit() is used by facanal() with method:"mle" to compute the maximum
likelihood ML criterion F_ml(x), x = log(psi) being minimized, its
gradient dF_ml(x)/dx and the current loadings and eigenvalues.

result <- mlcrit(logpsi, k, structure(r, nfact [,del]) computes a
scalar, vector or structure, depending on the value of integer k.

r is a REAL symmetric covariance or correlation matrix, nfact > 0 is the
number of factors, and del >= 0 is a small REAL scalar (default is 1e-5)
controlling how the gradient is computed.

When del = 0 (value used by facanal()), the gradient is computed
analytically.

When del > 0, the gradient vector has elements gradient[i] = (F_ml(x) +
del*(run(p)==i)) - F_ml(x))/del.

@@@@result#Result
The result computed is as follows:
   k           result
   0           scalar F_ml(x), x = logpsi
   1           vector gradient dF_ml(x)/dx
  -2           structure(loadings, eigenvals)

When r is not positive definite, MISSING (k = 0) or rep(MISSING,
nrows(s)) (k = 1) is returned.

No argument checking is done.

@@@@see_also#Cross references
See also facanal(), glscrit(), ulscrit()

====mulvar_index*#
%%%%
Macros related to discriminant analysis
  discrim(), jackknife(), discrimquad(), probsquad()
Macros related to subset discriminant analysis
  dastepsetup(), dastepstatus(), dasteplook(), daentervar(),
  daremovevar(), davarid(), dastepselect(), dascreen()
Obsolete macros related to subset discriminant analysis
  compf()*, forstep()*, backstep()*
Macros related to Hotelling's T^2
  hotellval(), hotell2val()
Macros related to MANOVA hypothesis tests
  cumpillai(), cumtrace(), cumwilks(), seqF()
Macros related to factor analysis
  facanal(), glscrit(), mlcrit(), ulscrit(), glsfactor(), glsresids(),
  ulsfactor(), ulsresids(), stepgls(), stepml(), stepuls(), goodfit()
Data generation macro
  rmvnorm()
Macros related to plotting of multivariate data
  chiqqplot()
Other data analysis macros
  covar()*, distcomp(), groupcovar(), standardize()
Other macros
  mulvarhelp()
* Obsolete; retained for backward compatibility
%%%%
@@@@discriminant_analysis_macros%Discriminant analysis macros
Entries related to discriminant analysis
 discrim()     Macro to compute a linear discriminant function, assuming
               multivariate normal distributions with same variance-
               covariance matrices
 discrimquad() Macro to compute coefficients for quadratic
               discrimination, assuming multivariate normal
               distributions with differing variance-covariance matrices
 jackknife()   Macro to estimate mis-classification probabilities using
               the leave-one-out method
 probsquad()   Macro to compute posterior probabilities based on
               quadratic discriminant functions computed by macro
               discrimquad()

@@@@variable_selection%Selecting variables in discriminant analysis
 dastepselect()Macro to do forward, backward or mixed stepwise variable
               selection
 dascreen()    Macro to select "best" subset of variables
 daentervar()  Macro to enter a variable in stepwise variable selection
 daremovevar() Macro to remove a variable in stepwise variable selection
 dasteplook()  Macro to access components of structure _DASTEPSTATE
 dastepsetup() Macro to initialize structure _DASTEPSTATE
 dastepstatus()Macro to summarize the contents of structure _DASTEPSTATE
 davarid()     Utility macro to decode variable number and name
 _DASTEPSTATE  Description of structure _DASTEPSTATE which summarizes
               the current state of the variable selection process

@@@@obsolete_discriminant_macros%Obsolete discriminant macros
The following are obsolete but are retained for backward compatibility
 backstep()    Macro to do single backward step in variable selection
 compf()       Macro to compute F-to-enter in stepwise variable
               selection
 forstep()     Macro to do single forward step in variable selection

@@@@hotellings_t_squared_macros%Computing Hotelling's T^2
Macros related to Hotelling's T^2
 hotellval()   Macro to compute single sample Hotelling's T^2
 hotell2val()  Macro to compute two sample Hotelling's T^2

@@@@manova_tests%Tests of MANOVA hypothesis
 cumpillai()   Macro to compute P-values for Pillai's trace statistic
 cumtrace()    Macro to compute P-values for Hotellings' and Pillai's
               trace statistics
 cumwilks()    Macro to compute P-values for likelihood ratio test
               statistic
 seqF()        Macro to compute sequential F-statistics

@@@@factor_analysis_macros%Macros related to factor analysis
 glsfactor()   Macro to do generalized least squares (GLS) factor
               extraction using a method based on a non-linear least
               squares computation
 glsresid()    Macro used by glsfactor()
 ulsfactor()   Macro to do unweighted least squares (ULS) factor
               extraction using a method based on a non-linear least
               squares computation
 ulsresids()   Macro used by ulsfactor()
 goodfit()     Macro that computes several descriptive measures of close
               ness of the sample correlation or covariance matrix to an
               estimate based on the factor analytic model
 stepgls()     Macro to compute one step of an iterative method of
               generalized least squares (GSL) factor extraction.
 stepml()      Macro to compute one step of an iterative method of
               maximum likelihood (ML) factor extraction.
 stepuls()     Macro to compute one step of an iterative method of
               unweighted least squares (ULS) factor extraction.  This
               is essentially the iterated principal factor method.

@@@@random_multivariate_normal%Data generation macro
 rmvnorm()     Macro to compute a pseudorandom sample from a
               multivariate normal distribution

@@@@other_data_analysis_macros%Other data analysis macros
 chiqqplot()   Macro to make a chi-squared or sqrt(chi-squared) Q-Q
               plot using generalized distances, usually computed
               from residuals from a multivariate linear model
 covar()       Macro to compute the sample mean and covariance matrix
               from a multivariate sample
 distcomp()    Macro to compute generalized distances from of the
               rows of a matrix from the vector of column means
 groupcovar()  Macro to compute the group means and pooled variance-
               covariance matrix from multi-group data
 standardize() Macro to recenter and scale columns of data matrix; by
               by default, result has means = 0 and standard deviations
               = 1

@@@@other_macros%Other macros
 mulvarhelp()  Macro to print help and usage information related to
               macros in file Mulvar.mac.

====jackknife()#discrimination,classification
%%%%
probs <- jackknife(groups,y [,prior:P]), factor or vector of positive
  integers groups, REAL matrix y and positive vector P with no MISSING
  elements
%%%%
@@@@usage#Usage
jackknife(G, y), where G is a factor or vector of positive integers of
length n and y is a REAL n by p matrix with no MISSING elements, carries
out a jackknife validation of linear discriminant functions designed to
discriminate among the g groups defined by the levels of G.

@@@@what_it_does#What jackknife() does
When you try to estimate the error rate of a classification method by
counting the errors it makes in classifying the cases in the "training
sample", the data set you are using to estimate the method, your
estimate is biased in an optimistic direction.  That is, the proportion
of cases misclassified in the training sample tends to be smaller than
the proportion of cases misclassified in new sample (validation sample).
jackknife() attempts to avoid this bias by classifying each case in the
training sample with linear discriminant functions computed from all the
other cases in the training sample.  This is the "leave-one-out" method,
sometimes called the Lachenbruch-Mickey method.

@@@@value_returned#Value returned
Macro jackknife() returns a n by g+1 matrix probs.

probs[i,j], for j = 1,...,g is an estimate of the posterior probability
that the data in y[i,] were derived from population j.

probs[i,g+1] is an integer from 1 to g indicating the population in
which the case should be classified, that is the population for which
the posterior probability is largest.

@@@@posterior_probabilities#Posterior probabilities
For each i, 1 <= i <= n, the posterior probabilities probs[i,j], j =
1,..., g are computed as follows.

The linear discriminant function based on y[-i,], that is using all the
data except row i, and the discriminant functions scores for the data
in y[i,] are computed.  From these the posterior probabilities are
computed assuming equal prior probability 1/g for each of the groups.
Each group is assumed to be multivariate normal with the same variance-
covariance matrix in each group.

Because the discriminant functions used to classify y[i,] are computed
without using y[i,], the method is close to unbiased.

@@@@prior_keyword#Keyword prior
jackknife(G,y,prior:P), where P is a REAL vector of length n with no
MISSING elements, does the same except the posterior probabilities are
computed using P.

@@@@example#Example
Here is how you might use jackknife() to estimate the expected
probability of misclassification, assuming the prior probability that a
randomly selected case comes from population j is P[j].

  Cmd> probs <- jackknife(G, y, prior:P)

  Cmd> n <- tabs(,G,count:T) # vector of sample sizes

  Cmd> misclassprob <- tabs(,G,probs[,g+1],count:T)/n

  Cmd> misclassprob[hconcat(run(g),run(g))] <- 0# set diags to 0

  Cmd> sum(misclassprob' %*% P)

The last line computes the estimated probability case randomly selected
from a group with prior probabilitys P will be misclassified by
linear discriminant functions estimated from y.  misclassprob[i,j] with
i != j is an estimate that a case from population i is misclassifed as
population j.
@@@@______

This version of jackknife() is relatively fast since it computes the
successive leave-one-out discriminant functions by modification of the
discriminant functions using all the data, rather than starting fresh.

====mulvarhelp()#
%%%%
mulvarhelp(topic1 [, topic2 ...] [,usage:T])
mulvarhelp(index:T)
%%%%
@@@@usage#Usage
mulvarhelp(topicname) prints help on a topic related to file
mulvar.mac.  Usually topicname is the name of a macro in the file.

When quoted, topicname may contain "wildcard" characters "*" and "?".
You can also use help keyword 'key'.  See help() for details.

mulvarhelp(topicname1, topicname2, ...) prints help on more than one
topic.

mulvarhelp(topicname1 [, topicname2 ...], usage:T) prints just a brief
summary of usage for the each topic.
@@@@______

====mvngen()#
%%%%
mvngen() has been renamed rmvnorm().  Type 'usage(rmvnorm)'.
%%%%
@@@@usage#Usage
mvngen() has been renamed rmvnorm().  Type 'help(rmvnorm)'.

====probsquad()#classification,discrimination
%%%%
probsquad(x,structure(Q,L,addcon,grandmean) [, prior:P]), REAL matrix x
  with no MISSING elements, argument 2 a structure as returned by macro
  discrimquad(), REAL vector P of prior probabilities with P[i] > 0,
  default = rep(1/g,g)
%%%%
@usage
Suppose info = structure(Q,L,addcon,grandmean) as computed by
discrimquad(groups, y) summarizes quadratic discriminant functions
estimated from a n by p matrix of data from g groups.

probsquad(x, info), where x a REAL vector of length p, returns the
length g vector of posterior probabilities, assuming multivariate
normality of each groups and equal prior probability 1/g for each group.

probsquad(x, info, prior:P), where P is a length g vector with P[j] >
0, does the same using P[j]/sum(P) as prior probability of group j.

In both these usages x can be a m by p matrix, each row of which is an
observation and probsquad() returns a m by g matrix of posterior
probabilities.

@@@@see_also#Cross reference
See probsquad() for information about the form of info.
@@@@______

====rmvnorm()#random numbers
%%%%
rmvnorm(n , p), integers n > 0, p > 1
rmvnorm(n, sigma), REAL positive definite symmetric matrix sigma,
  nrows(sigma) > 1
rmvnorm(n, p or sigma, mu), REAL row or column vector mu with length(mu)
  = p or ncols(sigma)
%%%%
@@@@usage#Usage
rmvnorm(n, p), where n > 0 and p > 0 are integer scalars, returns an n
by p matrix whose rows are a random sample from the standard
multivariate normal distribution MVN(0, I_p), where I_p is the p by p
identity matrix.

rmvnorm(n, p, mu) where mu is a REAL row or column vector of length p
does the same, except the rows of the result are MVN(mu, I_p).

rmvnorm(n, sigma [,mu]), where sigma is a positive definite p by p
symmetric matrix with p > 1, does the same, except the rows of the
result are MVN(0, sigma) or MVN(mu, sigma).

When p = 1, use mu + sd*rnorm(n) or mu + sd*rmvnorm(n,1) to generate a
normal sample with mean mu and standard deviation sd.

@@@@see_also#Cross reference
See also rnorm().

@@@@______

====seqF()#manova,hypothesis tests
%%%%
result <- seqF(hypTerm [,errTerm] [,model:Model] [,order:Jorder]),
  hypTerm and errTerm positive integer or CHARACTER scalars, Model a
  CHARACTER scalar and Jorder integer vector containin permutation of
  run(p), p = number of variables.
%%%%
@@@@usage#Usage
result <- seqF(hypTerm,errTerm,model:Model,order:Jorder) computes
sequential F-tests of a MANOVA hypothesis.  Model is a CHARACTER scalar
or quoted string specifying a MANOVA GLM.  hypTerm and errTerm are
integer term numbers or CHARACTER scalar term names identifying the term
being tested and the appropriate error term.  Jorder is a permutation of
1, 2, ..., p = number of columns in the response matrix.  It specifies
the order the columns are to be tested in.

The only required argument is hypTerm.  When errTerm is omitted, the
last term is assumed to be the error term.  When model:Model is omitted,
variable STRMODEL is used, if it is defined.  When order:Jorder is
omitted, Jorder = run(p).

result is structure(f:fstats, fh:hypDf, fe:errorDf) where fstats, fh,
and fe are length p vectors of F-statistics, numerator degrees of
freedom and denominator degrees of freedom, respectively.

fstats[1] is an ANOVA F-statistic for y[,Jorder[1]].  For i > 1,
fstats[i] is an analysis of covariance F-statistic for y[,Jorder[i]]
with the columns of y[,Jorder[run(i-1)]] as covariates.

@@@@example#Example
  Cmd> list(variety,iris)
  iris            REAL   150   4     (labels)
  variety         REAL   150   1     FACTOR with 3 levels (labels)

  Cmd> seqF("variety",model:"iris=variety")
  component: f
        SepLen      SepWid      PetLen      PetWid
        119.26       94.13      310.26      24.904
  component: fh
        SepLen      SepWid      PetLen      PetWid
             2           2           2           2
  component: fe
        SepLen      SepWid      PetLen      PetWid
           147         146         145         144

  Cmd> result <- seqF("variety",model:"iris=variety",order:vector(2,1,4,3))

  Cmd> result
  component: f
        SepWid      SepLen      PetWid      PetLen
         49.16      189.65      272.24       35.59
  component: fh
        SepWid      SepLen      PetWid      PetLen
             2           2           2           2
  component: fe
        SepWid      SepLen      PetWid      PetLen
           147         146         145         144

  Cmd> cumF(result$f,result$fh,result$fe,upper:T) # P-values
  (1)   4.492e-17  2.5579e-41  8.0581e-50  2.7562e-13

@@@@see_also#Cross references
See also manova(), cumF().

%%%%

====standardize()#transformations
%%%%
ynew <- standardize(y [,locs [,scales]]), n by p REAL matrix y, locs and
  scales REAL scalars or row or column vectors of length p, defaults
  mean and standard deviations
%%%%
@@@@usage#Usage
ynew <- standardize(y), where y is a n by p REAL matrix, creates a n by
p REAL matrix ynew with ynew[i,j] = (y[i,j] - ybar[j])/sd[j], where ybar
and sd are vectors of column means and standard deviations of y.  ynew
will have column means = 0 and column standard deviations = 1.

ynew <- standardize(y,locs), where locs is a REAL row or column vector
of length p, does the same except ynew[i,j] = (y[i,j] - locs[j])/sd[j].

ynew <- standardize(y,locs,scales), where scales is a REAL row or column
vector of length p, does the same except ynew[i,j] = (y[i,j] -
locs[j])/scales[j].

@@@@scalar_locs_or_scales#Scalar locs or scales
If locs is a scalar, it is expanded to rep(locs, p).

If scales is a scalar, it is expanded to rep(scales, p).

@@@@missing_values#Missing values
Any means or standard deviations needed are computed by describe(y).
MISSING elements are ignored except for a warning message.

If any elements of locs are MISSING, they are replaced by the mean of
the corresponding columns of y.

If any elements of scales are MISSING, they are replaced by the standard
deviations of the corresponding columns of y.

@@@@see_also#Cross reference
See also describe().

====stepgls()#factor analysis,iteration
%%%%
factorstats <- stepgls(s, psi, m [,nsteps:n] [, print:T]), REAL positive
  definite symmetric matrix s, REAL vector psi or structure with psi[1]
  a REAL vector, integers m > 0, n > 0
%%%%
@@@@usage#Usage
factorstats <- stepgls(r, psi, m) performs one step of an iteration
which attempts to extract m factors in factor analyis.  r is a p by p
correlation or variance-covariance matrix, psi is a REAL vector of trial
uniquenesses (specific variances) of length p with min(psi) > 0 and m >
0 is an integer.

The returned value factorstats = structure(psi:newpsi, loadings:L,
crit:criterion), where newpsi is an updated vector of uniquenesses, L a
p by m matrix of factor loadings satisfying L' %*% L is diagonal, and
criterion is a value of the criterion being minimized; see below.

factorstats <- stepgls(r, psi, m, nsteps:n), n > 0 an integer scalar,
does the same except n steps of the iteration are performed.

@@@@side_effect_variables#Side Effect Variables
In addition, stepgls() creates side effect variabls PSI, LOADINGS, and
CRITERION containing newphi, L and criterion.

Actually L is the loadings that minimize the criterion for argument
phi, not the output newphi, and criterion is the criterion assocated
with phi and L.

@@@@what_it_does#What stepgls() does
Each step reduces the generalized least squares (GLS) criterion =
trace((I - solve(r) %*% rhat) %*% (I - solve(r) %*% rhat)), where rhat
has the form rhat = dmat(psi) + V where V is the unique rank m matrix
that minimizes this criterion for given psi.

@@@@side_effect_variables#Side effect variables
stepgls() also creates side effect variabls PSI, LOADINGS, and CRITERION
containing newphi, L and criterion.

@@@@print_keyword#Keyword print
factorstats <- stemgsl(r, psi, m [,nsteps:n], print:T) does the same
except the updated psi and the criterion value are printed on each step.
To have updated results printed less frequently, say, every 50 steps,
use a for-loop which takes one step on each trip through the loop:

  Cmd> psi <- psistart # set starting value

  Cmd> for (i,1,500){psi <- stemgls(r,psi,m,print:i %% 50 == 0);;}

prints out psi and and the criterion every 50 steps.

@@@@caveats#Caveats
The iteration performed by stepgls() is similar but not identical to
what is known Principal Factor Iteration.  It is not unusual for it to
converge very slowly.  In addition, it is possible for a step to produce
a psi with min(psi) < 0 or r - dmat(psi) not non-negative definite.
When this happens stepgls() aborts.  You can retrieve the most recent
uniquenesses and loadings from variables PSI and LOADINGS.

It is usually preferable to use a method that more directly minimizes
the criterion, such as macro glsfactor().  In addition, one of the
function minimization macros, dfp, bfs or broyden could be used to
minimize the criterion as a function of psi or log(psi).
@@@@______

====stepml()#factor analysis,iteration
%%%%
factorstats <- stepml(s, psi, m [,nsteps:n] [, print:T]), REAL positive
  definite symmetric matrix s, REAL vector psi or structure with psi[1]
  a REAL vector, integers m > 0, n > 0
%%%%
@@@@usage#Usage
factorstats <- stepml(r, psi, m) performs one step of an iteration which
attempts to extract m factors in factor analyis.  r is a p by p
correlation or variance-covariance matrix, psi is a REAL vector of trial
uniquenesses (specific variances) of length p with min(psi) > 0 and m >
0 is an integer.

The returned value factorstats = structure(psi:newpsi, loadings:L,
crit:criterion), where newpsi is an updated vector of uniquenesses, L a
p by m matrix of factor loadings satisfying L' %*% L is diagonal, and
criterion is a value of the criterion being minimized; see below.

factorstats <- stepml(r, psi, m, nsteps:n), n > 0 an integer scalar,
does the same except n steps of the iteration are performed.

Actually L is the loadings that minimize the criterion for argument
phi, not the output newphi, and criterion is the criterion assocated
with phi and L.

@@@@what_it_does#What stepml() does
Each step reduces the likelihood (ML) criterion = log(det(rhat)) -
trace(r %*% solve(rhat)) - log(det(r)) - p, where rhat has the form
rhat = dmat(psi) + V where V is the unique rank m matrix that minimizes
this criterionfor given psi.

@@@@side_effect_variables#Side effect variables
stepml() also creates side effect variabls PSI, LOADINGS, and CRITERION
containing newphi, L and criterion.

@@@@print_keyword#Keyword print
factorstats <- stepml(r, psi, m [,nsteps:n], print:T) does the same
except the updated psi and the criterion value are printed on each step.
To have updated results printed less frequently, say, every 50 steps,
use a for-loop which takes one step on each trip through the loop:

  Cmd> psi <- psistart # set starting value

  Cmd> for (i,1,500){psi <- stepml(r,psi,m,print:i %% 50 == 0);;}

prints out psi and and the criterion every 50 steps.

@@@@caveats#Caveats
The iteration performed by stepml() is similar but not identical to what
is known Principal Factor Iteration.  It is not unusual for it to
converge very slowly.  In addition, it is possible a step to produce a
psi with min(psi) < 0 or r - dmat(psi) not non-negative definite.  When
this happens stepml() aborts.  You can retrieve the most recent
uniquenesses and loadings from variables PSI and LOADINGS.

It is usually preferable to use a method that more directly minimizes
the criterion such as facanal().

@@@@see_also#Cross references
See also facanal(), stepgls(), stepuls()
@@@@______

====stepuls()#factor analysis,iteration
%%%%
factorstats <- stepml(s, psi, m [,nsteps:n] [, print:T]), REAL positive
  definite symmetric matrix s, REAL vector psi or structure with psi[1]
  a REAL vector, integers m > 0, n > 0
%%%%
@@@@usage#Usage
factorstats <- stepuls(r, psi, m) performs one step of an iteration
which attempts to extract m factors in factor analyis.  r is a p by p
correlation or variance-covariance matrix, psi is a REAL vector of trial
uniquenesses (specific variances) of length p with min(psi) > 0 and m >
0 is an integer.

The returned value factorstats = structure(psi:newpsi, loadings:L,
crit:criterion), where newpsi is an updated vector of uniquenesses, L a
p by m matrix of factor loadings satisfying L' %*% L is diagonal, and
criterion is a value of the criterion being minimized; see below.

factorstats <- stepuls(r, psi, m, nsteps:n), n > 0 an integer scalar,
does the same except n steps of the iteration are performed.

Actually L contains the loadings that minimize the criterion for
argument phi, not the output newphi, and criterion is the criterion
assocated with phi and L.

@@@@what_it_does#What stepuls() does
Each step reduces the unweighted least squares (ULS) criterion =
trace((r - rhat) %*% (r - rhat)), where rhat has the form rhat =
dmat(psi) + V where V is the unique rank m matrix that minimizes this
criterion for given psi.

@@@@side_effect_variables#Side effect variables
stepuls() also creates side effect variabls PSI, LOADINGS, and CRITERION
containing newphi, L and criterion.

@@@@print_keyword#Keyword print
factorstats <- stepuls(r, psi, m [,nsteps:n], print:T) does the same
except the updated psi and the criterion value are printed on each step.
To have updated results printed less frequently, say, every 50 steps,
use a for-loop which takes one step on each trip through the loop:

  Cmd> psi <- psistart # set starting value

  Cmd> for (i,1,500){psi <- stepuls(r,psi,m,print:i %% 50 == 0);;}

prints out psi and and the criterion every 50 steps.

@@@@caveats#Caveats
The iteration performed by stepuls() is also known as Principal Factor
Iteration.  It is not unusual for it to converge very slowly.  In
addition, it is possible a step to produce a psi with min(psi) < 0 or r
- dmat(psi) not non-negative definite.  When this happens stepuls()
aborts.  You can retrieve the most recent uniquenesses and loadings
from variables PSI and LOADINGS.

It is usually preferable to use a method that more directly minimizes
the criterion, such as macro facanal().

@@@@see_also#Cross reference
See also facanal(), stepml(), stepgls()

@@@@______

====ulscrit()#factor analysis
%%%%
result <- ulscrit(logpsi, k, structure(r, nfact [, del])), REAL vector
  logpsi, integer scalar k, integer nfact > 0, small REAL scalar del
%%%%
@@@@usage#Usage
ulscrit() is used by facanal() with method:"uls" to compute the
unweighted least squares (ULS) criterion F_uls(x), x = log(psi) being
minimized, its gradient dF_uls(x)/dx and the current loadings and
eigenvalues.

result <- ulscrit(logpsi, k, structure(r, nfact [,del]) computes a
scalar, vector or structure, depending on the value of integer k.

r is a REAL symmetric covariance or correlation matrix, nfact > 0 is the
number of factors, and del >= 0 is a small REAL scalar (default is 1e-5)
controlling how the gradient is computed.

When del = 0 (value used by facanal()), the gradient is computed
analytically.

When del > 0, the gradient vector has elements gradient[i] = (F_uls(x) +
del*(run(p)==i)) - F_uls(x))/del.

@@@@value_returned#Value returned
The result computed is as follows:
   k           result
   0           scalar F_uls(x), x = logpsi
   1           vector gradient dF_uls(x)/dx
  -2           structure(loadings, eigenvals)

When r is not positive definite, MISSING (k = 0) or rep(MISSING,
nrows(s)) (k = 1) is returned.

No argument checking is done.

@@@@see_also#Cross references
See also facanal(), glscrit(), mlcrit().

====ulsfactor()#factor analysis
%%%%
ulsfactor(s, m [, quiet:T,silent:T,start:psi0,uselogs:T], maxit:nmax,
  minit:nmin, print:T,crit:crvec]), REAL symmetric matrix s with no
  MISSING values, integers nmax > 0 and nmin > 0, crvec =
  vector(numsig, nsigsq, delta)
%%%%
@@@@usage#Usage
You use ulsfactor() to find ULS (unweighted least squares) estimates of
the vector psi of factor analysis uniquenesses and the loading matrix L
by means of nonlinear least squares.  It uses macro ulsresids() to
compute residuals required by non-linear least squares macro levmar().

ulsfactor(r, m) estimates uniquenesses and loadings for a m-factor
analysis based on p by p symmetric correlation or variance-covariance
matrix r.  m > 0 must be an integer < (2*p + 1 - sqrt(8*p+1))/2.  The
default starting value for the uniquenesses is psi0 = 1/diag(solve(r)).

@@@@printed_output#Printed output
ulsfactor() prints the estimated uniquenesses psihat and loading matrix,
the minimized criterion value, and vals, a vector of numbers computed
from certain eigenvalues; see below.

@@@@value_returned#Value returned
The result is an "invisible" (assignable but not automatically printed)
variable of the form
   structure(psihat:psi,loadings:L,criterion:crit,eigenvals:vals,\
      status:iconv,iter:n)
The values of these components are as follows

   psi     length p vector of estimated uniquenesses
   L       p by m matrix of estimated loadings satisfying
           L' %*% L is diagonal
   crit    sum of squared "residuals" = vector(r - rhat)
   vals    v, where v is the vector of eigenvalues of r - dmat(psi)
   iconv   integer >= 0 indicating convergence status, with 0 and 4
           indicating non-convergence
   n       number of iterations

@@@@keywords#Keyword phrase arguments
There are several keyword phrases that can be used as arguments to
control the behavior of ulsfactor().
  start:psi0     psi0 a REAL vector of length p with min(psi0) > 0
                 to be used as starting values for the iteration
  uselogs:T      log(psi) is used in the iteration instead of psi;
                 this ensures psi does not become negative
  maxit:nmax     No more than nmax iterations are to be used; the
                 default is 30
  minit:nmin     At least nmin iterations will be performed; the
                 default is 1
  crit:crvec     crvec = vector(numsig, nsigsq, delta), 3 criteria
                 for convergence; a negative criterion is ignored;
                 see macro levmar() for details
  quiet:T        uniquenesses, loadings and vals are not printed;
                 only certain warning messages are printed
  silent:T       nothing is printed except error messages
  print:T        macro levmar() will print a status report on every
                 iteration

@@@@convergence_status_indicator#Convergence status indicator
Values of iconv indicate the following situations
 iconv  Meaning
   0    Not converged
   1    Converged with relative change in psi or log(psi) < 10^-numsig
   2    Converged with relative change in rss <= 10^-nsigsq
   3    Converged with ||gradient|| < delta
   4    Interation stopped; could not reduce rss.

@@@@caveat#Caveat
With or without uselogs:T, there is no guarantee that the solution
found is admissible in the sense that r - dmat(psihat) is positive
semi-definite.  A warning is printed when is not.  With uselogs:T this
situation is sometimes indicated by an 'argument to solve() singular'
error message.

Without uselogs:T, there is no guarantee that min(psihat) > = 0.  A
warning is printed if it is not.

@@@@invisible_variable#Invisible variable _PSIGLS
Every time macro ulsresids() is called by levmar(), it saves a copy of
the current value of psi in invisible variable _PSIGLS.  Type
print(_PSIGLS) to see this value.
@@@@______

ulsfactor() is very much a work in progress that may in fact be
abandoned in favor of other methods using optimization macros in macro
file Math.mac distributed with MacAnova.

====ulsresids()#factor analysis
%%%%
ulsresids() is used by macro ulsfactor() to do ULS factor extraction by
a nonlinear least squares method.
%%%%
@@@@usage#Usage
ulsresid(theta,0,vector(r),m) returns the residual matrix r - rhat as a
vector of length p^2, where r is a REAL p by p correlation or variance-
covariance matrix.  m > 0 is an integer specifying the number of
factors to extract.

rhat = dmat(psi) + V, where V is the best rank m approximation to r -
dmat(psi) in the least squares sense.

ulsresids() is used by macro ulsfactor() to do unweighted least squares
(ULS) factor extraction.

@@@@use_of_logs#Use of logs
When __USELOGSULS, a LOGICAL scalar defined by ulsfactor(), has value
False, argument theta is interpreted as psi, the length(p) REAL vector
of uniquenesses.  When __USELOGSULS is True, theta is interpreted as
log(psi).  ulsresids() saves a copy of psi in variable _PSIULS.
@@@@______

_E_O_F_#This should be the last line, an internal End Of File marker