function [OptBet,LOF]=CV(Rgh,LklFnc,CVnGrp,Y,Bss,Pnl,Bet0,SclFct,Tau)
%**************************************************************************
%Illustration of non-stationary extreme value modelling in MATLAB accompanying 
%Threshold modeling of nonstationary extremes
%Paul Northrop, Philip Jonathan and David Randell
%in 
%Extreme Value Modeling and Risk Analysis: Methods and Applications
%Dipak Dey and Jun Yan (Editors)
%2015
%**************************************************************************
%
%**************************************************************************
%CROSS-VALIDATION FOR PARAMETER ESTIMATION
%
%Rgh      1 x 1 Roughness coefficient (lambda in: PNLL = NLL + lambda * (beta'*Pnl*beta))
%LklFnc   text  Desciption of likelihood function to use (QR, Pss or GP)
%CVnGrp   1 x 1 Number of cross-validation groups
%Y        n x 1 Responses
%Bss      n x p Spline basis evaluated at the covariate values corresponding to responses
%Pnl      p x p Spline penality matrix
%Bet0     p x 1 Starting value for spline parameter vector
%SclFct   1 x 1 Scale factor for roughness of GP scale
%Tau      1 x 1 Non-exceedance probability for extreme value threshold 
%
%OptBet   p x 1 Optimal spline parameter vector
%LOF      p x 1 Cross-validated lack of fit (log likelihood summed over "left out" observations)
%
%Note that p for GP is double the size of p for QR and Pss
%**************************************************************************

%Set up
opts=optimset('Display','off','LargeScale','off');  %turn off output from fminsearch (or other optimisers if available)
nRgh=size(Rgh,1); %number of roughness coefficients
LOF=zeros(nRgh,1); %initialise output

%Allocate observations to cross-validation groups
ICV=CVGrp(length(Y),CVnGrp,'Random');

%Determine if MATLAB Optimization toolbox is available. If so, use better optimiser.
if license('test','Optimization_Toolbox')==1;
    fminfun=@fminunc; 
else;
    fminfun=@fminsearch; %simplex
end;

%Perform cross-validation
for iR=1:nRgh; %loop over roughness values
    fprintf('.'); %progress dot to standard output
    for iCV=1:CVnGrp; %loop over cross-validation groups
        tI=ICV==iCV; %assess goodness for tI, %estimate model for ~tI
        if exist('Tau','var')==0; %Pss, GP and LS
            OptBet=fminfun(@(Bet)LklFnc(Bet,Y(~tI),Bss(~tI,:))+NLLPnl(Bet,Rgh(iR),Pnl,SclFct),Bet0,opts); %estimate parameters using training data
            LOF(iR)=LOF(iR)+LklFnc(OptBet,Y(tI),Bss(tI,:));  %NLL for prediction
        else; %QR
            OptBet=fminfun(@(Bet)LklFnc(Bet,Y(~tI),Bss(~tI,:),Tau)+NLLPnl(Bet,Rgh(iR),Pnl,SclFct),Bet0,opts); %estimate parameters using training data
            LOF(iR)=LOF(iR)+LklFnc(OptBet,Y(tI),Bss(tI,:),Tau);  %NLL for prediction
        end;               
    end;
end;
fprintf('\n'); %newline to standard output
[~,I]=min(LOF); %locate minimum lack of fit
optRgh=Rgh(I); %specify optimal roughness 

%Estimate parameters for full data and optimal roughness
if exist('Tau','var')==0; %Pss, GP and LS
    OptBet=fminfun(@(Bet)LklFnc(Bet,Y,Bss)+NLLPnl(Bet,optRgh,Pnl,SclFct),Bet0,opts);    
else; %QR
    OptBet=fminfun(@(Bet)LklFnc(Bet,Y,Bss,Tau)+NLLPnl(Bet,optRgh,Pnl,SclFct),Bet0,opts);    
end;

return;