R/monomialpen.R

In drtagkim/mcgillfdar: Functional Data Analysis

monomialpen <- function(basisobj, Lfdobj=int2Lfd(2),
                        rng=basisobj$rangeval)
{
#  MONOMIALPEN  Computes a monomial basis penalty matrix.
#
#  Warning:  This version is incomplete in the sense that
#    it only works with LFDOBJ = D^m
#
#  Arguments:
#  BASISOBJ ... a monomial basis object
#  LFDOBJ   ... either the order of derivative or a
#               linear differential operator to be penalized.
#  RNG      ... a range over which the penalty matrix is computed.
#  Returns a list the first element of which is the basis matrix
#   and the second element of which is the diagonal of the penalty matrix.

#  Last modified:  26 October 2005

#  check BASISOBJ

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

type <- basisobj$type
if (!(type == "monom")) stop("BASISOBJ not of type monom")

#  check LFDOBJ

Lfdobj <- int2Lfd(Lfdobj)

#  check whether LFDOBJ is of form D^m

nderiv <- Lfdobj$nderiv

#  get basis information

nbasis    <- basisobj$nbasis
exponents <- basisobj$params

#  check exponents

for (ibasis in 1:nbasis) {
	ideg <- exponents[ibasis]
	if (ideg-floor(ideg) != 0) stop(
		"An exponent is not an integer.")
}

#  compute the penalty matrix

penaltymat <- matrix(0,nbasis,nbasis)
for (ibasis in 1:nbasis) {
    ideg <- exponents[ibasis]
    ifac <- 1
    if (nderiv > 0) {
	     ifac <- ideg
	     if (nderiv > 1) for (k in 2:nderiv) ifac <- ifac*(ideg - k + 1)
    }
    for (jbasis in 1:ibasis) {
        jdeg <- exponents[jbasis]
        jfac <- 1
        if (nderiv > 0) {
	         jfac <- jdeg
            if (nderiv > 1) for (k in 2:nderiv) jfac <- jfac*(jdeg - k + 1)
        }
        if ((ideg >= nderiv) & (jdeg >= nderiv)) {
	        ipow <- ideg+jdeg-2*nderiv+1
            penaltymat[ibasis,jbasis] <-
                (rng[2]^ipow - rng[1]^ipow)*ifac*jfac/ipow
            penaltymat[jbasis,ibasis] <- penaltymat[ibasis,jbasis]
        }
    }
}
penaltymat
}