Get the spline design matrix in Rcpp

[mengqi00]

I am trying to use the 'splines2' package in Rcpp, but I am confused about how to get the same design matrix as the splines::spineDesign() in r? I can only find the basis matrix. Are there functions that directly calculate the design matrix using x, knot sequence and order? Thank you very much!

[wenjie2wang]

Thanks for posting the question here. A simple example is as follows:

// example.cpp

#include <RcppArmadillo.h>
#include <splines2Armadillo.h>
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(splines2)]]

// [[Rcpp::export]]
arma::mat splineDesign2(const arma::vec& knot_seq,
                        const arma::vec& x,
                        const unsigned int ord = 4,
                        const unsigned int derivs = 0)
{
    splines2::BSpline obj { x, ord - 1, knot_seq };
    if (derivs > 0) {
        return obj.derivative(derivs);
    }
    return obj.basis();
}

# example.R

Rcpp::sourceCpp("example.cpp")

set.seed(123)
x <- runif(100, 0, 10)
ord <- 4
knots <- c(1, 3, 7)
knot_seq <- c(rep(0, ord), knots, rep(10, ord))

## basis
res1 <- splines::splineDesign(knot_seq, x, ord)
res2 <- splineDesign2(knot_seq, x, ord)
all.equal(res1, res2) # expect to be TRUE

## derivatives
res1 <- splines::splineDesign(knot_seq, x, ord, 1)
res2 <- splineDesign2(knot_seq, x, ord, 1)
all.equal(res1, res2) # expect to be TRUE

[mengqi00]

This makes sense! Thank you very much!

[wenjie2wang]

You are welcome.