info   MACRO
) This file contains macros related to categorical data analysis
) Currently it contains the following:
)
) binomlimits(x,n,level [,type])       Compute confidence limits for
)                                      pi = P(success)
) chisq2way(twoWayTable)               Chi-squared test of independence
) condistrib(table,J)                  Find distributions conditional on
)                                      specified marginals
) contable()                           Chi-squared test of independence
)   contable(twoWayTable [,print:T])     Based on two-way table of counts
)   contable(a, b, [,print:T])           Based on factors a and b
)   contable(a, b, counts [,print:T])    Based on factors a and b and vector
)                                        counts
) catbargraph(y [,catnames:catNames] [barwidth:w] [,valuelab:T] \
)     [,graphics keywords])            Bar graph of y vs a categorical variable
) maketable(x [, ncols] [,rowlabs:labs1] [,collabs:labs2])
)                                      Create matrix with row and column
)                                      labels from data in x
) readtable(filename, ncols [,rowlabs:labs1] [,collabs:labs2] \
)     [, vecread keywords])            Read two way table into matrix,
)                                      adding row and column labels
) File version 010215  
%info%

maketable   MACRO DOLLARS
) Macro to create two-way table from a vector or matrix as
) as a matrix with row and column labels.
)
) Usage:
)  table <- maketable(x [, ncols] [,rowlabs:labs1] [,collabs:labs2])
)  x               REAL vector or matrix containing the data to go in the
)                  table.  If x is a matrix, the data is taken a
)                  column at a time
)  ncols           positive integer which must divide n = length(x).
)                  When x is a matrix, you can omit ncols in which
)                  case the default is ncols = ncols(x).
)  labs1, labs2    CHARACTER scalar or vector; if a vector, length must
)                  match n/ncols or ncols, n = number of data items.  The 
)                  defaults for labs1 and labs2 are "Row " and "Col ",
)                  respectively.  When scalars, they are expanded by
)                  appending the row or column number
))
)) Written by C. Bingham, kb@stat.umn.edu, 030416
)) Version 030416
# $S(x [, ncols] [,rowlabs:labs1] [,collabs:labs2])
@x <- argvalue($1,"data","real matrix")
if (isvector(@x) && $v < 2) {
	error("ncols is required when argument 1 is a vector")
}
@n <- length(@x)
@ncols <- if ($v == 1) {
	ncols(@x)
} else {
	argvalue($2,"ncols","positive count")
}
if (@n %% @ncols != 0) {
	error("number of columns does not divide number of data values")
}
@nrows <- @n/@ncols

@rowlabs <- keyvalue($K,"rowlab*","character vector",default:"Row ")
@collabs <- keyvalue($K,"collab*","character vector", default:"Col ")

if (!isscalar(@rowlabs) && length(@rowlabs) != @nrows){
	error("Number of row labels doesn't match number of rows")
}
if (!isscalar(@collabs) && length(@collabs) != @ncols){
	error("Number of column labels doesn't match number of columns")
}
matrix(delete(@x, return:T),@nrows,labels:structure(@rowlabs,@collabs))
%maketable%

readtable   MACRO DOLLARS
) Macro to read two-way table into a matrix with row and column labels.
) The file should be in format readable by vecread()
) Usage:
)  table <- readtable(fileName, ncols [,rowlabs:labs1] [,collabs:labs2]\
)    [,vecread keywords])
)  fileName        CHARACTER scalar containing the file name or
)                  string:String where String is a CHARACTER scalar or
)                  vector to be "read"
)  ncols           positive integer; must divide number of data items
)  labs1, labs2    CHARACTER scalar or vector; if a vector, length must
)                  match n/ncols or ncols, n = number of data items.  The 
)                  defaults for labs1 and labs2 are "Row " and "Col ",
)                  respectively.  When scalars, they are expanded by
)                  appending the row or column number
) vecread() is used to read the file.  The default value of 'quiet' is F
) except with string:String
)
)) At some point, this should be upgraded so that it gets the number
)) of columns from the file and column labels from the first line in the file.
)) Written by C. Bingham, kb@stat.umn.edu, 010214
)) Version 010215
# $S(fileName, ncolumns [,rowlabs:labs1] [,collabs:labs2]\
#     [,vecread keywords])
@string <- keyvalue($K,"string","string")
if (isnull(@string) && $v != 2 || !isnull(@string) && $v != 1){
	error("usage: $S(fileName, ncolumns [,vecread keywords])",macroname:F)
}
@ncols <- argvalue($2, "number of columns","positive count")
@rowlabs <- keyvalue($K,"rowlab*","character vector",default:"Row ")
@collabs <- keyvalue($K,"collab*","character vector", default:"Col ")
if (isnull(@string)){
	@name <- argvalue($1,"file name","string")
	@y <- if ($k != 0) {
		vecread(@name, $K, quiet:F, badkeyok:T)
	} else {
		vecread(@name, quiet:F)
	}
	delete(@name)
} else {
	@y <- vecread($K, badkeyok:T)
}
delete(@string)
@n <- length(@y)
if (length(@y) %% @ncols != 0){
	error("number of columns does not divide number of data values")
}
@nrows <- @n/@ncols
if (!isscalar(@rowlabs) && length(@rowlabs) != @nrows){
	error("Number of row labels doesn't match number of rows")
}
if (!isscalar(@collabs) && length(@collabs) != @ncols){
	error("Number of column labels doesn't match number of columns")
}
matrix(delete(@y, return:T),@ncols,labels:structure(@collabs,@rowlabs))'
%readtable%

chisq2way   MACRO DOLLARS
) Simple macro to do test of independence of r x c contingency table
) Usage:
)       result <- chisq2way(n) where n is a r by c table (matrix) of counts
) result is set to structure(n,expect,chisq,df)
) Access each part by result$n, result$expect, result$chisq and
) result$df
)
))Version of 990923
# $S(n) where n is a matrix of non-negative integers
if ($v != 1){
	error("usage: $S(n), n a matrix of non-negative integers")
}
@n <- argvalue($1,"table of counts","nonnegative integer matrix")
if(!haslabels(@n)){
	@n <- matrix(@n,labels:structure("Row ","Col "))
}
@colsums <- sum(@n)
@rowsums <- sum(@n')'
@N <- sum(vector(@n))
@e <- matrix(@rowsums * @colsums/@N,labels:getlabels(@n))
@chisq <- sum(vector((@n - @e)^2/@e))
@df <- (nrows(@n) - 1)*(ncols(@n) - 1)
structure(n:delete(@n,return:T),expect:delete(@e,return:T),\
	chisq:delete(@chisq,return:T), df:delete(@df,return:T))
%chisq2way%

contable      MACRO DOLLARS
) Macro to do test of independence of r x c contingency table
) Usage: contable(n [,print:T]) where n is a r by c matrix of counts >= 0
)        contable(a,b,n [,print:T]), where a and b are factors
)          of the same length and n is a vector of counts >= 0
)          of the same length as a and b
)        contable(a,b [,print:T]) means contable(a,b,rep(1,length(a)))
)    With print:T, A table of relative frequencies, row probabilities column
)    probabilities, expectations, totals, chi-squared and P-vallue is printed.
)    The value returned is structure(counts expect,probs,chisq,df,p_value)
)    This will not be printed but can be assigned
)) Written by C. Bingham, kb@stat.umn.edu
)) 010214 Version
# $S(n [,print:T]), n non-negative integer matrix
# $S(rowNumber,rowNumber [,counts] [,print:T]), rowNumber,rowNumber factors,
#    counts non-negative integer vector
if ($v > 3){
	error("usage: $S(n) or $S(rowNumber,colNumber [,counts])", macroname:F)
}
if ($v == 1){
	@n <- matrix(argvalue($1,"table of counts","nonnegative integer matrix"))
}else{
	@a <- factor(argvalue($1,"row factor","positive integer vector"))
	@b <- factor(argvalue($2,"column factor", "positive integer vector"))
	if (length(@a) != length(@b)){
		error("row and column factors have different lengths")
	}
	@counts <- if ($v == 2){
		rep(1,length(@a))
	}else{
		argvalue($3,"counts","nonnegative integer vector")
	}
	if (length(@a) != length(@counts)){
		error("counts and factors have different lengths")
	}
	@n <- tabs(delete(@counts,return:T),@a,@b,sum:T)
}
@print <- keyvalue($K,"print","TF",default:F)
if (!haslabels(@n)){
	@n <- matrix(@n,labels:structure("Row ","Col "))
}
@colsums <- sum(@n)
@rowsums <- sum(@n')'
@N <- sum(@rowsums)
if (@N == 0){
	error("all counts are 0")
}
@I <- run(nrows(@n))[@rowsums != 0]
@J <- run(ncols(@n))[@colsums' != 0]
if (length(@I) < nrows(@n)){
	print("WARNING: one or more all zero rows omitted")
	@n <- @n[@I,]
	@rowsums <- @rowsums[@I]
}
if (length(@J) < ncols(@n)){
	print("WARNING: one or more all zero columns omitted")
	@n <- @n[,@J]
	@colsums <- @colsums[@J]
}
@nrows <- nrows(@n)
@ncols <- ncols(@n)
if (@ncols == 1 || @nrows == 1){
	error("only one row or column with non-zero counts")
}
@df <- (@nrows-1)*(@ncols-1)
@e <- @colsums * @rowsums / @N
@chisq <- sum(vector((@n - @e)^2/@e))

@rowlabs <- getlabels(@n,1)
if (@print){
	print("         Contingency Table Analysis")
	print(paste(charwidth:9,"", charwidth:8,getlabels(@n,2),\
		charwidth:1,"|",charwidth:8,"Total",justify:"R"))
	for (@i, run(@nrows)){
		print(paste(charwidth:9,@rowlabs[@i],intwidth:8,@n[@i,],\
			charwidth:1,"|",@rowsums[@i]))
		print(paste(charwidth:9,justify:"R","Expect  ",format:"8.5g",@e[@i,],\
			charwidth:1,"|"))
		print(paste(charwidth:9,justify:"R","Rel Freq",format:"8.5f",@n[@i,]/@N,\
			charwidth:1,"|",@rowsums[@i]/@N))
		print(paste(charwidth:9,justify:"R","Row Prob",format:"8.5f",\
			@n[@i,]/@rowsums[@i],charwidth:1,"|"))
		print(paste(charwidth:9,justify:"R","Col Prob",format:"8.5f",\
			@n[@i,]/@colsums,	charwidth:1,"|"))
		print(paste(charwidth:9,"",charwidth:8,rep("--------",@ncols),\
			charwidth:1,"+","--------"))
	}
	print(paste(charwidth:9,"Total",intwidth:8,@colsums,\
		charwidth:1,"|",@N))
	print(paste(charwidth:9,"Rel Freq",format:"8.5f",@colsums/@N,\
		charwidth:1,"|",@N/@N))
	print(paste(sep:"","Chi-squared(",@df,") = ",@chisq,\
		", P-value = ",1 - cumchi(@chisq,@df)))
	delete(@i)
} #if (@print)
@result <- structure(counts:@n,expect:@e,probs:@n/@N,chisq:@chisq,\
df:@df,p_value:1 - cumchi(@chisq,@df))
delete(@n,@rowsums,@colsums,@e,@N,@chisq,@rowlabs,@J,@nrows,@ncols)
delete(@result,return:T,invisible:delete(@print,return:T))
%contable%

catbargraph   MACRO DOLLARS
) Macro to make a bar plot of a vector y against a categorical variable
) Usage:
)   catbargraph(y [,catnames:catNames] [barwidth:w]\
)               [,valuelab:T] [,graphics keywords])
)     y            non-MISSING REAL vector, usually counts or percentages
)     catNames     CHARACTER scalar or vector of length length(y)
)     w            REAL scalar 0 < w <= 1, width of bars, default = .75
)     valuelab:T   y[i] will be printed just above bar i
) 
) Let N = length(y).  Then N bars of width w and height y[i] will be drawn
) centered at i, i = 1, ..., N.
) 
) Without catnames:catNames, bars will be labelled with row labels from y,
) if any, or "Cat_1", "Cat_2", ... otherwise.
) 
) With catnames:catNames, bar i will be labeled with catNames[i] unless
) catNames is a scalar when it will be labelled with paste(catNames, i,
) sep:"").
)) 010214 written by C. Bingham (kb@stat.umn.edu)
)) Version 010214 
# $S(y [,catnames:charVec] [barwidth:w] [,valuelab:T] [,graphics keywords])
@Y <- argvalue($1,"y","nonmissing real")
if (!isvector(@Y)){
	@Y <- @Y'
	if (!isvector(@Y)){
		error("y is not row vector or column vector")
	}
}
@ncat <- nrows(@Y)
@keys <- if ($k > 0){
	structure($K)
} else {
	structure(notakey:NULL)
}

@catnames <- keyvalue(@keys,"catname*","character vector")
if (isnull(@catnames)){
	@catnames <- if (haslabels(@Y)) {
		getlabels(@Y,1)
	} else {
		getlabels(vector(run(@ncat),labels:"Cat_"))
	}
} elseif (isscalar(@catnames)){
	@catnames <- getlabels(vector(run(@ncat),labels:@catnames))
} elseif (length(@catnames) != @ncat){
	error("number of category names different from number of categories")
}
@valuelab <- keyvalue($K,"valuelab*","TF",default:F)
@width <- keyvalue(@keys,"barwidth","positive number", default:.75)
if (@width > 1){
	error("value for 'barwidth' must be <= 1")
}
@x <- run(@ncat) + vector(-1,-1,1,1)'*@width/2
@x <- vector(hconcat(@x,rep(?,@ncat))')[-5*@ncat]
@y <- vector(0,@Y' * vector(1,1,0,0,0))[run(5*@ncat-1)]

if (!@valuelab){
	lineplot(@x,@y,xticks:run(@ncat),xticklabs:@catnames,xticklen:0,\
		  $K,yaxis:F,xlab:" ",ylab:"$1",xmin:0,xmax:@ncat+1,ymax:1.1*max(@y),\
		  title:"Bar graph of $1")
} else {
	lineplot(@x,@y,xticks:run(@ncat),xticklabs:@catnames,xticklen:0,\
		  yaxis:F,xlab:" ",ylab:"$1",xmin:0,xmax:@ncat+1,ymax:1.1*max(@y),\
		  title:"Bar graph of $1",show:F)
	for (@i,1,@ncat){
		@catnames[@i] <- paste(@Y[@i],format:".3g")
	}
	delete(@i)
	addstrings(run(@ncat),@Y + .025*max(@Y),strings:@catnames,$K)
}
delete(@x,@y,@Y,@width,@catnames,@ncat)
%catbargraph%

binomlimits  macro  dollars
) Macro to to compute binomial confidence limits
) Usage:
) limits <- binomlimits(x,n,level [,type])
)   x        nonnegative integer scalar or vector
)   n        positive integer scalar or vector
)   level    REAL scalar between 0 and 1
)   type     one of "usual" or "wald", "score", "agresti", "wilson",
)            "exact","onestep"
) 
) If neither x or n are scalars, they must have the same length
) If both are scalars, limits = vector(lowerLimit,upperLimit)
) Otherise, limits = hconcat(lowerLimits,upperLimits)
) C. Bingham kb@stat.umn.edu
)) Version 000628
)) 020806 use upper:T on invF()
)) 030411 Added "wilson" as an equivalent to "agresti".  As a
))        result you now need to specify at least "wa" for Wald
))        and "wi" for Wilson
)) 030412 uses upper:T on invF() only for version >= 4.13
# usage: limits <- $S(x,n,level [,type])
@x <- argvalue($1,"x","nonnegative integer vector")
@n <- argvalue($2,"n","positive integer vector")
@coverage <- argvalue($3,"confidence level","positive scalar")

@type <- if($v >= 4) {
	argvalue($4,"confidence limit type","string")
} else {
	"usual" # same as "wald"
}						  

if (@coverage >= 1) {
	error(paste("confidence level",@coverage,">= 1"))
}
if (@coverage <= .5) {
	print(paste("WARNING: confidence level",@coverage,"<= .5"))
}
if (!isscalar(@x) && !isscalar(@n) && length(@x) != length(@n)) {
	error("x and n are vectors but don't have the same length")
}
if (sum(@x > @n) > 0) {
	error("Some x > n")
}
@N <- max(length(@x),length(@n))
@alpha <- 1 - delete(@coverage,return:T)
@z <- invnor(1 - @alpha/2)

@phat <- @x/@n
@lower <- 0*@phat
@upper <- @lower + 1
@x <-+ @lower
@n <-+ @lower
@alpha <-+ @lower

if (match("e*",@type,0,exact:F) != 0) {#exact
	@J0 <- @x != 0
	@Jn <- @x != @n
	@J <- @J0 && @Jn
	
	if (sum(@J) > 0) {
		@alpha[@J] <- @alpha[@J]/2
	}
# This uses 1 sided intervals for 0 and n;
# replace preceding by @alpha <- @alpha/2 if that's not what you want

	if (sum(@Jn) > 0) { # x != n
		@x1 <- @x[@Jn]
		@n1 <- @n[@Jn]
		@fval <- if (vecread(string:VERSION,silent:T)[1] >= 4.13) {
			invF(@alpha[@Jn],2*@x1+2,2*(@n1-@x1),upper:T)
		} else {
			invF(1 - @alpha[@Jn],2*@x1+2,2*(@n1-@x1))
		}
		@u <- (@x1+1)*@fval/(@n1-@x1)
		@upper[@Jn] <- @u/(1 + @u)
	}
	if (sum(@J0) > 0) { # x != 0
		@x1 <- @x[@J0]
		@n1 <- @n[@J0]
		@u <- @x1*invF(@alpha[@J0],2*@x1,2*(@n1-@x1+1))/(@n1-@x1+1)
		@lower[@J0] <- @u/(1 + @u)
	}
	delete(@J,@J0,@Jn,@x1,@n1,@u)
} elseif (match("u*",@type,0,exact:F) != 0 ||\
		  match("wa*",@type,0,exact:F) != 0) {#usual/Wald
	@d <- @z*sqrt(@phat*(1-@phat)/@n)
	@lower <- @phat - @d
	@upper <- @phat + @d
	delete(@d)
} elseif (match("a*",@type,0,exact:F) != 0 ||\
			match("wil*",@type,0,exact:F) != 0) {#Agresti/Wilson
	@phat <- (@x+2)/(@n+4)
	@d <- @z*sqrt(@phat*(1-@phat)/(@n + 4))
	@lower <- @phat - @d
	@upper <- @phat + @d
	delete(@d)
} elseif (match("s*",@type,0,exact:F) != 0) {#Score
	@d <- @z^2/@n
	@s <- @z*sqrt(@phat*(1-@phat)/@n + .25*@d/@n)/(1 + @d)
	@lower <- (@phat + @d/2)/(1 + @d) - @s
	@upper <- @lower + 2*@s
	delete(@d,@s)
} elseif (match("o*",@type,0,exact:F) != 0) {#Approximate
	@d <- @z*sqrt(@phat*(1-@phat)/@n)
	@lower <- @phat - @d
	@upper <- @phat + @d
	@J <- @lower > 0
	if (sum(@J) > 0) {
		@lower[@J] <- @phat[@J] - @z*sqrt(@lower[@J]*(1-@lower[@J])/@n[@J])
	}
	
	@J <- @upper < 1
	if (sum(@J) > 0) {
		@upper[@J] <- @phat[@J] + @z*sqrt(@upper[@J]*(1-@upper[@J])/@n[@J])
	}
		
	delete(@d,@J)
} else {
	error(paste("\"",@type,"\" is not a recognized confidence limit type",\
				sep:""))
}
@lower[@lower < 0] <- 0
@upper[@upper > 1] <- 1
@limits <- if(@N == 1) {
	vector(@lower,@upper)
} else {
	hconcat(@lower,@upper)
}
delete(@alpha,@z,@x,@n,@lower,@upper)
delete(@limits,return:T)
%binomlimits%

binomcover  macro  dollars
) Macro to to compute coverage of binomial confidence limits
) Usage:
) limits <- binomcover(p,n,level [,type] [,limits:T])
)   p        REAL scalar or vector between 0 and 1
)   n        positive integer scalar
)   level    REAL scalar between 0 and 1
)   type     one of "usual" or "wald", "score", "agresti" or "wilson", "exact"
)
) result is structure(p:p,coverage:c)
)  where p is the argument and c is vector the same length as p
)  with c[i] = P(lower_limit <= p[i] <= upper_limit)
)
) With limits:T, the result is structure(p:p,coverage:c,lower:l, upper:u)
)  where l and u are REAL vectors of length n+1 with l[i] and u[i] the lower
)  and upper limits when x = i+1.
)
) You can plot coverage vs p for 95% exact limits, say, simply by
)   lineplot(binomcover(p,n,.95,"exact")
) C. Bingham kb@stat.umn.edu
)) Version 000629
)) Version 030411 uses upper:T on invF() only for version 4.13 or later
))   accepts wilson as synonym of agresti
# usage: coverage <- $S(p,n,level [,type] [,limits:T])
@p <- argvalue($1,"p","positive vector")
@n <- argvalue($2,"n","positive integer scalar")
@coverage <- argvalue($3,"confidence level","positive scalar")

@type <- if($v >= 4) {
	argvalue($4,"confidence limit type","string")
} else {
	"u"
}						  
@limits <- keyvalue($K,"limit*","TF",default:F)
if (max(@p) >= 1) {
	error("p >= 1")
}

if (@coverage >= 1) {
	error(paste("confidence level",@coverage,">= 1"))
}
if (@coverage <= .5) {
	print(paste("WARNING: confidence level",@coverage,"<= .5"),\
		macroname:F)
}

@x <- run(0,@n)

@alpha <- 1 - delete(@coverage,return:T)
@z <- invnor(1 - @alpha/2)

@phat <- @x/@n
@lower <- 0*@phat
@upper <- @lower + 1

if (match("e*",@type,0,exact:F) != 0) {#exact
# This uses 1 sided intervals for 0 and n;
# replace preceding by @alpha <- @alpha/2 if that's not what you want
	@lower[1] <- 0
	@upper[1] <- 1 - @alpha^(1/@n)
	@lower[@n+1] <- @alpha^(1/@n)
	@upper[@n+1] <- 1
	if (@n > 2) {
		@x1 <- run(1,@n-1)
		@fval <- if (vecread(string:VERSION,silent:T)[1] >= 4.13) {
			invF(@alpha/2,2*@x1+2,2*(@n-@x1),upper:T)
		} else {
			invF(1 - @alpha/2,2*@x1+2,2*(@n-@x1))
		}
		@u <- (@x1+1)*@fval/(@n-@x1)
		@upper[@x1+1] <- @u/(1 + @u)
		@u <- @x1*invF(@alpha/2,2*@x1,2*(@n-@x1+1))/(@n-@x1+1)
		@lower[@x1+1] <- @u/(1 + @u)
		delete(@u,@x1)
	}
} elseif (match("u*",@type,0,exact:F) != 0 ||\
		  match("wa*",@type,0,exact:F) != 0) {#usual or Wald
	@d <- @z*sqrt(@phat*(1-@phat)/@n)
	@lower <- @phat - @d
	@upper <- @phat + @d
	delete(@d)
} elseif (match("a*",@type,0,exact:F) != 0 ||\
		  match("wi*",@type,0,exact:F) != 0) {#Agresti or Wilson
	@phat <- (@x+2)/(@n+4)
	@d <- @z*sqrt(@phat*(1-@phat)/(@n + 4))
	@lower <- @phat - @d
	@upper <- @phat + @d
	delete(@d)
} elseif (match("s*",@type,0,exact:F) != 0) {#Score
	@d <- @z^2/@n
	@s <- @z*sqrt(@phat*(1-@phat)/@n + .25*@d/@n)/(1 + @d)
	@lower <- (@phat + @d/2)/(1 + @d) - @s
	@upper <- @lower + 2*@s
	delete(@d,@s)
} elseif (match("o*",@type,0,exact:F) != 0) {#Approximate
	@d <- @z*sqrt(@phat*(1-@phat)/@n)
	@lower <- @phat - @d
	@upper <- @phat + @d
	@J <- @lower > 0
	if (sum(@J) > 0) {
		@lower[@J] <- @phat[@J] - @z*sqrt(@lower[@J]*(1-@lower[@J])/@n)
	}
	@J <- @upper < 1
	if (sum(@J) > 0) {
		@upper[@J] <- @phat[@J] + @z*sqrt(@upper[@J]*(1-@upper[@J])/@n)
	}
	delete(@d,@J)
} else {
	error(paste("\"",@type,"\" is not a recognized confidence limit type",\
				sep:""))
}
@result <- @p
for(@i,1,length(@p)) {
	@J <- run(0,@n)[@lower <= @p[@i] && @upper >= @p[@i]]
	@result[@i] <- cumbin(max(@J),@n,@p[@i]) - cumbin(min(@J)-1,@n,@p[@i])
}
delete(@alpha,@z,@x,@J)
if (delete(@limits,return:T)){
	structure(p:delete(@p,return:T),coverage:delete(@result,return:T),\
			lower:delete(@lower,return:T),upper:delete(@upper,return:T))
} else {
	delete(@lower,@upper)
	structure(p:delete(@p,return:T),coverage:delete(@result,return:T))
}
%binomcover%


condistrib  MACRO DOLLARS
) Macro to compute conditional distributions based on an array of
) counts.
) Usage:
) newtable <- condistrib(table, J)
)   table         Non-negative array with ndim(table) > 1, usually of counts
)   J             Vector of unique positive integers <= ndims(table),
)                 length(J) < ndims(table), specifying dimension on
)                 which to condition
)   newtable      Non-negative array the same size and shape as table
)                 such that the sum over all dimensions not in J
)				 is 1 for each combination of levels of the dimensions
)                 in J
) Currently (010214) only ndims(table) = 2 is supported for which J
) must be either 1 or 2 and the row or column conditional
) distributions are computed
)
)) Written by C. Bingham (kb@stat.umn.edu) 010214
)) Version 010214
# $S(table,J), table an array of counts, J integer scalar or vector
@table <- argvalue($1,"table","nonnegative")
@J <- argvalue($2,"conditioning variables","positive integer vector")
@ndims <- ndims(@table)
if (@ndims == 1){
	error("Argument 1 must have at least 2 dimensions")
}
if (length(@J) != length(unique(@J)) || max(@J) > @ndims){
	error("duplicate or too large variable numbers")
}
if (length(@J) == @ndims){
	error("you are trying to condition on all variables")
}
if (@ndims > 2){
	error("$S can handle only two-way tables at present")
}
@K <- vector(@J,run(@ndims)[-@J])
@table <- t(@table,@K)
@table <-/ sum(@table')'
@table <- t(@table,grade(@K))
@table
%condistrib%