wenjie2wang / splines2

Regression Spline Functions and Classes
https://wwenjie.org/splines2
GNU General Public License v3.0
42 stars 3 forks source link

B/M/C/I-splines: different results for basis in R and Cpp #21

Closed ccardehu closed 7 months ago

ccardehu commented 7 months ago

Hello,

I was wondering if anyone has found differences between the R and Rcpp basis outputs of the XX-splines family of functions. Below is a (reproducible) toy example of an Rcpp implementation of the Isplines function and the corresponding R function. As you can see, the output is different. The differences are smaller as the input vector "x" grows, but still.

Many thanks for your time and for any comment or pointer.

Best, Camilo

Rcpp code (saved as test.cpp):

#include <RcppArmadillo.h>
#include <splines2Armadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(splines2)]]

// [[Rcpp::export]]
arma::mat testCpp(arma::vec& x, arma::vec& knots, const unsigned int deg){
  const arma::vec bouK = {0.0,1.0};
  splines2::ISpline iSp(x,knots,deg,bouK);
  arma::mat out = iSp.basis(false);
  return(out);
}

R Code:

set.seed(1234)
X = runif(10)
knots = c(0.3,0.5,0.7)
degree = 2
resR = splines2::isp(X, knots = knots, degree = degree, intercept = F,Boundaryknots. = c(0,1))
Rcpp::sourceCpp("test.cpp")
resCpp = testCpp(X,knots,degree)

Output:

head(resR,5):
             1          2         3          4 5
[1,] 0.1850983 0.01150104 0.0000000 0.00000000 0
[2,] 1.0000000 0.99150787 0.6193809 0.07021538 0
[3,] 1.0000000 0.98648150 0.5688592 0.05008625 0
[4,] 1.0000000 0.99185709 0.6235199 0.07209209 0
[5,] 1.0000000 1.00000000 1.0000000 1.00000000 1

head(resRcpp,5):
          [,1]      [,2]      [,3]       [,4]      [,5]
[1,] 0.2063023 0.0140001 0.0000000 0.00000000 0.0000000
[2,] 1.0000000 0.9916231 0.5063851 0.03658500 0.0000000
[3,] 1.0000000 0.9866649 0.4630147 0.02609692 0.0000000
[4,] 1.0000000 0.9919676 0.5099709 0.03756284 0.0000000
[5,] 1.0000000 1.0000000 0.9743759 0.70877256 0.1543224

Session info:

R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22621)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8    LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
[5] LC_TIME=Spanish_Colombia.utf8    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_4.2.2           tools_4.2.2              RcppArmadillo_0.12.8.0.0 Rcpp_1.0.12              splines2_0.5.1  
wenjie2wang commented 7 months ago

There is a typo in your example:

set.seed(1234)
X = runif(10)
knots = c(0.3,0.5,0.7)
degree = 2
## The `Boundary.knots` was not correctly specified in the following comment.
## resR = splines2::isp(X, knots = knots, degree = degree, intercept = F,Boundaryknots. = c(0,1))
resR = splines2::isp(X, knots = knots, degree = degree, intercept = F, Boundary.knots = c(0,1))
Rcpp::sourceCpp("test.cpp")
resCpp = testCpp(X,knots,degree)

You should be able to get the same results after fixing the typo.

ccardehu commented 7 months ago

There is a typo in your example:

set.seed(1234)
X = runif(10)
knots = c(0.3,0.5,0.7)
degree = 2
## The `Boundary.knots` was not correctly specified in the following comment.
## resR = splines2::isp(X, knots = knots, degree = degree, intercept = F,Boundaryknots. = c(0,1))
resR = splines2::isp(X, knots = knots, degree = degree, intercept = F, Boundary.knots = c(0,1))
Rcpp::sourceCpp("test.cpp")
resCpp = testCpp(X,knots,degree)

You should be able to get the same results after fixing the typo.

Omg, I'm sorry! Thank you very much for your time (silly question!)