% This live script will cover some selected topics in splines.  Generally speaking, when trying to fit a dataset (x,y) under the assumption that y=f(x)+\epsilon, we can choose to represent the functionf\left(x\right) in many ways.  Many models, such a ordinary least squares, are additive models, where our function takes the form
% f(x)=\beta_0+\beta_1b_1(x)+\beta_2b_2(x)+\ldots+\beta_Kb_K(x)
% Here, the various \beta_i are coefficients we will try to learn from the data, while the functions b_i \left(x\right) are called basis functions. These basis functions can really be just about anything we want, from as simple as just linear functions (linear regression) to things like sine and cosine, polynomials, etc.  What functions we might want to use as basis functions for a given dataset will depend on the nature of the data, and any existing relationship we think might hold between x and y.
% One such choice of basis functions are splines.  In general, splines attempt to represent the function f\left(x\right) as a sum of several piecewise defined functions, typically polynomials.  When discussing splines (which again, are piecewise defined), we need to introduce the idea of knots: these are just locations x=c where the piecewise defined splines switch from one definition to another.  The choice of knots will be elaborated upon more below, but for now just note that if we use more knots, we will have more flexibility in our fit, since we will be essentially including more and more basis functions in our representation of f\left(x\right).
% As any calculus 1 student will know, a piecewise defined function can quite easily be discontinuous or non-differentiable at knots.  This is where regression splines come in.  The goal of regression splines is to choose a series of basis functions so that, given the knots you have chosen, the resulting f\left(x\right) is continuous and possibly also differentiable (up to maybe several orders).
% For example, suppose we wanted to use cubic polynomials for our basis functions, and a single knot located at some value in the interior of all the x values of our dataset.  Each cubic polynomial has 4 parameters, which are the coefficients of the various terms: so, the set of two such polynomials (left and right of the single knot) have 8 total parameters that can be adjusted to fit data.  But, if we further want to make it so that the resulting f\left(x\right) is continuous, that constraint will effectively remove one of these parameters, leaving 7.  We could further enforce differentiability, removing another parameter to leave 6.  Or we could enforce twice differentiability, to leave us with 5 parameters.  Of course, if we want to enforce thrice differentiability we would be left with only 4 parameters and a single cubic polynomial (the knot would do nothing), so we won't go that far.  Plus, a curve that is twice differentiable "looks" pretty nice and smooth, so we might be happy with that.  Then using cubic polynomials for our regression splines and enforcing twice differentiability with one knot has more parameters to fit to the data than a single cubic polynomial (plus we are free to choose where the knot should be).  If you imagine that we increase the number of knots from there, every additional knot will only add one more effective parameter, as the new region and associated cubic that is introduced each time has to be twice differentiable again, removing three of the four parameters that it brings with it.  In the end, then, with K knots, we end up with 4+K parameters that can be used to fit this function to the data.
% There are a variety of ways to write cubic splines that behave as indicated above.  One such representation is used in the code below, which defines basis functions h\left(x,c\right), where c is the location of a knot.  Below these basis functions are used to fit a 1-d dataset using least squares to find the coefficients for the basis functions.

h=@(x,c) (max(x-c,zeros(1,size(x,2)))).^3; %defining our basis function prototype.  c should be the location of a knot
%below we construct a dataset from some nonlinear function
xa=0; %this is the left hand side of our x domain
xb=9; %this is the right hand side of our x domain
n=50; % the number of data points
sig2=.5;
f=@(x) 4+15*x.*exp(-x); %this is just an example - try some other function if you want!
xvals=sort(xa+(xb-xa)*rand(n,1));
ytrue=f(xvals);
epsilon=normrnd(0,sqrt(sig2),n,1);
yvals=ytrue+epsilon;

K=10; %number of knots, which should be no more than n-4
cvals=quantile(xvals,K); %this just places the knots "evenly" over the range of xvals
xmat=[ones(n,1),xvals,xvals.^2,xvals.^3]; %these are the first 4 basis functions in this representation, which also correspond to the four terms in the first cubic polynomial
for knot=1:K
    xmat=[xmat,h(xvals,cvals(knot))]; %add the other basis functions, each of which is only nonzero to the right of a particular knot
end
beta=xmat\yvals; %ols using our basis functions
resid=yvals-xmat*beta; %residuals
mse=resid'*resid/n %MSE

%visualize
xplot=linspace(xa,xb,100)';
xplot=sort([xplot;xvals]);
yplot=beta(1)+beta(2)*xplot+beta(3)*xplot.^2+beta(4)*xplot.^3;
for knot=1:K
    yplot=yplot+beta(4+knot)*h(xplot,cvals(knot));
end
plot(xplot,yplot,'b-')
hold on
plot(xvals,yvals,'ro');

%do polynomial with same number of free coefficients for comparison.  Note
%this isn't great if K is too big.
degree=K+3;
betap=polyfit(xvals,yvals,degree);
plot(xplot,polyval(betap,xplot),'g-')
hold off

%
% If you run the code above several times, you will likely observe some instances where the spline fit seems to reach very large y values at the ends of the x values (or sometimes elsewhere).  This is an unfortunate side effect of the method.  The code below can be run in this case to make it easier to see what is happening in the rest of the region.

goodmin=min(yvals)-3;
goodmax=max(yvals)+3;
axis([xa,xb,goodmin,goodmax])

%
% However, there is a way to counteract this.  Instead of settling for the regression splines above, we can enforce two more conditions: demand that our function is linear in the region before the first knot and also linear in the region after the last knot.  This should largely remove this unwanted behavior, as a linear function is much less likely to grow to very large values in these regions.  Note that by adding these new constraints (really 4 constraints, since the linearity demands that the quadratic and cubic terms are both zero in both intervals), we end up with only K parameters to estimate.  This new set of basis functions are called natural (cubic) splines.  Again, there are multiple ways to represent these, but the code below uses one of them and does the same thing as above:

d=@(x,c,xK) (h(x,c)-h(x,xK))./(xK-c); % a building block of the natural spline basis functions.  xK is the location of the last knot.
basis=@(x,c,ind,K) (ind==1)*ones(1,size(x,2))+(ind==2)*x+(ind>2)*(d(x,c(max(ind-2,1)),c(K))-d(x,c(K-1),c(K))); %natural spline basis functions.  need to provide the index of the basis function you want from 1 to K

K=22; %number of knots, which should be no more than n
cvals=quantile(xvals,K); %this just places the knots "evenly" over the range of xvals
xmat=[basis(xvals,cvals,1,K)];
for baseind=2:K
    xmat=[xmat,basis(xvals,cvals,baseind,K)];
end
beta=xmat\yvals; %ols using our basis functions
resid=yvals-xmat*beta; %residuals
mse=resid'*resid/n %MSE

%visualize
xplot=linspace(xa,xb,100)';
yplot=beta(1);
for baseind=2:K
    yplot=yplot+beta(baseind)*basis(xplot,cvals,baseind,K);
end
hold on
plot(xplot,yplot)
hold off

% Upon examining the results of the code above, you should see that the curve is linear in the region before the first knot and after the last, as desired.  You should also note that, for a given value of K, the regression splines above will have a smaller (but maybe only slightly) MSE than the natural splines.  This is because the regression splines have more parameters available for a given K in order to make the data fit well to the function.
% Here is an interesting fact: given that natural cubic splines have K parameters to fit, and that these K parameters will directly effect the value of the overall curve at the K knots, if one chooses to use K=n, and place the knots at the location of the x values in the dataset (assuming these values are all unique), then the resulting curve will perfectly interpolate the dataset, leading to 0 (or some equivalent machine precision number) training MSE!  Try this above.
%
% Another, seemingly unrelated type of spline is called a smoothing spline. The way smoothing splines are defined is similar to the way we set up ridge regression and LASSO early in the semester.  That is, they are written as a minimization problem, where one term in the objective function that we are to minimize is the RSS (this is the fidelity or loss term, as its purpose is to make the function f similar to the data) and the other is some kind of penalty or regularizer term that tries to favor certain kinds of solutions over others, irrespective of their fidelity.  One such objective function is as follows:
% M=\sum_{i=1}^n(y_i-f(x_i))^2 + \lambda \int(f''(t))^2dt
% The problem is then to find the function f that minimizes the above objective - this is then your fit to the data.  If we examine this objective function, the RSS is the first term, as advertized: this term tries to make our function fit the data well by passing as close as possible to the datapoints in the RSS sense.  The second term in this case is asking us to penalize by some user-chosen non-negative parameter \lambda (this sets how strongly we care about the penalty term) multiplied by the integral (over the domain of the data, at least) of the second derivative of the function squared.  Note that this penatly term is inherently restricting us to consider only functions that are twice differentiable (or else this term makes no sense).   Also note that the penalty term is definitely non-negative due to the fact that the second derivative is being squared (but other ways of doing this, like using absolute value, can and are also considered, with different results).
% Let us now consider what kind of function f might be found that minimizes the objective above.  If we were to ignore the RSS term completely, the penalty term by itself would be most happy with any function that has zero second derivative everywhere: a straight line.  If we take \lambda to be very large, then we will end up with somethign that is pretty much a straight line (since we mostly care about the penalty term), but of all straight lines the one with the lowest RSS will be the final answer: this is just standard OLS regression.  On the other hand, if \lambda =0, this objective will just favor any function that completely interpolates the data (passes through every datapoint).  Interestingly, one can show that for any finite \lambda, the function that minimizes the above objective is a natural cubic spline (which is twice differentiable) with K=n knots located at each of the x values of our dataset!  Part of the proof of this is noting, as stated above, that such a spline can perfectly interpolate a sequence of points.  Hence, for any proposed function that minimizes the objective above on a dataset (x_i,y_i) by passing through the sequence of points (x_i,\hat{y}_i), a natural cubic spline can be used that passes through the same points, thus giving the same RSS.  The other part involves showing that no other function could possibly have a smaller penalty term than the natural cubic spline: see course notes for that part.
% Given this, one can rewrite the minimization problem in matrix vector form.  Let us call the natrual cubic spline basis functions N_j \left(x\right), for j=1,2,\ldots,n, and define the matrix \mathbf{N} such that entry \mathbf{N}_{ij}=N_j(x_i) (this is equivalent to our usual X matrix).  Then, given that our fitted function f shall be f(x)=\sum_j\beta_jN_j(x), the penalty term natrually induces a matrix \Omega, where \Omega_{ij}=\int N_i''(t)N_j''(t)dt.  Then the objective is to select the coefficient vector \beta to minimize:
% M=(y-\mathbf{N}\beta)^T(y-\mathbf{N}\beta)+\lambda\beta^T\Omega\beta
% Note that this objective function is extremely similar to that of ridge regression, except for the introduction of the matrix \Omega, which does not depend on the coefficient vector.  Then we can solve for the minimum analytically as we did in ridge regression, to find that the best fitting coefficients are given by:
% \hat{\beta}=\left(\mathbf{N}^T\mathbf{N}+\lambda\Omega)^{-1}\mathbf{N}^Ty
% The code below implements this particular smoothing spline.

%the below are basis function or their building blocks, plus second
%derivatives of such
h=@(x,c) (max(x-c,zeros(1,size(x,2)))).^3;
d2h=@(x,c) 6*max(x-c,zeros(1,size(x,2)));
d=@(x,c,xK) (h(x,c)-h(x,xK))./(xK-c);
d2d=@(x,c,xK) (d2h(x,c)-d2h(x,xK))./(xK-c);
basis=@(x,c,ind,K) (ind==1)*ones(1,size(x,2))+(ind==2)*x+(ind>2)*(d(x,c(max(ind-2,1)),c(K))-d(x,c(K-1),c(K)));
d2basis=@(x,c,ind,K) (ind<=2)*zeros(1,size(x,2))+(ind>2)*(d2d(x,c(max(ind-2,1)),c(K))-d2d(x,c(K-1),c(K)));


K=n;
nmat=zeros(K); %storage for our matrix N, then fill it in
for j=1:K
    tmp=basis(xvals',xvals,j,K);
    nmat(:,j)=tmp';
end
omega=zeros(K); %storage for our matrix Omega, then fill it in
for i=1:K
    for j=1:K
        omega(i,j)=integral(@(x) d2basis(x,xvals,i,K).*d2basis(x,xvals,j,K),0,10);
    end
end

lambda=.01; %pick whatever lambda you like
A=nmat'*nmat+lambda*omega;
beta=A\(nmat'*yvals); %best fitting coefficients
yhat=nmat*beta;
Smat=nmat*inv(A)*nmat'; %this matrix allows us to determine the effective degress of freedom for this smoothing spline
dof=trace(Smat) %effective degrees of freedom

%visualize
xeven=xa:.01:xb;
bplot=zeros(size(xeven,2),K);
for j=1:K
    tmp=basis(xeven,xvals,j,K);
    bplot(:,j)=tmp';
end
yspline=bplot*beta;
plot(xvals,yvals,'bo',xeven,yspline)

% Try some different values for \lambda above.  As you increase it, you will obtain functions that are more and more linear looking, while decreasing it will cause your function to get closer to interpolating the datapoints.
%
% The below section tries several different lambda values and computes the LOOCV error for each.  You can use this to pick the best value!

luse=0.001:.005:.1;
nl=size(luse,2);
cverr=zeros(1,nl);
dofs=cverr;
for j=1:nl
    lambda=luse(j);
    A=nmat'*nmat+lambda*omega;
    beta=A\(nmat'*yvals);
    yhat=nmat*beta;
    Smat=nmat*inv(A)*nmat';
    dofs(j)=trace(Smat);
    tmp=(yvals-yhat)./(1-diag(Smat));
    cverr(j)=tmp'*tmp/K;
end
plot(luse,cverr)
plot(luse,dofs)