TABLE OF CONTENTS

CmakeLikelihood/SmoothnessPenalty [ Functions ]

NAME

    SmoothnessPenalty --- compute smoothness penalties

FUNCTION

Compute the smoothness penalty for a spline CCurve as part of making likelihoods.

The penalty can be either a penalty on the sum of squared second differences or the integrated squared second derivative, or a Gaussian penalty. The penalty matrix and penalty type stored in the CCurve structure should be of matching type.

SYNOPSIS

975 double SmoothnessPenalty(curveP theCurve)

INPUTS

    theCurve  CCurve structure

OUTPUTS

    pen   the smoothness penalty, already divided by 2*(prior variance)
          that is, the entire smoothness penalty term.

SOURCE

 979 {
 980     double pen;
 981     int c1 = 1;
 982     // Gaussian penalty
 983     if(theCurve->SplinePenaltyType == pnone){ 
 984         pen = pow(F77_CALL(dnrm2)(&(theCurve->nj), theCurve->SplinePar, &c1),2.0);
 985     }else{
 986         double * temp = (double *) calloc(theCurve->nj, sizeof(double));
 987         double * par;
 988         // 2diff penalty works directly with parameters, 2deriv works with exp(theta)
 989         if(theCurve->SplinePenaltyType == pdiff) par = theCurve->SplinePar;
 990         else par = theCurve->SplineEPar;
 991         double c1d = 1.0;
 992         char uplo = 'U';
 993         // compute par %*% P %*% par
 994         F77_CALL(dsymv)( &uplo, &(theCurve->nj), &c1d, theCurve->SplinePenaltyMatrix,
 995                 &(theCurve->nj), par, &c1, &c1d, temp, &c1);
 996         pen = F77_CALL(ddot)( &(theCurve->nj), temp, &c1, par, &c1);
 997         if(theCurve->SplinePenaltyType == plogderiv) pen=log(pen+1);
 998         free(temp);
 999     }
1000     // normalize by prior variance
1001     pen = pen / (2 * theCurve->SplinePriorvar[0]);
1002     return pen;
1003 }