File: mSpline.R

package info (click to toggle)
r-cran-splines2 0.4.1-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 1,512 kB
  • sloc: cpp: 1,988; ansic: 165; sh: 13; makefile: 2
file content (166 lines) | stat: -rw-r--r-- 6,859 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
##
## R package splines2 by Wenjie Wang and Jun Yan
## Copyright (C) 2016-2021
##
## This file is part of the R package splines2.
##
## The R package splines2 is free software: You can redistribute it and/or
## modify it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or any later
## version (at your option). See the GNU General Public License at
## <https://www.gnu.org/licenses/> for details.
##
## The R package splines2 is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
##

##' M-Spline Basis for Polynomial Splines
##'
##' Generates the basis matrix of regular M-spline, periodic M-spline, and the
##' corresponding integrals and derivatives.
##'
##' This function contains an implementation of the close form M-spline basis
##' based on the recursion formula given by Ramsay (1988).  For monotone
##' regression, one can use I-splines (see \code{\link{iSpline}}) instead of
##' M-splines.
##'
##' @inheritParams bSpline
##'
##' @param df Degree of freedom that equals to the column number of the returned
##'     matrix.  One can specify \code{df} rather than \code{knots}.  For
##'     M-splines, the function chooses \code{df - degree -
##'     as.integer(intercept)} internal knots at suitable quantiles of \code{x}
##'     ignoring missing values and those \code{x} outside of the boundary.  For
##'     periodic spline based on M-spline (\code{periodic = TRUE}), \code{df -
##'     as.integer(intercept)} internal knots will be chosen instead and the
##'     number of internal knots must be greater or equal to the specified
##'     \code{degree}.  If internal knots are specified via \code{knots}, the
##'     specified \code{df} will be ignored.
##' @param knots The internal breakpoints that define the splines.  The default
##'     is \code{NULL}, which results in a basis for ordinary polynomial
##'     regression.  Typical values are the mean or median for one knot,
##'     quantiles for more knots.  For periodic splines (\code{periodic =
##'     TRUE}), the number of knots must be greater or equal to the specified
##'     \code{degree}.
##' @param Boundary.knots Boundary points at which to anchor the splines.  By
##'     default, they are the range of \code{x} excluding \code{NA}.  If both
##'     \code{knots} and \code{Boundary.knots} are supplied, the basis
##'     parameters do not depend on \code{x}. Data can extend beyond
##'     \code{Boundary.knots}.  For periodic splines (\code{periodic = TRUE}),
##'     the specified \code{Boundary.knots} define the cyclic period.
##' @param periodic A logical value.  If \code{TRUE}, the periodic splines will
##'     be returned instead of regular M-splines.  The default value is
##'     \code{FALSE}.
##' @param derivs A non-negative integer specifying the order of derivatives of
##'     M-splines. The default value is \code{0L} for M-spline bases.
##' @param integral A logical value.  If \code{TRUE}, the corresponding
##'     integrals of spline bases will be returned.  The default value is
##'     \code{FALSE}.  For periodic spline, the integral of each basis is
##'     integrated from the left boundary knot.
##'
##' @return A numeric matrix of \code{length(x)} rows and \code{df} columns if
##'     \code{df} is specified.  If \code{knots} are specified instead, the
##'     output matrix will consist of \code{length(knots) + degree +
##'     as.integer(intercept)} columns if \code{periodic = FALSE}, or
##'     \code{length(knots) + as.integer(intercept)} columns if \code{periodic =
##'     TRUE}.  Attributes that correspond to the arguments specified are
##'     returned for usage of other functions in this package.
##'
##' @references
##' Ramsay, J. O. (1988). Monotone regression splines in action.
##' \emph{Statistical science}, 3(4), 425--441.
##'
##' @example inst/examples/ex-mSpline.R
##'
##' @seealso
##' \code{\link{bSpline}} for B-splines;
##' \code{\link{iSpline}} for I-splines;
##' \code{\link{cSpline}} for C-splines.
##'
##' @export
mSpline <- function(x, df = NULL, knots = NULL, degree = 3L,
                    intercept = FALSE, Boundary.knots = NULL,
                    periodic = FALSE, derivs = 0L, integral = FALSE,
                    ...)
{
    ## check inputs
    if ((derivs <- as.integer(derivs)) < 0) {
        stop("The 'derivs' must be a non-negative integer.")
    }
    if ((degree <- as.integer(degree)) < 0)
        stop("The 'degree' must be a nonnegative integer.")
    if (is.null(df)) {
        df <- 0L
    } else {
        df <- as.integer(df)
        if (df < 0) {
            stop("The 'df' must be a nonnegative integer.")
        } else if (periodic && df < degree) {
            stop("The 'df' must be >= 'degree' for periodic spline basis.")
        }
    }
    knots <- null2num0(knots)
    Boundary.knots <- null2num0(Boundary.knots)
    ## take care of possible NA's in `x`
    nax <- is.na(x)
    if (all(nax)) {
        stop("The 'x' cannot be all NA's!")
    }
    ## remove NA's in x
    xx <- if (nas <- any(nax)) {
              x[! nax]
          } else {
              x
          }
    ## call the engine function
    out <- if (periodic) {
               rcpp_periodic_mSpline(
                   x = xx,
                   df = df,
                   degree = degree,
                   internal_knots = knots,
                   boundary_knots = Boundary.knots,
                   derivs = derivs,
                   integral = integral,
                   complete_basis = intercept
               )
           } else {
               rcpp_mSpline(
                   x = xx,
                   df = df,
                   degree = degree,
                   internal_knots = knots,
                   boundary_knots = Boundary.knots,
                   derivs = derivs,
                   integral = integral,
                   complete_basis = intercept
               )
           }
    ## throw warning if any x is outside of the boundary
    b_knots <- attr(out, "Boundary.knots")
    if (! periodic && any((xx < b_knots[1L]) | (xx > b_knots[2L]))) {
        warning(wrapMessages(
            "Some 'x' values beyond boundary knots",
            "may cause ill-conditioned bases."
        ))
    }
    ## keep NA's as is
    if (nas) {
        nmat <- matrix(NA, length(nax), ncol(out))
        nmat[! nax, ] <- out
        saved_attr <- attributes(out)
        saved_attr$dim[1] <- length(nax)
        out <- nmat
        attributes(out) <- saved_attr
        attr(out, "x") <- x
    }
    ## add dimnames for consistency
    name_x <- names(x)
    if (! is.null(name_x)) {
        row.names(out) <- name_x
    }
    ## add class
    class(out) <- c("matrix", "mSpline")
    out
}