densityest   MACRO DOLLARS
) Macro to compute a density estimate by the kernel estimate
) f(x) = sum((W((x - y[i])/h)/h),i=1,...,n) where W(z) is a
) density and y is a data vector of length n
) Usage:
)   densy <- densityest(y,x [,h])
)     y           REAL vector of data]
)     x           REAL vector with no MISSING values at which density
)                 of y is to be estimated
)     h           positive scalar, default h = sy/4, sy = SD of y
)     densy       structure(x:x,density:densityEst,h:h)
)     densEst     REAL vector with densEst[j] = estimated density
)                 at x[j], using W(z) = exp(-.5*z^2)/sqrt(2*PI)
)   densy <- densityest(y, x [, h] , window:W)
)     W           macro such that W(z) defines a density
)) C. Bingham (kb@stat.umn.edu)
)) 020304
#$S(data,x [,h] [,window:windMacro])
@data <- argvalue($1,"data","real vector")
@x <- argvalue($2,"x","real vector nonmissing")
if (anymissing(@data)){
	@data <- @data[!ismissing(@data)]
	if (isnull(@data)){
		error("no non-MISSING data")
	}
}
@h <- if ($v <= 2){
	.25*describe(@data,stddev:T)
} else {
	argvalue($3,"h","positive number")
}
@wind <- keyvalue($K,"window","macro")
if (isnull(@wind)){
	@wind <- macro("exp(-.5*(\$1)^2)/sqrt(2*PI)")
}
@density <- 0*@x
@n <- length(@data)

for (@i,1,@n){
	@density <-+ @wind((@x - @data[@i])/@h)
}
@density <-/ @n*@h
delete(@data,@wind,@i)
structure(x:delete(@x,return:T),density:delete(@density,return:T),\
		h:delete(@h,return:T))
%densityest%

_E_N_D_O_F_M_A_C_R_O_S_

Help file for densityest.mac for MacAnova
(C) 2003 by Gary W. Oehlert and Christopher Bingham
Updated 030204 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
Density
???? Ending marker for keys

030204 created

====densityest()#density
%%%%
densy <- densityest(y,xvals [,h]), REAL vector y of data, REAL vector
  xvals of value where density is to be estimated, optional REAL
  scalar h > 0 standard deviation of Normal kernel function.  densy is
  structure(x:xvals,density:densityhat, h:h)
%%%%
@@@@introduction
densityest() estimates a density from a random zample using the kernal
estimate f_hat(x) = sum(W((x-y[i])/h))/(n*h) where W(x) is a
probability density.  The default for W(x) is the standard normal
density.

@@@@usage
densy <- densityest(y,xvals) computes an estimated density based on data
in y at values xval.  y and xvals are REAL vectors.  Often xvals will
have the form run(xmin,xmax,step).

The result, densy, is structure(x:xvals, density:densityhat, h:h) where
densityhat is a vector with densityhat[i] = estimated density at
xvals[i] and h is the default smoothing parameter = SD[y]/4.

densy <- densityest(y,xvals,h) where h > 0 is a REAL scalar, does the
same except h is the smoothing parameter instead of SD[y]/4.

densy <- densityest(y,x [,h], window:W) where W() is a macro such that
W(x) defines a density does the same, except W(x) is used instead of the
standard normal density.

@@@@graphing_density
You can directly graph the estimated density by, for example

  Cmd> lineplot(densityest(y, x),title:"Estimated density")
  
You can make a histogram overlaid with a density estimate by, for example,

  Cmd> hist(y,vector(0,15), show:F) # histogram not yet displayed
  
  Cmd> addlines(densityest(y,run(0,75,2.5)))

@@@@______

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