TABLE OF CONTENTS

makeLikelihood/smoothpen [ Functions ]

NAME

    smoothpen --- compute the smoothness penalty for a curve

FUNCTION

Penalize lack of smoothness of a B-spline curve as part of a likelihood computation

SYNOPSIS

1796 smoothpen <- function(curve, der = 0)

INPUTS

    curve      an RCurve structure
    der        the derivative of the smoothness penalty to compute

OUTPUTS

    value of the smoothness penalty

SOURCE

1799 {
1800     type <- curve$spline.penalty
1801     name <- curve$name
1802     theta <- curve$spline.par
1803     # extract the penalty matrix
1804     P <- curve$spline.penaltymatrix
1805     sigma2 <- curve$spline.priorvar
1806     if(der >= 2) stop("second derivative not implemented")
1807     # second difference penalty
1808     if(type == "2diff"){
1809         if(der == 0) return(max( t(theta)%*%P%*%theta / (2 * sigma2), 0))
1810         if(der == 1) return( P%*%theta /sigma2)
1811     }
1812     # second derivative penalty
1813     if(type == "2deriv"){
1814         et <- exp(theta)
1815         if(der == 0) return(max( t(et)%*%P%*%et / (2 * sigma2), 0))
1816         if(der == 1) return( mdiag(et)%*%P%*%et / sigma2 )
1817     }
1818     # penalty on the log 2nd derivative
1819     if(type == "log2deriv"){
1820         et <- exp(theta)
1821         ePe <- as.numeric(t(et)%*%P%*%et) 
1822         if(der == 0) return(max(log(ePe + 1)/ (2 * sigma2), 0))
1823         if(der == 1) return( mdiag(et)%*%P%*%et / sigma2 /(ePe + 1))
1824     }
1825     # gaussian "penalty", which isn't really a smoothness penalty at all.
1826     if(type == "none"){
1827         if(der == 0) return(theta%*%theta / (2 * sigma2))
1828         if(der == 1) return( theta / sigma2)
1829     }
1830 }