R/polygpen.R

In fda: Functional Data Analysis

polygpen <- function(basisobj, Lfdobj=int2Lfd(1))
{

#  Computes the polygonal penalty matrix.
#  Arguments:
#  BASISOBJ ... a basis object of type "polyg"
#  LFDOBJ   ... either the order of derivative or a
#               linear differential operator to be penalized.
#          The highest derivative must be either 0 or 1.
#  Returns the penalty matrix.

#  Last modified 3 January 2008

if (!(inherits(basisobj, "basisfd"))) stop(
    "First argument is not a basis object.")

Lfdobj <- int2Lfd(Lfdobj)

type <- basisobj$type
if (type != "polyg" && type != "polygonal") 
stop("BASISOBJ not of type polyg or polygonal")

nderiv <- Lfdobj$nderiv

if (nderiv > 1) stop(
    "Derivative greater than 1 cannot be taken for polygonal basis.")

bwtlist <- Lfdobj$bwtlist
isintLfd <- TRUE
if (nderiv > 0) {
  for (ideriv in 1:nderiv) {
    fdj <- bwtlist[[ideriv]]
    if (!is.null(fdj)) {
      if (any(fdj$coefs != 0)) {
        isintLfd <- FALSE
        break
      }
    }
  }
}

if (isintLfd) {
    args    <- basisobj$params
    n       <- length(args)
    argdiff <- diff(args)
    penaltymat <- diag(rep(1,n))
    if (nderiv == 0) {
        penaltymat[1,1] = argdiff[  1]/3
        penaltymat[n,n] = argdiff[n-1]/3
        indx = 2:(n-1)
        for (i in indx)
            penaltymat[i,i] = (argdiff[i]+argdiff[i-1])/3
        indx = 2:n
        for (i in indx) {
            penaltymat[i,i-1] = argdiff[i-1]/6
            penaltymat[i-1,i] = argdiff[i-1]/6
        }
    } else {
        argdiff = 1/argdiff
        penaltymat[1,1] = argdiff[  1]
        penaltymat[n,n] = argdiff[n-1]
        indx = 2:(n-1)
        for (i in indx)
            penaltymat[i,i] = argdiff[i]+argdiff[i-1]
        indx = 2:n
        for (i in indx) {
            penaltymat[i,i-1] = -argdiff[i-1]
            penaltymat[i-1,i] = -argdiff[i-1]
        }
    }
    penaltymatrix <- penaltymat
} else {
    penaltymatrix <- inprod(basisobj, basisobj, Lfdobj, Lfdobj)
}

return( penaltymatrix )
}