The generated basis slightly differs from one created in R

[martinmodrak]

Hi, your code has been super helpful in getting me started for a more complex model, huge thanks. I however noticed a minor issue that kept me scratching my head for a while: if I fit coefficients with your Stan code and then use these coefficients with basis generated in R with a call to bs, the results slightly differ.

The culprit seems to be the last line of the transformed data block in both b_spline.stan and b_spline_penalized.stan which reads

B[num_knots + spline_degree - 1, num_data] = 1; 

This seems to assume that the last element of X is exactly the last knot, which does not work as intended when X is not ordered (there is a weird kink in the resulting function). And even when X is ordered, the resulting basis differs in this position from the one generated in R.

When the last element of X is exactly the last knot, this returns the same basis as the R function, but that seems to be more of a quirk of the bs function - in all definitions I found, B-spline should be zero/undefined at the last knot: it is defined on left-closed, right-open interval [knot_0, knot_n). Further, I have not found any mention of such adjustment in any description of the recursive definition of B-splines.

When I remove this line, the generated basis is exactly the same as the basis generated in R for all points within the knot interval and all my inferences work smoothly.

[martinmodrak]

To make it work for all cases, the last line of transformed data needs to be replaced with

for(i in 1:num_data) {
    if(X[i] == ext_knots[num_knots + 2*spline_degree - 1]) {
      B[num_basis, i] = 1;
    }
}

this probably does not need to be put in your original code, but putting it here, in hopes that someone who wants to adapt the code will find it helpful.