info     0
) Macros for computing probability distributions
)  density(name, x, parameters [,keyword phrases])
)    name       one of "normal", "chisquare", "student", "F",
)               "beta", or "gamma"
)    x          REAL variable of arguments
)    parameters REAL scalar or vector
) The following are all interfaces that call density
)  densnor(x [, mu [,sigma]] [,keyword phrases])
)  densstu(x, df [,keyword phrases])
)  denschi(x, df [,keyword phrases])
)  densF(x, df1,df2 [,keyword phrases])
)  densbeta(x, nu1, nu2 [,keyword phrases])
)  densgamma(x, alpha [, mu] [,keyword phrases])
)
) Keyword phrases all:T, x:T or F, parameters:T or F, name:T or F
) Without keyword phrases, macros return density in the same
) size and shape as x.
) Otherwise they return structure([x ,] density [,parameters] [,name])
)) Written by C. Bingham (kb@stat.umn.edu)
))991030 Changed all to use DOLLARS, getmacros(density), added keywords
))010113 Corrected bug in Student's t

density   MACRO DOLLARS
) Usage:
)    density(name, x, parameters [,keyword phrases])
) or density(name, x [,keyword phrases])
)   name must be "normal", "chisquare", "student", "F", "beta", "gamma".
)   x is a REAL variable
)   parameters is a REAL scalar or vector
)   Keyword phrases all:T, x:T or F, parameters:T or F, name:T or F
) Without keyword phrases, returns density in the same size and
) shape as x
) Otherwise returns structure([x ,] density [,parameters] [,name])
)
) For normal, parameters = vector(mu [, sigma]).  If not provided,
)    mu = 0, sigma = 1 are assumed
) For chisquare, parameters = df > 0
) For student, parameters = df > 0
) For f, parameters = vector(df1 > 0, df2 > 0)
) For gamma, parameters = vector(alpha [, mu]), default mu = alpha
) For beta, parameters = vector(alpha, beta)
)) Written by C. Bingham (kb@stat.umn.edu)
)) 991030 version
)) 010113 version; corrected bug in Student's t
# usage: $S(name, x [,parameters])
@usage <- "$S(name, x [,parameters])"

@nargs <- if ($k != 0 || $v == 0){
	$v
}else{
	$v - ($A[$v] == "")
}
if (@nargs < 2 || @nargs > 3){
	error(@usage)
}

@name <- argvalue($01, "name", "string")
@x <- argvalue($02, "x", "nonmissing real vector")
@p <- if (@nargs > 2){
	argvalue($03, "$3", "nonmissing real vector")
}else{
	NULL
}
delete(@nargs)

@iname <- match(@name,\
	vector("normal", "chisquare", "student", "F", "beta","gamma","f"), 0)
if (@iname == 0){
	error(paste("\"",@name,"\""," is not a legal distribution name",sep:""))
}
if (@iname == 7){
	@iname <- 4
}

@all <- keyvalue($K, "all", "TF", default:F)
@xval <- keyvalue($K, "x", "TF", default:@all)
@params <- keyvalue($K, "param*", "TF", default:@all)
@densname <- keyvalue($K, "name", "TF", default:@all)

@name <- paste(@name,"density")
if (@iname == 1){#normal
	if(isnull(@p)){
		@p <- 0
	}
	if (length(@p) == 1){
		@p <- vector(@p, 1)
	}
	@mu <- @p[1]
	@sigma <- @p[2]
	if (@sigma <= 0){
		error("standard deviation not positive")
	}
	@result <- exp(-.5*((@x - @mu)/@sigma)^2)/(@sigma*sqrt(2*PI))
	delete(@mu, @sigma)
}elseif (@iname == 2){#chisquare
	if (isnull(@p)){
		error(paste("degrees of freedom required for",@name))
	}
	@df1 <- @p[1]/2
	if (@df1 <= 0){
		error(paste("degrees of freedom not positive for",@name))
	}
	@result <- 0*@x
	@J <- vector(@x > 0)
	if (sum(@J) > 0){
		@result[@J] <-\
			exp(-@df1*log(2) - lgamma(@df1) - @x[@J]/2 + (@df1 - 1)*log(@x[@J]))
	}
	delete(@df1, @J)
}elseif (@iname == 3){#student
	if (isnull(@p)){
		error(paste("degrees of freedom required for",@name))
	}
	@df1 <- @p[1]
	if (@df1 <= 0){
		error(paste("degrees of freedom not positive for",@name))
	}
	@result <- \
		exp(-(.5*log(@df1)+.5*log(PI)+lgamma(@df1/2)-lgamma((@df1+1)/2) +\
			((@df1+1)/2)*log(1 + @x^2/@df1)))
	delete(@df1)
}elseif (@iname == 4){#F distribution
	if (length(@p) != 2){
		error(paste("df1 and df2 required for",@name))
	}
	if (min(@p) <= 0){
		error(paste("df1 or df2 <= 0 for",@name))
	}
	@df1 <- @p[1]
	@df2 <- @p[2]

	@result <- 0*@x
	@J <- vector(@x > 0)
	if (sum(@J) > 0){
		@result[@J] <-\
			exp(@df1*log(@df1)/2+@df2*log(@df2)/2 +\
			lgamma((@df1 + @df2)/2) - lgamma(@df1/2) - lgamma(@df2/2) +\
			(@df1/2-1)*log(@x[@J])-((@df1 + @df2)/2)*log(@df2 + @df1*@x[@J]))
	}
	delete(@df1, @df2, @J)
}elseif (@iname == 5){#beta
	if (length(@p) != 2){
		error(paste("two shape parameters required for",@name))
	}
	if (min(@p) <= 0){
		error(paste("shape parameters < 0 is illegal for", @name))
	}
	@df1 <- @p[1]
	@df2 <- @p[2]

	@result <- 0*@x
	@J <- vector(@x > 0 && @x < 1)
	if (sum(@J) > 0){
		@result[@J] <- exp(lgamma(@df1+@df2) - lgamma(@df1) - lgamma(@df2) +\
			(@df1 - 1)*log(@x[@J]) + (@df2 - 1)*log(1 - @x[@J]))
	}
	delete(@df1, @df2, @J)
}elseif (@iname == 6){#gamma
	if (isnull(@p)){
		error("shape parameter required for", @name)
	}
	if (length(@p) == 1){
		@p <- rep(@p,2)
	}
	if(min(@p) <= 0){
		error(paste("parameters for",@name,"must be positive"))
	}
	@df1 <- @p[1]
	@theta <- @p[2]/@df1 #@p[2] is mu = E[x]

	@result <- 0*@x
	@J <- vector(@x > 0)
	if (sum(@J) > 0){
		@result[@J] <-\
			exp(-lgamma(@df1) - @df1*log(@theta) - @x[@J]/@theta +\
			(@df1 - 1)*log(@x[@J]))
	}
	delete(@df1, @theta, @J)
}
if (@densname || @params || @xval){
	@result <- structure(density:@result)
	if (@xval){
		@result <- strconcat(x:@x, @result)
	}
	if (@params){
		@result <- strconcat(@result,parameters:@p)
	}
	if (@densname){
		@result <- strconcat(@result, name:@name)
	}
}
delete(@x, @p, @name, @iname)
delete(@result, return:T)
%density%

densnor       MACRO DOLLARS
) densnor(x [,mu [, sigma]] [,keyword phrases])
)  computes normal density
) Keyword phrases all:T, x:T or F,parameters:T or F, name:T or F
) Without keyword phrases, returns density in the same size and
) shape as x
) Otherwise returns structure([x ,] density [,parameters] [,name])
)) Written by C. Bingham (kb@stat.umn.edu)
)) 991030 version
# usage: $S(x [,mu [,sigma]])
if ($v == 0){
	error("usage: $S(x [,mu [,sigma]])")
}
@mu <- if($v > 1){
	argvalue($02,"mu","number")
}else{
	0
}
@sigma <- if($v > 2){
	argvalue($03,"sigma", "positive scalar")
}else{
	1
}
if (!ismacro(density)){
	getmacros(density, quiet:T)
	if (!ismacro(density)){
		error("can't locate macro \"density\"")
	}
}
density("normal",$01,vector(@mu,@sigma), $K)
%densnor%

densstu      MACRO DOLLARS
) densstu(x, df [,keyword phrases])
)  computes density of Student's t
) Keyword phrases all:T, x:T or F,parameters:T or F, name:T or F
) Without keyword phrases, returns density in the same size and
) shape as x
) Otherwise returns structure([x ,] density [,parameters] [,name])
)) Written by C. Bingham (kb@stat.umn.edu)
)) 991030 version
# usage: $S(x , df)
if ($v != 2){
	error("usage: $S(x, df")
}
@df <- argvalue($02,"degrees of freedom","positive number")
if (!ismacro(density)){
	getmacros(density, quiet:T)
	if (!ismacro(density)){
		error("can't locate macro \"density\"")
	}
}

density("student",$01,@df, $K)
%densstu%'

denschi     MACRO DOLLARS
) denschi(x, df [,keyword phrases])
)   computes density of chi-squared distribution
) Keyword phrases all:T, x:T or F,parameters:T or F, name:T or F
) Without keyword phrases, returns density in the same size and
) shape as x
) Otherwise returns structure([x ,] density [,parameters] [,name])
)) Written by C. Bingham (kb@stat.umn.edu)
)) 991030 version
# usage: $S(x , df)
if ($v != 2){
	error("usage: $S(x, df")
}
@df <- argvalue($02,"degrees of freedom","positive number")
if (!ismacro(density)){
	getmacros(density, quiet:T)
	if (!ismacro(density)){
		error("can't locate macro \"density\"")
	}
}

density("chisquare",$01,@df, $K)
%denschi%

densF     MACRO DOLLARS
) densF(x, df1, df2 [, keyword phrases])
)  computes density of F distribution
) Keyword phrases all:T, x:T or F,parameters:T or F, name:T or F
) Without keyword phrases, returns density in the same size and
) shape as x
) Otherwise returns structure([x ,] density [,params] [,name])
)) Written by C. Bingham (kb@stat.umn.edu)
)) 991030 version
# usage: $S(x , df1, df2 [, keyword phrases])
if ($v != 3){
	error("usage: $S(x, df1, df2")
}
@df1 <- argvalue($02,"numerator d.f.","positive number")
@df2 <- argvalue($03,"denominator d.f.","positive number")
if (!ismacro(density)){
	getmacros(density, quiet:T)
	if (!ismacro(density)){
		error("can't locate macro \"density\"")
	}
}

density("F",$01,vector(@df1, @df2), $K)
%densF%

densbeta  MACRO DOLLARS
) densbeta(x,nu1,nu2 [,keyword phrases])
)   computes density of beta distribution
) Keyword phrases all:T, x:T or F,parameters:T or F, name:T or F
) Without keyword phrases, returns density in the same size and
) shape as x
) Otherwise returns structure([x ,] density [,params] [,name])
) Density is x^(nu1-1)*(1-x)^(nu2-1))/beta(nu1,nu2)
)) Written by C. Bingham (kb@stat.umn.edu)
)) 991030 version
# usage: $S(x , nu1, nu2)
if ($v != 3){
	error("usage: $S(x, nu1, nu2")
}
@df1 <- argvalue($02,"$2","positive number")
@df2 <- argvalue($03,"$3","positive number")
if (!ismacro(density)){
	getmacros(density, quiet:T)
	if (!ismacro(density)){
		error("can't locate macro \"density\"")
	}
}
if (!ismacro(density)){
	getmacros(density, quiet:T)
	if (!ismacro(density)){
		error("can't locate macro \"density\"")
	}
}

density("beta",$01,vector(@df1, @df2), $K)
%densbeta%

densgamma MACRO DOLLARS
) densgamma(x, alpha [, mu] [,keyword phrases])
)  computes density of gamma distribution
) Keyword phrases all:T, x:T or F,parameters:T or F, name:T or F
) Without keyword phrases, returns density in the same size and
) shape as x
) Otherwise returns structure([x ,] density [,params] [,name])
) Density is
)   ((alpha/mu)^alpha/gamma(alpha)) * exp(-(alpha/mu)*x)*x^(alpha-1)
) Default for mu is alpha
)) Written by C. Bingham (kb@stat.umn.edu)
)) 991030 version
# usage: $S(x , alpha [, mu])
if ($v < 2 || $v > 3){
	error("usage: $S(x, alpha [, mu])")
}
@df1 <- argvalue($02,"shape parameter","positive number")
if ($v == 3){
	@df1 <- vector(@df1, argvalue($03,"mu","positive number"))
}
if (!ismacro(density)){
	getmacros(density, quiet:T)
	if (!ismacro(density)){
		error("can't locate macro \"density\"")
	}
}

density("gamma",$01,@df1, $K)
%densgamma%