Is my Restricted (Natural) Cubic Spline equation wrong?

Question

I am trying to fit a restricted cubic spline (natural cubic spline) with 4 knots to toy data, attempting to follow Hastie, Tibshirani, Friedman 2nd ed. 5.2.1 p.144-146, Eqs 5.4 and 5.5. Data: Is basically a transposed ‘S’ shape. R-code is:

n <© 100
x <- (1:n)/n
true <- ((exp(1.2*x)+1.5*sin(7*x))-1)/3
noise <- rnorm(n, 0, 0.15)
y <- true + noise
plot(x,y)

I set knots at: {.2, .4, .6, .8} and am fitting using the non-linear NLS() function in R, but I can’t get the S-shape of the data no matter what I try.

My equations is wrong ? Or I am completely off-base in my approach? Any suggestions?

(Book-excerpt, my equation, and data-plot posted below)

A natural cubic spline with $K$ knots is represented by $K$ basis functions. One can start from a basis for cubic splines, and derive the reduced basis by imposing the boundary constraints. For example, starting from the truncated power series basis described in section 5.2, we arrive at where $$d_k=\frac{(X-\xi_k)^3_+ - (X-\xi_K)^3_+}{\xi_K - \xi_k} \enspace .$$ Each of these basis functions can be seend to have zero second and third derivative for $X\geq \xi_K$. $$y = \beta_0 + \beta_1 x + \beta_2\left(\left[\frac{(x-k_1)^3_+ - (x-k_4)^3_+}{k_4 - k_1}\right] - \left[\frac{(x-k_3)^3_+ - (x-k_4)^3_+}{k_4 - k_3}\right] \right) + \beta_3\left(\left[\frac{(x-k_2)^3_+ - (x-k_4)^3_+}{k_4 - k_2}\right] - \left[\frac{(x-k_3)^3_+ - (x-k_4)^3_+}{k_4 - k_3}\right] \right)$$

---

Can I ask a simpler question: Sites/books say: for my natural cubic spline approach (i.e. restricted cubic spline) w/ 4 knots, I need 4 basis functions. Is it Beta_0, Beta_1*x, and '4 more' ? Or is indeed just 4 betas (conceptually, as I have above) ? Thank you.

---

Thank you for your guidance. I am fitting the spline, along with many other modeling co-variates, i.e. explanatory variables. So, a simple canned package for the spline alone is not sufficient.
The actual shape I expect, when I fit to my modeling data, is a very positive skewed distributional form (I was trying to fit the fancy shaped example data, with an inflection, thinking it would be a good test of my coding of the spline functional-form).

I thought about snitching the functional form and calibrated-parameterization (from your Python above or from R) - but its a cubic-spline, not a natural cubic spline.

I get how my ftn is linear to the LHS of first knot. I guess next step is for me to see that various terms cancel, and indeed I'd be linear to the RHS of the right-most knot too. Also look for a 4-term basis that is said to be 'stable'.

Answer 1

Answering your second question first: for $K$ knots, natural cubic splines really have only $K$ degrees of freedom/basis functions.

This is quite easy to see if you count degrees of freedom: $K$ knots split the real line into $K+1$ intervals, and a piecewise cubic function in these intervals has $4(K+1)$ degrees of freedom. By matching function value, first and second derivative in each knot, we get $3K$ conditions, hence $K+4$ degrees of freedom are left. Finally, the "natural" spline condition says that we want the spline to be linear in the first and last interval, which kills another $2+2$ degrees of freedom (the quadratic and cubic coefficients in these intervals). All in all, we are left with $K$ degrees of freedom.

For your first question, it's really hard to tell what's going wrong without seeing your code. In general, I wouldn't recommend implementing spline routines yourself unless you have to; the basis you have used is very numerically unstable, and in general the code is a bit finicky to get right.

However, any serious scientific computing environment has these things built in, so generally you don't need to write it yourself. Here is your example in Python/scipy:

In [1]: %pylab
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

In [2]: n = 100

In [3]: x = linspace(0, 1, n)

In [10]: y = (((exp(1.2*x) + 1.5*sin(7*x))-1)/3) + normal(0, 0.15, size=n)

In [13]: from scipy.interpolate import LSQUnivariateSpline

In [15]: t = [0.2, 0.4, 0.6, 0.8]

In [16]: spl = LSQUnivariateSpline(x, y, t)

In [21]: plt.style.use('bmh')

In [22]: plot(x, spl(x), '-', x, y, 'o')

In [23]: plot(x, (((exp(1.2*x) + 1.5*sin(7*x))-1)/3), '--')

The blue line is the smoothing B-spline with four knots, the dashed line the original function without noise.

Answer 2

FYI, the formula from Hastie et al is correct. I just implemented it myself in R and compared the result to what I get from splines::ns() from R. The basis is not identical, but it spans the same space because the fitted values are exactly the same.

There is an error somewhere in your implementation. I followed the dumb approach, having implemented seven small functions for four knots: d1-d3, N1-N4, and it worked fine.

Answer 3

Here's an implementation of restricted cubic spline (based on this matlab code)

import scipy.linalg as lin

def rcs(x,y,knots):
    n = len(y)
    k = knots
    X1 = x
    q = len(k)-1
    myX=np.zeros((n,len(knots)-2))

    for j in range(q-1):
        tmp1 = (x-k[j])**3 * (x>k[j])
    tmp2 = (x-k[q-1])**3 * (x>k[q-1])*(k[q]-k[j])
    XX= tmp1-tmp2/(k[q]-k[q-1])
        tmp1 = (x-k[q])**3 * (x>k[q])
        tmp2 = (k[q-1]-k[j])
    XX = XX+tmp1*tmp2/(k[q]-k[q-1])
    myX[:,j]=XX

    X = np.hstack( (np.ones((n,1)),np.reshape(X1,(n,1)),myX) )
    bhat = np.linalg.lstsq(X,y)[0]
    bhatt = np.zeros(len(knots)+1)
    bhatt[len(bhat)] = (bhat[2:]*(k[0:-2]-k[-1])).sum()
    bhatt[len(bhat)] = bhatt[len(bhat)] / (k[-1]-k[-2])
    bhatt = np.hstack([bhatt, 0])    
    bhatt[-1] = (bhat[2:]*(k[0:-2]-k[-2])).sum()
    bhatt[-1] = bhatt[-1] / (k[-2]-k[-1])
    bhat = np.hstack((bhat, bhatt[-2:]))
    return bhat

def speval(x,coefs,knots):
    tmp = coefs[0] + coefs[1]*x
    for k in range(len(knots)): tmp = tmp + coefs[k+2]*((x-knots[k])**3)*(x>knots[k])
    return tmp



import pandas as pd
x = np.random.randn(300)*np.sqrt(2)
e = np.random.randn(300)*np.sqrt(0.5)
y = np.sin(x)+e
df = pd.DataFrame([x,y]).T
df.columns = ['x','y']
df = df.sort_index(by='x')
print df.head()
knots=np.array([-5.5938, -3.7732, -1.9526, -0.1320, 1.6886, 3.5092, 5.3298]);
bhat = rcs(df.x,df.y,knots)
print bhat
df['spline'] = speval(df.x, bhat, knots)
df2 = df.set_index('x')
df2[['y','spline']].plot()
plt.hold(True)
for k in knots: plt.plot(k,speval(k,bhat,knots),'rd')