Many thanks to Prof. Douglas Bates for a fast and detailed solution to
my problem of generating a B-spline basis for estimation of a periodic
function.
I have included the response from Prof. Bates below. Also, I have
included a couple of functions I wrote to generate the basis. The
functions have been tested using spline.des() from the splines package
on StatLib, but I think that spline.des() in S+3.4 (HP-UX) is identical.
Note that the functions I wrote are experimental; if any of you find
any errors or implement improvements I am very interested (one
improvement would be the ability to exclude the intercept from the
basis).
Regards,
Henrik
--------------------------------------------------------------------------
Henrik Aalborg Nielsen
Department of Mathematical Modelling
Section for Statistics Fax: +45 4588 1397
Time Series Group Phone: +45 4525 3418
Technical University of Denmark E-mail: han@imm.dtu.dk
Building 321, 2800 Lyngby, Denmark URL: http://www.imm.dtu.dk/~han
--------------------------------------------------------------------------
Douglas Bates <bates@stat.wisc.edu> writes:
> Henrik Aalborg Nielsen <han@imm.dtu.dk> writes:
>
> > Hi Douglas,
> >
> > Thank you for your very fast reply!
> >
> > I am actually looking at your "splines" package and the function
> > "periodic.spline", but I cannot figure how to select the knots and modify
> > the output from "spline.des", when I don't want a knot at every data
> > point.
>
> spline.des allows you to specify the knot pattern that you want to
> use. The thing that periodic.spline does is to re-arrange the entries
> in that matrix to enforce the periodicity.
>
> > To be more specific: My data are recorded at every full hour of the day
> > (so x allways one of 0,1,...,23) and I want to model a diurnal variation
> >using maybe 6 degrees of freedom.
>
> Suppose then that you make the knots every four hours and you have a
> fourth order spline (order is degree + 1 so a cubic spline is fourth
> order). Your knots over the range of the data will be at 0, 4, 8, 12,
> 16, and 20. To fit a cubic spline (i.e. order = 4 since order =
> degree + 1) you need three knots to the left and three to the right of
> the data. Your total knot vector will be -12 -8 -4 0 4 8 12 16 20 24
> 28 32 36. To fit the coefficients for the periodic spline you add the
> last three columns onto the first three.
>
> I did a quick version using the splines package in R (similar to the
> S-PLUS version but with slight changes in the names in use)
>
> R> tm <- trunc(runif(40, min = 0, max = 24))
> R> tm # sample time data
> [1] 13 14 15 14 9 18 1 5 17 1 9 13 0 3 9 4 15 17 22 0
> [21] 2 1 22 23 5 13 9 10 7 1 14 17 13 2 9 11 2 12 9 22
> R> table(tm)
> 0 1 2 3 4 5 7 9 10 11 12 13 14 15 17 18 22 23
> 2 4 3 1 1 2 1 6 1 1 1 4 3 2 3 1 3 1
> R> seq(from = -12, to = 36, by = 4)
> [1] -12 -8 -4 0 4 8 12 16 20 24 28 32 36
> R> dd <- splineDesign(knots = seq(from = -12, to = 36, by = 4), x = tm)
> R> dim(dd)
> [1] 40 9
> R> dd[, 1:3] <- dd[, 1:3] + dd[, 7:9]
> R> dd <- dd[, 1:6]
> R> options(digits = 5)
> R> dd
> [,1] [,2] [,3] [,4] [,5] [,6]
> [1,] 0.0026042 0.0000000 0.0000000 0.0703125 0.6119792 0.3151042
> [2,] 0.0208333 0.0000000 0.0000000 0.0208333 0.4791667 0.4791667
> [3,] 0.0703125 0.0000000 0.0000000 0.0026042 0.3151042 0.6119792
> [4,] 0.0208333 0.0000000 0.0000000 0.0208333 0.4791667 0.4791667
> [5,] 0.0000000 0.0000000 0.0703125 0.6119792 0.3151042 0.0026042
> [6,] 0.4791667 0.0208333 0.0000000 0.0000000 0.0208333 0.4791667
> [7,] 0.0703125 0.6119792 0.3151042 0.0026042 0.0000000 0.0000000
> [8,] 0.0000000 0.0703125 0.6119792 0.3151042 0.0026042 0.0000000
> [9,] 0.3151042 0.0026042 0.0000000 0.0000000 0.0703125 0.6119792
> [10,] 0.0703125 0.6119792 0.3151042 0.0026042 0.0000000 0.0000000
> [11,] 0.0000000 0.0000000 0.0703125 0.6119792 0.3151042 0.0026042
> [12,] 0.0026042 0.0000000 0.0000000 0.0703125 0.6119792 0.3151042
> [13,] 0.1666667 0.6666667 0.1666667 0.0000000 0.0000000 0.0000000
> [14,] 0.0026042 0.3151042 0.6119792 0.0703125 0.0000000 0.0000000
> [15,] 0.0000000 0.0000000 0.0703125 0.6119792 0.3151042 0.0026042
> [16,] 0.0000000 0.1666667 0.6666667 0.1666667 0.0000000 0.0000000
> [17,] 0.0703125 0.0000000 0.0000000 0.0026042 0.3151042 0.6119792
> [18,] 0.3151042 0.0026042 0.0000000 0.0000000 0.0703125 0.6119792
> [19,] 0.4791667 0.4791667 0.0208333 0.0000000 0.0000000 0.0208333
> [20,] 0.1666667 0.6666667 0.1666667 0.0000000 0.0000000 0.0000000
> [21,] 0.0208333 0.4791667 0.4791667 0.0208333 0.0000000 0.0000000
> [22,] 0.0703125 0.6119792 0.3151042 0.0026042 0.0000000 0.0000000
> [23,] 0.4791667 0.4791667 0.0208333 0.0000000 0.0000000 0.0208333
> [24,] 0.3151042 0.6119792 0.0703125 0.0000000 0.0000000 0.0026042
> [25,] 0.0000000 0.0703125 0.6119792 0.3151042 0.0026042 0.0000000
> [26,] 0.0026042 0.0000000 0.0000000 0.0703125 0.6119792 0.3151042
> [27,] 0.0000000 0.0000000 0.0703125 0.6119792 0.3151042 0.0026042
> [28,] 0.0000000 0.0000000 0.0208333 0.4791667 0.4791667 0.0208333
> [29,] 0.0000000 0.0026042 0.3151042 0.6119792 0.0703125 0.0000000
> [30,] 0.0703125 0.6119792 0.3151042 0.0026042 0.0000000 0.0000000
> [31,] 0.0208333 0.0000000 0.0000000 0.0208333 0.4791667 0.4791667
> [32,] 0.3151042 0.0026042 0.0000000 0.0000000 0.0703125 0.6119792
> [33,] 0.0026042 0.0000000 0.0000000 0.0703125 0.6119792 0.3151042
> [34,] 0.0208333 0.4791667 0.4791667 0.0208333 0.0000000 0.0000000
> [35,] 0.0000000 0.0000000 0.0703125 0.6119792 0.3151042 0.0026042
> [36,] 0.0000000 0.0000000 0.0026042 0.3151042 0.6119792 0.0703125
> [37,] 0.0208333 0.4791667 0.4791667 0.0208333 0.0000000 0.0000000
> [38,] 0.0000000 0.0000000 0.0000000 0.1666667 0.6666667 0.1666667
> [39,] 0.0000000 0.0000000 0.0703125 0.6119792 0.3151042 0.0026042
> [40,] 0.4791667 0.4791667 0.0208333 0.0000000 0.0000000 0.0208333
bs.per.ek <- function(x, period=24, degree=3, nknots=6)
{
## EXPERIMENTAL
## The strange name means: B-spline, periodic, equidistant knots
## x must lie in [0 ; period) (if x == period it should be recorded as 0).
if(any(x < 0 | x >= period)) stop("x out of bounds")
if(degree < 1) stop("degree must a positive integer")
## Construct knots
knots.x <- seq(0, period, length=nknots+1)
delta <- diff(knots.x[1:2])
knots.left <- seq(from = -delta, by = -delta, length = degree)
knots.right <- seq(from = period + delta, by = delta, length = degree)
knots <- c(knots.left, knots.x, knots.right)
## Design matrix:
X <- spline.des(x=x, knots=knots, ord=degree+1)$design
X[, 1:degree] <- X[, 1:degree] + X[, ncol(X) - (degree-1):0]
X <- X[, 1:(ncol(X)-degree)]
return(X)
}
bs.per <- function(x, period=24, degree=3, knots=6)
{
## EXPERIMENTAL
## The name means: B-spline basis, periodic
## x must lie in [0 ; period) (if x == period it should be recorded as 0).
## Especially the non-equidistant placement of knots is somewhat experimental.
## If the konts argument has length one it is intrepreted as
## the number of equidistant knots:
if(length(knots) == 1)
return(bs.per.ek(x=x, period=period, degree=degree, nknots=knots))
## Some checks:
if(degree < 1) stop("degree must a positive integer")
if(any(x < 0 | x >= period)) stop("x out of bounds")
if(any(knots < 0 | knots >= period)) stop("knots out of bounds")
if(any(diff(knots) <= 0)) stop("knots must be strictly increasing")
if(length(knots) < degree + 1)
stop("don't know what to do when the No. of knots is less than degree+1")
## Construct knots
if(knots[1] > 0) {
knots <- knots - knots[1]
x <- x - knots[1]
x[x < 0] <- x[x < 0] + period
}
knots.x <- knots
knots.left <- knots[length(knots) - (degree-1):0] - period # last 'degree'
knots
knots.right <- knots[1:(degree+1)] + period # first 'degree + 1' knots
knots.all <- c(knots.left, knots.x, knots.right)
## Design matrix:
X <- spline.des(x=x, knots=knots.all, ord=degree+1)$design
X[, 1:degree] <- X[, 1:degree] + X[, ncol(X) - (degree-1):0]
X <- X[, 1:(ncol(X)-degree)]
return(X)
}
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
|