TABLE OF CONTENTS

splineUtils/splineconv [ Functions ]

NAME

    splineconv --- compute the convolution of two splines

FUNCTION

This is needed to construct a penalty on the integrated squared second derivative. This function computes the convolution of two spline basis functions, that is, 

       int_0^infty ( x^k B_1(x) * B_2(x) ) dx

where the splines may be of different orders, but are defined on the same set of knots.

SYNOPSIS

1572 splineconv <- function(k, n1, j1, n2, j2, knots)

INPUTS

    k      the power of x used in the convolution
    n1     order of the first spline
    j1     largest knot of the first spline
    n2     order of the second spline
    j2     largest knot of the second spline
    knots  set of knots on which the splines are defined

OUTPUTS

    The k-th order convolution of the splines defined by the input parameters

SOURCE

1575 {
1576     # if the splines don't overlap, the convolution is 0
1577     if(j1 - n1 >= j2 | j2 - n2 >= j1) return(0)
1578     # if both splines are first-order, the convolution is trivial
1579     if(n1 == 1 & n2 == 1){
1580         out <- 1 / (k + 1) * (knots[ki(knots, j1)]^(k + 1) - knots[ki(knots, j1 - 1)]^(k + 1))
1581         return(out)
1582     }
1583     # By symmetry, we can assume that n1>n2 wlog. If this is not the case, switch them around
1584     if(n2 > n1){
1585         n3 <- n1; n1 <- n2; n2 <- n3
1586         j3 <- j1; j1 <- j2; j2 <- j3
1587     }
1588     # use a magic formula that can be derived by integration by parts
1589     out <- 0
1590     denom1 <- knots[ki(knots, j1 - 1)] - knots[ki(knots, j1 - n1)]
1591     denom2 <- knots[ki(knots, j1)] - knots[ki(knots, j1 - n1 + 1)]
1592     if(denom1 > 0){
1593         out <- out + 1 / denom1 * splineconv(k + 1, n1 - 1, j1 - 1, n2, j2, knots)
1594         out <- out - knots[ki(knots, j1 - n1)] / denom1 * 
1595             splineconv(k, n1 - 1, j1 - 1, n2, j2, knots)
1596     }
1597     if(denom2 > 0){
1598         out <- out + knots[ki(knots, j1)] / denom2 * 
1599             splineconv(k, n1 - 1, j1, n2, j2, knots)
1600         out <- out - 1 / denom2 * splineconv(k + 1, n1 - 1, j1, n2, j2, knots)
1601     }
1602     return(out)
1603 }