fda (version 6.2.0)

density.fd: Compute a Probability Density Function

Description

Like the regular S-PLUS function density, this function computes a probability density function for a sample of values of a random variable. However, in this case the density function is defined by a functional parameter object WfdParobj along with a normalizing constant C.

The density function $p(x)$ has the form p(x) = C exp[W(x)] where function $W(x)$ is defined by the functional data object WfdParobj.

Usage

# S3 method for fd
density(x, WfdParobj, conv=0.0001, iterlim=20,
           active=1:nbasis, dbglev=0, ...)

Value

a named list of length 4 containing:

Wfdobj

a functional data object defining function $W(x)$ that that optimizes the fit to the data of the monotone function that it defines.

C

the normalizing constant.

Flist

a named list containing three results for the final converged solution: (1) f: the optimal function value being minimized, (2) grad: the gradient vector at the optimal solution, and (3) norm: the norm of the gradient vector at the optimal solution.

iternum

the number of iterations.

iterhist

a iternum+1 by 5 matrix containing the iteration history.

Arguments

x

a set observations, which may be one of two forms:

1. a vector of observations $x_i$

2. a two-column matrix, with the observations $x_i$ in the first column, and frequencies $f_i$ in the second.

The first option corresponds to all $f_i = 1$.

WfdParobj

a functional parameter object specifying the initial value, basis object, roughness penalty and smoothing parameter defining function $W(t).$

conv

a positive constant defining the convergence criterion.

iterlim

the maximum number of iterations allowed.

active

a logical vector of length equal to the number of coefficients defining Wfdobj. If an entry is TRUE, the corresponding coefficient is estimated, and if FALSE, it is held at the value defining the argument Wfdobj. Normally the first coefficient is set to 0 and not estimated, since it is assumed that $W(0) = 0$.

dbglev

either 0, 1, or 2. This controls the amount information printed out on each iteration, with 0 implying no output, 1 intermediate output level, and 2 full output. If levels 1 and 2 are used, it is helpful to turn off the output buffering option in S-PLUS.

...

Other arguments to match the generic function 'density'

Details

The goal of the function is provide a smooth density function estimate that approaches some target density by an amount that is controlled by the linear differential operator Lfdobj and the penalty parameter. For example, if the second derivative of $W(t)$ is penalized heavily, this will force the function to approach a straight line, which in turn will force the density function itself to be nearly normal or Gaussian. Similarly, to each textbook density function there corresponds a $W(t)$, and to each of these in turn their corresponds a linear differential operator that will, when apply to $W(t)$, produce zero as a result. To plot the density function or to evaluate it, evaluate Wfdobj, exponentiate the resulting vector, and then divide by the normalizing constant C.

References

Ramsay, James O., Hooker, Giles, and Graves, Spencer (2009), Functional data analysis with R and Matlab, Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2005), Functional Data Analysis, 2nd ed., Springer, New York.

Ramsay, James O., and Silverman, Bernard W. (2002), Applied Functional Data Analysis, Springer, New York.

See Also

intensity.fd

Examples

#  set up range for density
rangeval <- c(-3,3)
#  set up some standard normal data
x <- rnorm(50)
#  make sure values within the range
x[x < -3] <- -2.99
x[x >  3] <-  2.99
#  set up basis for W(x)
basisobj <- create.bspline.basis(rangeval, 11)
#  set up initial value for Wfdobj
Wfd0 <- fd(matrix(0,11,1), basisobj)
WfdParobj <- fdPar(Wfd0)
#  estimate density
denslist <- density.fd(x, WfdParobj)
#  plot density
oldpar <- par(no.readonly=TRUE)
xval <- seq(-3,3,.2)
wval <- eval.fd(xval, denslist$Wfdobj)
pval <- exp(wval)/denslist$C
plot(xval, pval, type="l", ylim=c(0,0.4))
points(x,rep(0,50))
par(oldpar)