Functions of splines2::mSpline do not sum to for periodic = TRUE and non-equidistant knots

[BerriJ]

Hi @wenjie2wang,

Sorry to bother you again concerning splines2::mSpline. I noticed that the spline functions do not sum to 1 if periodic = TRUE and non-equidistant knots are used. Consider the following reprex:

# install.packages("pbs")
# install.packages("splines2")

x <- 0:1000 / 1000
order <- 2
deg <- order - 1

knots <- c(0.00, 0.082, 0.23, 0.47, 1.00)

par(mfrow = c(1,2))

ts.plot(
    pbs::pbs(x,
        # pbs add boundary knots automatically
        knots = knots[c(-1, -5)],
        degree = deg, intercept = TRUE
    ),
    col = seq_along(knots),
    lwd = 2
)

ts.plot(
    splines2::mSpline(x,
        knots = knots[c(-1, -5)],
        degree = deg,
        Boundary.knots = c(0, 1),
        periodic = TRUE,
        intercept = TRUE
    ),
    col = seq_along(knots),
    lwd = 2
)

Do you know if this is the supposed behavior? For my application, I need the functions to sum to one.

Thank you very much. Best, BerriJ

[wenjie2wang]

Do you know if this is the supposed behavior?

Yes, it is expected behavior for M-spline basis functions regardless of the argument periodic, which is the essential difference between B-spline and M-spline basis functions. The M-spline basis functions are normalized so that their integrals are one. Consider the following example:

library(splines2)

x <- 0:1000 / 1000
order <- 2
deg <- order - 1
knots <- c(0.00, 0.082, 0.23, 0.47, 1.00)

res <- mSpline(
    x,
    knots = knots[c(-1, -5)],
    degree = deg,
    Boundary.knots = c(0, 1),
    periodic = TRUE,
    intercept = TRUE,
    integral = TRUE # check the integral
)
matplot(x, res, type = "l")

^{Created on 2023-03-27 with reprex v2.0.2}

[BerriJ]

Thank you for this insight. So, unfortunately, it's not possible to create periodic B-Splines using this package, right? Would you consider adding this functionality?

[wenjie2wang]

No problem. I would refer you to Section 2.2 of 10.6339/21-JDS1020 and see if you can apply (1) to obtain the periodic B-splines from the M-spline counterparts. I will consider adding it to the package.

[BerriJ]

Thank you very much. I was unaware of (1), and it works fine (although one has to tweak the simple knot sequence a bit to reflect the periodicity). I'm busy now, but I'll add some code to this issue next week.

If you like, I can also add a PR to integrate this into splines2. Ideas I came up with are:

• Adding a periodic argument to the bspline function
• As a method for the m-spline object for example, to_bspline()

[wenjie2wang]

Glad to hear that (1) works!

Thank you for expressing interest in submitting a PR! I have very similar ideas regarding this feature request. We can construct periodic B-splines based on M-splines inside bSpline(). However, it is just a quick workaround and can be inefficient. In fact, it is straightforward to implement a PeriodicBSpline class in C++ following the existing PeriodicMSpline class. I think it better to implement a template class for them. A PR is welcome if you feel comfortable to implement the template class in C++.

[BerriJ]

So, unfortunately, I have barely touched template classes yet, so implementing one for the spline functions is out of scope for me now. However, I made a wrapper for the m-spline class so that it applies (1). Maybe this saves you from some indexing troubles when creating the weight vector.

// periodic_bsplines.cpp

#include <RcppArmadillo.h>
// include header file from splines2 package
#include <splines2Armadillo.h>

// [[Rcpp::export]]
arma::mat splines2_periodic(const arma::vec &x,
                            const arma::vec &knots,
                            const unsigned int deg,
                            const bool &intercept = true)
{

    unsigned int order = deg + 1;
    arma::uvec inner_idx = arma::regspace<arma::uvec>(order,
                                                      knots.n_elem - order - 1);
    arma::uvec bound_idx = {deg, knots.n_elem - order};

    // Create periodic mspline object
    splines2::PeriodicMSpline ps(x, knots(inner_idx), deg, knots(bound_idx));
    arma::mat ps_mat = ps.basis(true);

    if (!intercept)
        ps_mat.shed_col(0);

    // We will use https://doi.org/10.6339/21-JDS1020 eq. (1) to convert
    // the mspline basis to a bspline basis

    // We need this sequence to calculate the weights
    arma::vec knots_ext = knots.subvec(bound_idx(0), bound_idx(1));
    knots_ext = join_cols(knots_ext,
                          knots(inner_idx.head(deg)) + knots(bound_idx(1)));

    for (unsigned int i = 0; i < ps_mat.n_cols; i++)
    {
        double w = knots_ext(1 - intercept + i + order) -
                   knots_ext(1 - intercept + i);
        ps_mat.col(i) *= w / order;
    }

    return ps_mat;
}

/*** R

x <- 0:100 / 100
deg <- 2

knots <- c(-0.2, -0.1, 0.0, 0.1, 0.2, 0.5, 1.0, 1.5, 2.1)

res <- splines2_periodic(x, knots, deg, TRUE)
ts.plot(res)
rowSums(res)
*/

The above can be executed using Rcpp::sourceCpp("periodic_bsplines.cpp"). I also uploaded the above as gist.

[wenjie2wang]

Thanks for the reference code! I will let you know after I implement the PeriodicBSpline, which should be straightforward enough without the need of scaling.

[wenjie2wang]

I added periodicBSpline mainly in this commit for your reference: https://github.com/wenjie2wang/splines2/commit/f27b4ed715957bc62615f96ada6c6676c8172bc6.