function [CV,g,s2,edf] = smu_CV(x,y,alpha);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Compute cross validation score for a smoothing spline fit for a given

%  value of the smoothing parameter, alpha.

%

%  Use the Hutchinson and de Hoog algorithm described in Green & Silverman

%  (1994), sec 3.2.2  

%

%  INPUTS:

%     x        n x 1 vector    Values of explanatory variable.

%                              Must be distinct and ordered x_i < x_{i+1}

%                              Require n >= 3

%     y        n x 1 vector    Values of response variable.

%     alpha    scalar          Given value of the smoothing parameter

%

%  OUTPUTS:

%     CV        scalar         The cross-validation score

%     g         n x 1 vector   Values of spline at x:  g(x)

%     s2        scalar         Estimate of sigma^2, the residual variance

%     edf       scalar         The effective df = tr(I - A)

%

%  Written by Karen Chiswell on 23 March 2006

%  Modified:         27 March 2006              Karen Chiswell

%                          Compute estimate of sigma^2 and the effective DF

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



n = length(x);



%         Get smoothing spline fit for given value of alpha

[g,gam,h,ssM,sM,dM] = smu_spline(x,y,alpha);





%         Get the central bands of the inverse of R + alpha*Q'Q

[ssI,sI,dI] = inv_band5(ssM,sM,dM);



%         Compute the diagonal elements of Q*inv(R + alpha*Q'Q)*Q'

%         See equation (3.12) in Green and Silverman

ImA = zeros(n,1);

ImA(1) =  dI(1) / (h(1)*h(1));

ImA(2) =  dI(1) * (-1/h(1) - 1/h(2)) * (-1/h(1) - 1/h(2)) ...

           + dI(2) / (h(2)*h(2)) ...

           + 2 * sI(1) * (-1/h(1) - 1/h(2))  / h(2);

ImA(n-1) =  dI(n-3) / (h(n-2)*h(n-2)) ...

           + dI(n-2) * (-1/h(n-2) - 1/h(n-1)) * (-1/h(n-2) - 1/h(n-1)) ...

           + 2 * sI(n-3) * (-1/h(n-2) - 1/h(n-1)) / h(n-2);

ImA(n) =  dI(n-2) / (h(n-1)*h(n-1));

for i = 3:n-2;

    j = i - 1;

    ImA(i) =   dI(j-1) / (h(i-1)*h(i-1)) ...

             + dI(j) * (-1/h(i-1) - 1/h(i)) * (-1/h(i-1) - 1/h(i)) ...

             + dI(j+1) / (h(i)*h(i)) ...

             + 2 * sI(j-1) * (-1/h(i-1) - 1/h(i)) / h(i-1) ...

             + 2 * ssI(j-1) / (h(i-1)*h(i)) ...

             + 2 * sI(j) * (-1/h(i-1) - 1/h(i))  / h(i);

end   % loop on i



%          Multiply by alpha:  I - A = alpha*[Q*inv(R + alpha*Q'Q)*Q']

ImA = alpha*ImA;    



%          Compute CV = sum_{i=1}^n [yi - g(i)/ImA(i)]^2 / n

ymg = (y - g) ./ ImA;

CV = (ymg' * ymg) / n;



%          Compute effective degress of freedom = tr(I - A)

edf = sum(ImA);



%          Compute estimate of residual variance

s2 = sum((y - g).^2);

s2 = s2/edf;



%    end of smu_CV.m