# R/density.fd.R In fda: Functional Data Analysis

#### Documented in density.fd

```density <- function(x, ...)UseMethod('density')

density.fd <- function(x, WfdParobj, conv=0.0001, iterlim=20,
active=1:nbasis, dbglev=0, ...) {
# DENSITYFD estimates the density of a sample of scalar observations.

#  These observations may be one of two forms:
#   1.  a vector of observatons x_i
#   2.  a two-column matrix, with the observations x_i in the
#       first column, and frequencies f_i in the second.
#   Option 1. corresponds to all f_i = 1.

#  Arguments are:
#  X         ... data value array, either a vector or a two-column
#                matrix.
#  WFDPAROBJ ... functional parameter object specifying the initial log
#              density, the linear differential operator used to smooth
#              smooth it, and the smoothing parameter.
#  CONV      ... convergence criterion
#  ITERLIM   ... iteration limit for scoring iterations
#  ACTIVE    ... indices among 1:NBASIS of parameters to optimize
#  DBGLEV    ... level of output of computation history

#  Returns:
#  A list containing
#  WFDOBJ ...   functional data basis object defining final density
#  C      ...   normalizing constant for density p = exp(WFDOBJ)/C
#  FLIST  ...   Struct object containing
#               FSTR\$f     final log likelihood
#               FSTR\$norm  final norm of gradient
#  ITERNUM   Number of iterations
#  ITERHIST  History of iterations

#  To plot the density function or to evaluate it, evaluate WFDOBJ,
#  exponentiate the resulting vector, and then divide by the normalizing
#  constant C.

#  check WfdParobj

if (!inherits(WfdParobj, "fdPar")) {
if (inherits(WfdParobj, "fd") || inherits(WfdParobj, "basisfd")) {
WfdParobj <- fdPar(WfdParobj)
} else {
stop("WFDPAROBJ is not a fdPar object")
}
}

#  set up WFDOBJ

Wfdobj   <- WfdParobj\$fd

#  set up LFDOBJ

Lfdobj <- WfdParobj\$Lfd
Lfdobj <- int2Lfd(Lfdobj)

#  set up BASIS

basisobj <- Wfdobj\$basis
nbasis   <- basisobj\$nbasis
rangex   <- basisobj\$rangeval

x    <- as.matrix(x)
xdim <- dim(x)
N    <- xdim[1]
m    <- xdim[2]

if (m > 2 && N > 2)
stop("Argument X must have either one or two columns.")

if ((N == 1 | N == 2) & m > 1) {
x <- t(x)
n <- N
N <- m
m <- n
}

if (m == 1) {
f <- rep(1,N)
} else {
f    <- x[,2]
fsum <- sum(f)
f    <- f/fsum
x    <- x[,1]
}
f = as.matrix(f)

#  check for values outside of the range of WFD0

inrng <- (1:N)[x >= rangex[1] & x <= rangex[2]]
if (length(inrng) != N) {
print(c(length(inrng), N))
print(c(rangex[1], rangex[2], min(x), max(x)))
warning("Some values in X out of range and not used.")
}
x     <- x[inrng]
f     <- f[inrng]
nobs  <- length(x)

#  set up some arrays

climit <- c(rep(-50,nbasis),rep(400,nbasis))
cvec0  <- Wfdobj\$coefs
dbgwrd <- dbglev > 1

zeromat <- zerobasis(nbasis)

#  initialize matrix Kmat defining penalty term

lambda <- WfdParobj\$lambda
if (lambda > 0) Kmat <- lambda*getbasispenalty(basisobj, Lfdobj)

#  evaluate log likelihood
#    and its derivatives with respect to these coefficients

result <- loglfnden(x, f, basisobj, cvec0)
logl   <- result[[1]]
Dlogl  <- result[[2]]

#  compute initial badness of fit measures

fun  <- -logl
gvec <- -Dlogl
if (lambda > 0) {
gvec <- gvec + 2*(Kmat %*% cvec0)
fun    <- fun + t(cvec0) %*% Kmat %*% cvec0
}
Foldstr <- list(f = fun, norm = sqrt(mean(gvec^2)))
gvec0 <- t(zeromat) %*% as.matrix(gvec)

#  compute the initial expected Hessian

hmat <- Varfnden(x, basisobj, cvec0)
if (lambda > 0) hmat <- hmat + 2*Kmat
hmat0 = t(zeromat) %*% hmat %*% zeromat

#  evaluate the initial update vector for correcting the initial bmat

deltac0  <- -solve(hmat0,gvec0)
deltac   <- zeromat %*% as.matrix(deltac0)
cosangle <- -sum(gvec0*deltac0)/sqrt(sum(gvec0^2)*sum(deltac0^2))

#  initialize iteration status arrays

iternum <- 0
status <- c(iternum, Foldstr\$f, -logl, Foldstr\$norm)
if (dbglev > 0) {
cat("Iteration  Criterion  Neg. Log L  Grad. Norm\n")
cat("      ")
cat(format(iternum))
cat("    ")
cat(format(status[2:4]))
cat("\n")
}
iterhist <- matrix(0,iterlim+1,length(status))
iterhist[1,]  <- status

#  quit if ITERLIM == 0

if (iterlim == 0) {
Flist     <- Foldstr
iterhist <- iterhist[1,]
C        <- normden.phi(basisobj, cvec0)
return( list(Wfdobj=Wfdobj, C=C, Flist=Flist, iternum=iternum,
iterhist=iterhist) )
}

#  -------  Begin iterations  -----------

STEPMAX <- 5
MAXSTEP <- 400
trial   <- 1
cvec    <- cvec0
linemat <- matrix(0,3,5)

for (iter in 1:iterlim) {
iternum <- iternum + 1
#  take optimal stepsize
dblwrd <- rep(FALSE,2)
limwrd <- rep(FALSE,2)
stpwrd <- 0
ind    <- 0
#  compute slope
Flist <- Foldstr
linemat[2,1] <- sum(deltac*gvec)
#  normalize search direction vector
sdg     <- sqrt(sum(deltac^2))
deltac  <- deltac/sdg
dgsum   <- sum(deltac)
linemat[2,1] <- linemat[2,1]/sdg
#  return with stop condition if (initial slope is nonnegative
if (linemat[2,1] >= 0) {
print("Initial slope nonnegative.")
ind <- 3
iterhist <- iterhist[1:(iternum+1),]
break
}
#  return successfully if (initial slope is very small
if (linemat[2,1] >= -1e-5) {
if (dbglev>1) print("Initial slope too small")
iterhist <- iterhist[1:(iternum+1),]
break
}
#  load up initial search matrix
linemat[1,1:4] <- 0
linemat[2,1:4] <- linemat[2,1]
linemat[3,1:4] <- Foldstr\$f
#  output initial results for stepsize 0
stepiter  <- 0
if (dbglev > 1) {
cat("              ")
cat(format(stepiter))
cat(format(linemat[,1]))
cat("\n")
}
ips <- 0
#  first step set to trial
linemat[1,5]  <- trial
#  Main iteration loop for linesrch
for (stepiter in 1:STEPMAX) {
#  ensure that step does not go beyond limits on parameters
limflg  <- 0
#  check the step size
result <- stepchk(linemat[1,5], cvec, deltac, limwrd, ind,
climit, active, dbgwrd)
linemat[1,5] <- result[[1]]
ind          <- result[[2]]
limwrd       <- result[[3]]
if (linemat[1,5] <= 1e-9) {
#  Current step size too small  terminate
Flist   <- Foldstr
cvecnew <- cvec
gvecnew <- gvec
if (dbglev > 1) print(paste("Stepsize too small:", linemat[1,5]))
if (limflg) ind <- 1 else ind <- 4
break
}
cvecnew <- cvec + linemat[1,5]*deltac
#  compute new function value and gradient
result  <- loglfnden(x, f, basisobj, cvecnew)
logl    <- result[[1]]
Dlogl   <- result[[2]]
Flist\$f <- -logl
gvecnew <- -as.matrix(Dlogl)
if (lambda > 0) {
gvecnew <- gvecnew + 2*Kmat %*% cvecnew
Flist\$f <- Flist\$f + t(cvecnew) %*% Kmat %*% cvecnew
}
gvecnew0 <- t(zeromat) %*% gvecnew
Flist\$norm <- sqrt(mean(gvecnew0^2))
#  compute new directional derivative
linemat[2,5] <- sum(deltac*gvecnew)
linemat[3,5] <- Flist\$f
if (dbglev > 1) {
cat("              ")
cat(format(stepiter))
cat(format(linemat[,5]))
cat("\n")
}
#  compute next step
result  <- stepit(linemat, ips, dblwrd, MAXSTEP)
linemat <- result\$linemat
ips     <- result\$ips
ind     <- result\$ind
dblwrd  <- result\$dblwrd
trial   <- linemat[1,5]
#  ind == 0 implies convergence
if (ind == 0 | ind == 5) break
#  end of line search loop
}

#  update current parameter vectors

cvec  <- cvecnew
gvec  <- gvecnew
gvec0 <- t(zeromat) %*% as.matrix(gvec)
Wfdobj\$coefs <- cvec
status <- c(iternum, Flist\$f, -logl, Flist\$norm)
iterhist[iter+1,] <- status
if (dbglev > 0) {
cat("      ")
cat(format(iternum))
cat("    ")
cat(format(status[2:4]))
cat("\n")
}

#  test for convergence

if (abs(Flist\$f-Foldstr\$f) < conv) {
iterhist <- iterhist[1:(iternum+1),]
C <- normden.phi(basisobj, cvec)
denslist <- list("Wfdobj" = Wfdobj, "C" = C, "Flist" = Flist,
"iternum" = iternum, "iterhist" = iterhist)
return( denslist )
}
if (Flist\$f >= Foldstr\$f) break
#  compute the Hessian
hmat <- Varfnden(x, basisobj, cvec)
if (lambda > 0) hmat <- hmat + 2*Kmat
hmat0 <- t(zeromat) %*% hmat %*% zeromat
#  evaluate the update vector
deltac0 <- -solve(hmat0,gvec0)
cosangle  <- -sum(gvec0*deltac0)/sqrt(sum(gvec0^2)*sum(deltac0^2))
if (cosangle < 0) {
if (dbglev > 1) print("cos(angle) negative")
deltac0 <- -gvec0
}
deltac <- zeromat %*% as.matrix(deltac0)
Foldstr <- Flist
#  end of iterations
}
#  compute final normalizing constant
C <- normden.phi(basisobj, cvec)
denslist <- list("Wfdobj" = Wfdobj, "C" = C, "Flist" = Flist,
"iternum" = iternum, "iterhist" = iterhist)
return( denslist )
}

#  -----------------------------------------------------------------------------

loglfnden <- function(x, f, basisobj, cvec=FALSE) {
#  Computes the log likelihood and its derivative with
#    respect to the coefficients in CVEC
N       <- length(x)
nbasis  <- basisobj\$nbasis
fmat    <- outer(f, rep(1,nbasis))
fsum    <- sum(f)
nobs    <- length(x)
phimat  <- getbasismatrix(x, basisobj)
Cval    <- normden.phi(basisobj, cvec, )
logl    <- sum((phimat %*% cvec) * f - fsum*log(Cval)/N)
EDw     <- expectden.phi(basisobj, cvec, Cval)
Dlogl   <- apply((phimat - outer(rep(1,nobs),EDw))*fmat,2,sum)
return( list(logl, Dlogl) )
}

#  -----------------------------------------------------------------------------

Varfnden <- function(x, basisobj, cvec=FALSE) {
#  Computes the expected Hessian
nbasis  <- basisobj\$nbasis
nobs    <- length(x)
Cval    <- normden.phi(basisobj, cvec)
EDw     <- expectden.phi(basisobj, cvec, Cval)
EDwDwt  <- expectden.phiphit(basisobj, cvec, Cval)
Varphi  <- nobs*(EDwDwt - outer(EDw,EDw))
return(Varphi)
}

#  -----------------------------------------------------------------------------

normden.phi <- function(basisobj, cvec, JMAX=15, EPS=1e-7) {

#  Computes integrals of
#      p(x) = exp phi"(x) %*% cvec
#  by numerical integration using Romberg integration

#  check arguments, and convert basis objects to functional data objects

if (!inherits(basisobj, "basisfd") )
stop("First argument must be a basis function object.")

nbasis <- basisobj\$nbasis
oneb   <- matrix(1,1,nbasis)
rng    <- basisobj\$rangeval

#  set up first iteration

width <- rng[2] - rng[1]
JMAXP <- JMAX + 1
h <- matrix(1,JMAXP,1)
h[2] <- 0.25
#  matrix SMAT contains history of discrete approximations to the integral
smat <- matrix(0,JMAXP,1)
#  the first iteration uses just the }points
x  <- rng
nx <- length(x)
ox <- matrix(1,nx,1)
fx <- getbasismatrix(x, basisobj)
wx <- fx %*% cvec
wx[wx < -50] <- -50
px <- exp(wx)
smat <- matrix(0,JMAXP,1)
smat[1]  <- width*sum(px)/2
tnm <- 0.5
j   <- 1

#  now iterate to convergence
for (j in 2:JMAX) {
tnm  <- tnm*2
del  <- width/tnm
if (j == 2) {
x <- (rng[1] + rng[2])/2
} else {
x <- seq(rng[1]+del/2, rng[2], del)
}
fx <- getbasismatrix(x, basisobj)
wx <- fx %*% cvec
wx[wx < -50] <- -50
px <- exp(wx)
smat[j] <- (smat[j-1] + width*sum(px)/tnm)/2
if (j >= 5) {
ind <- (j-4):j
result <- polintarray(h[ind],smat[ind],0)
ss  <- result[[1]]
dss <- result[[2]]
if (!any(abs(dss) >= EPS*max(abs(ss)))) {
#  successful convergence
return(ss)
}
}
smat[j+1] <- smat[j]
h[j+1]    <- 0.25*h[j]
}
warning(paste("No convergence after ",JMAX," steps in NORMDEN.PHI"))
return(ss)
}

#  -----------------------------------------------------------------------------

expectden.phi <- function(basisobj, cvec, Cval=1, nderiv=0,
JMAX=15, EPS=1e-7) {
#  Computes expectations of basis functions with respect to density
#      p(x) <- Cval^{-1} exp t(c)*phi(x)
#  by numerical integration using Romberg integration

#  check arguments, and convert basis objects to functional data objects

if (!inherits(basisobj, "basisfd"))
stop("First argument must be a basis function object.")

nbasis <- basisobj\$nbasis
rng    <- basisobj\$rangeval
oneb   <- matrix(1,1,nbasis)

#  set up first iteration

width <- rng[2] - rng[1]
JMAXP <- JMAX + 1
h <- matrix(1,JMAXP,1)
h[2] <- 0.25
#  matrix SMAT contains the history of discrete approximations to the integral
sumj <- matrix(0,1,nbasis)
#  the first iteration uses just the }points
x  <- rng
nx <- length(x)
ox <- matrix(1,nx,nx)
fx <- getbasismatrix(x, basisobj, 0)
wx <- fx %*% cvec
wx[wx < -50] <- -50
px <- exp(wx)/Cval
if (nderiv == 0) {
Dfx <- fx
} else {
Dfx <- getbasismatrix(x, basisobj, 1)
}
sumj <- t(Dfx) %*% px
smat <- matrix(0,JMAXP,nbasis)
smat[1,]  <- width*as.vector(sumj)/2
tnm <- 0.5
j   <- 1

#  now iterate to convergence

for (j in 2:JMAX) {
tnm  <- tnm*2
del  <- width/tnm
if (j == 2) {
x <- (rng[1] + rng[2])/2
} else {
x <- seq(rng[1]+del/2, rng[2], del)
}
nx <- length(x)
fx <- getbasismatrix(x, basisobj, 0)
wx <- fx %*% cvec
wx[wx < -50] <- -50
px <- exp(wx)/Cval
if (nderiv == 0) {
Dfx <- fx
} else {
Dfx <- getbasismatrix(x, basisobj, 1)
}
sumj <- t(Dfx) %*% px
smat[j,] <- (smat[j-1,] + width*as.vector(sumj)/tnm)/2
if (j >= 5) {
ind <- (j-4):j
temp <- smat[ind,]
result <- polintarray(h[ind],temp,0)
ss  <- result[[1]]
dss <- result[[2]]
if (!any(abs(dss) > EPS*max(abs(ss)))) {
#  successful convergence
return(ss)
}
}
smat[j+1,] <- smat[j,]
h[j+1]     <- 0.25*h[j]
}
warning(paste("No convergence after ",JMAX," steps in EXPECTDEN.PHI"))
return(ss)
}

#  -----------------------------------------------------------------------------

expectden.phiphit <- function(basisobj, cvec, Cval=1, nderiv1=0, nderiv2=0,
JMAX=15, EPS=1e-7) {

#  Computes expectations of cross product of basis functions with
#  respect to density
#      p(x) = Cval^{-1} int [exp t(c) %*% phi(x)] phi(x) t(phi(x)) dx
#  by numerical integration using Romberg integration

#  check arguments, and convert basis objects to functional data objects

if (!inherits(basisobj, "basisfd"))
stop("First argument must be a basis function object.")

nbasis <- basisobj\$nbasis
rng    <- basisobj\$rangeval
oneb   <- matrix(1,1,nbasis)

#  set up first iteration

width <- rng[2] - rng[1]
JMAXP <- JMAX + 1
h <- matrix(1,JMAXP,1)
h[2] <- 0.25
#  matrix SMAT contains history of discrete approximations to the integral
#  the first iteration uses just the }points
x  <- rng
nx <- length(x)
fx <- getbasismatrix(x, basisobj, 0)
wx <- fx %*% cvec
wx[wx < -50] <- -50
px <- exp(wx)/Cval
if (nderiv1 == 0) {
Dfx1 <- fx
} else {
Dfx1 <- getbasismatrix(x, basisobj, 1)
}
if (nderiv2 == 0) {
Dfx2 <- fx
} else {
Dfx2 <- getbasismatrix(x, basisobj, 2)
}
oneb <- matrix(1,1,nbasis)
sumj <- t(Dfx1) %*% ((px %*% oneb) * Dfx2)
smat <- array(0,c(JMAXP,nbasis,nbasis))
smat[1,,]  <- width*as.matrix(sumj)/2
tnm <- 0.5
j   <- 1

#  now iterate to convergence
for (j in 2:JMAX) {
tnm  <- tnm*2
del  <- width/tnm
if (j == 2) {
x <- (rng[1] + rng[2])/2
} else {
x <- seq(rng[1]+del/2, rng[2], del)
}
nx <- length(x)
fx <- getbasismatrix(x, basisobj, 0)
wx <- fx %*% cvec
wx[wx < -50] <- -50
px <- exp(wx)/Cval
if (nderiv1 == 0) {
Dfx1 <- fx
} else {
Dfx1 <- getbasismatrix(x, basisobj, 1)
}
if (nderiv2 == 0) {
Dfx2 <- fx
} else {
Dfx2 <- getbasismatrix(x, basisobj, 2)
}
sumj <- t(Dfx1) %*% ((px %*% oneb) * Dfx2)
smat[j,,] <- (smat[j-1,,] + width*as.matrix(sumj)/tnm)/2
if (j >= 5) {
ind <- (j-4):j
temp <- smat[ind,,]
result <- polintarray(h[ind],temp,0)
ss  <- result[[1]]
dss <- result[[2]]
if (!any(abs(dss) > EPS*max(max(abs(ss))))) {
#  successful convergence
return(ss)
}
}
smat[j+1,,] <- smat[j,,]
h[j+1] <- 0.25*h[j]
}
warning(paste("No convergence after ",JMAX," steps in EXPECTDEN.PHIPHIT"))
return(ss)
}

#  -----------------------------------------------------------------------------

polintarray <- function(xa, ya, x0) {
#  YA is an array with up to 4 dimensions
#     with 1st dim the same length same as the vector XA
n     <- length(xa)
nydim <- 1
} else {
}
if (yadim[1] != n) stop("First dimension of YA must match XA")
difx <- xa - x0
absxmxa <- abs(difx)
ns <- min((1:n)[absxmxa == min(absxmxa)])
cs <- ya
ds <- ya
if (nydim == 1) y <- ya[ns]
if (nydim == 2) y <- ya[ns,]
if (nydim == 3) y <- ya[ns,,]
if (nydim == 4) y <- ya[ns,,,]
ns <- ns - 1
for (m in 1:(n-1)) {
if (nydim == 1) {
for (i in 1:(n-m)) {
ho    <- difx[i]
hp    <- difx[i+m]
w     <- (cs[i+1] - ds[i])/(ho - hp)
ds[i] <- hp*w
cs[i] <- ho*w
}
if (2*ns < n-m) {
dy <- cs[ns+1]
} else {
dy <- ds[ns]
ns <- ns - 1
}
}
if (nydim == 2) {
for (i in 1:(n-m)) {
ho     <- difx[i]
hp     <- difx[i+m]
w      <- (cs[i+1,] - ds[i,])/(ho - hp)
ds[i,] <- hp*w
cs[i,] <- ho*w
}
if (2*ns < n-m) {
dy <- cs[ns+1,]
} else {
dy <- ds[ns,]
ns <- ns - 1
}
}
if (nydim == 3) {
for (i in 1:(n-m)) {
ho       <- difx[i]
hp       <- difx[i+m]
w        <- (cs[i+1,,] - ds[i,,])/(ho - hp)
ds[i,,] <- hp*w
cs[i,,] <- ho*w
}
if (2*ns < n-m) {
dy <- cs[ns+1,,]
} else {
dy <- ds[ns,,]
ns <- ns - 1
}
}
if (nydim == 4) {
for (i in 1:(n-m)) {
ho      <- difx[i]
hp      <- difx[i+m]
w       <- (cs[i+1,,,] - ds[i,,,])/(ho - hp)
ds[i,,,] <- hp*w
cs[i,,,] <- ho*w
}
if (2*ns < n-m) {
dy <- cs[ns+1,,,]

} else {
dy <- ds[ns,,,]
ns <- ns - 1
}
}
y <- y + dy
}
return( list(y, dy) )
}
```

## Try the fda package in your browser

Any scripts or data that you put into this service are public.

fda documentation built on May 29, 2024, 11:26 a.m.