density.fd <- function(x, WfdParobj, conv=0.0001, iterlim=20, active=2:nbasis, dbglev=1, ...) { # 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. # last modified 3 January 2008 by Jim Ramsay # Previously modified 2007 May 3 by Spencer Graves # 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] } # 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 hmat <- matrix(0,nbasis,nbasis) inactive <- rep(TRUE,nbasis) inactive[active] <- FALSE ninactive <- length((1:nbasis)[inactive]) dbgwrd <- dbglev > 1 # 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 f0 <- -logl gvec0 <- -Dlogl if (lambda > 0) { gvec0 <- gvec0 + 2*(Kmat %*% cvec0) f0 <- f0 + t(cvec0) %*% Kmat %*% cvec0 } if (ninactive > 0) gvec0[inactive] <- 0 Foldstr <- list(f = f0, norm = sqrt(mean(gvec0^2))) # compute the initial expected Hessian hmat0 <- Varfnden(x, basisobj, cvec0) if (lambda > 0) hmat0 <- hmat0 + 2*Kmat if (ninactive > 0) { hmat0[inactive,] <- 0 hmat0[,inactive] <- 0 hmat0[inactive,inactive] <- diag(rep(1,ninactive)) } # evaluate the initial update vector for correcting the initial bmat deltac <- -solve(hmat0,gvec0) cosangle <- -sum(gvec0*deltac)/sqrt(sum(gvec0^2)*sum(deltac^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) ) } else { gvec <- gvec0 hmat <- hmat0 } # ------- 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 <- c(0,0) limwrd <- c(0,0) 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 <- -Dlogl if (lambda > 0) { gvecnew <- gvecnew + 2*Kmat %*% cvecnew Flist$f <- Flist$f + t(cvecnew) %*% Kmat %*% cvecnew } if (ninactive > 0) gvecnew[inactive] <- 0 Flist$norm <- sqrt(mean(gvecnew^2)) linemat[3,5] <- Flist$f # compute new directional derivative linemat[2,5] <- sum(deltac*gvecnew) if (dbglev > 1) { cat(" ") cat(format(stepiter)) cat(format(linemat[,1])) cat("\n") } # compute next step result <- stepit(linemat, ips, dblwrd, MAXSTEP) linemat <- result[[1]] ips <- result[[2]] ind <- result[[3]] dblwrd <- result[[4]] 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 Wfdobj$coefs <- cvec status <- c(iternum, Flist$f, -logl, Flist$norm) iterhist[iter+1,] <- status 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 if (ninactive > 0) { hmat[inactive,] <- 0 hmat[,inactive] <- 0 hmat[inactive,inactive] <- diag(rep(1,ninactive)) } # evaluate the update vector deltac <- -solve(hmat,gvec) cosangle <- -sum(gvec*deltac)/sqrt(sum(gvec^2)*sum(deltac^2)) if (cosangle < 0) { if (dbglev > 1) print("cos(angle) negative") deltac <- -gvec } 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) { # 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 <- outer(rep(1,nobs),expectden.phi(basisobj, cvec, cval)) Dlogl <- apply((phimat - EDW)*fmat,2,sum) return( list(logl, Dlogl) ) } # --------------------------------------------------------------- Varfnden <- function(x, basisobj, cvec) { # Computes the expected Hessian nbasis <- basisobj$nbasis nobs <- length(x) cval <- normden.phi(basisobj, cvec) EDw <- outer(rep(1,nobs),expectden.phi(basisobj, cvec, cval)) EDwDwt <- nobs*expectden.phiphit(basisobj, cvec, cval) Varphi <- EDwDwt - crossprod(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 # Arguments: # BASISOBJ ... Basis function object with basis functions phi. # CVEC ... coefficient vector defining density, of length NBASIS # MU ... mean values to be subtracted from variates # SIGMA ... standard deviation to define u = (x - mu)/sigma # RNG ... vector of length 2 giving the interval over which the # integration is to take place. Multiply a standard interval # like (-5,5) by sigma to make it scale free # JMAX ... maximum number of allowable iterations # EPS ... convergence criterion for relative stop # Return: # The integral of the function. # 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 the 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[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, rng=rangeval, 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 # Arguments: # BASISOBJ ... A basis function object object. # CVEC ... coefficient vector defining density, of length NBASIS # CVAL ... normalizing constant defining density # MU ... mean value to be subtracted from variates # SIGMA ... standard deviation to define u = (x - mu)/sigma # RNG ... vector of length 2 giving the interval over which the # integration is to take place # NDERIV ... order of derivative required for basis function expectation # UWRD ... if (TRUE, expectation is of (D PHI)*U # JMAX ... maximum number of allowable iterations # EPS ... convergence criterion for relative stop # Return: # A vector SS of length NBASIS of integrals of functions. # 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) rangeval <- 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 the history of discrete approximations to the integral smat <- matrix(0,JMAXP,nbasis) 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) 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[1,] <- width*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) 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*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, rng=rangeval, JMAX=15, EPS=1e-7) { # Computes expectations of cross product of basis functions with # respect to density # p(x) = Cval^{-1} exp t(c) %*% phi(x) # by numerical integration using Romberg integration # Arguments: # BASISOBJ ... A basis function object. # CVEC ... coefficient vector defining density # CVAL ... normalizing constant defining density # RNG ... vector of length 2 giving the interval over which the # integration is to take place # JMAX ... maximum number of allowable iterations # EPS ... convergence criterion for relative stop # Return: # A matrix of order NBASIS of integrals of functions. # 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) rangeval <- 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 the history of discrete approximations to the integral smat <- array(0,c(JMAXP,nbasis,nbasis)) # the first iteration uses just the }points x <- rng nx <- length(x) fx <- getbasismatrix(x, basisobj) 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, 1) } oneb <- matrix(1,1,nbasis) sumj <- t(Dfx1) %*% ((px %*% oneb) * Dfx2) smat[1,,] <- width*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) 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, 1) } sumj <- t(Dfx1) %*% ((px %*% oneb) * Dfx2) smat[j,,] <- (smat[j-1,,] + width*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) yadim <- dim(ya) if (is.null(yadim)) { yadim <- n nydim <- 1 } else { nydim <- length(yadim) } 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) ) }