function [Bet0,LOFShp,LOFScl]=GPSS(Bss,Pnl,Rgh,YExc,XExc,HlfWdt,CVnGrp);
%**************************************************************************
%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
%**************************************************************************

%**************************************************************************
%GP starting solution
%
%Bss      nExc x p basis matrix
%Pnl         p x p penatly matrix   Pnl=(D' * D)
%Rgh         r x 1 values of parameter roughness to consider
%YExc     nExc x 1 responses (threshold exceedances only)
%XExc     nExc x 1 covariate values corresponding to YExc   
%HlfWdt      1 x 1 half-width (on covariate scale) for selection of subset to fit
%CVnGrp      1 x 1 number of CV groups
%
%Bet0       2p x 1 parameter estimate (shape then scale) for starting solution
%LOFShp      r x 1 predictive lack of fit for shape estimate
%LOFScl      r x 1 predictive lack of fit for scale estimate
%**************************************************************************

%Sort sample in ascending values of XExc
[~,tOrd]=sort(XExc); YExc=YExc(tOrd);
n=length(YExc); %sample size

%Stationary GP fits
GP0=NaN(n,2); %initialise parameter array for fits to subsets
for i=1:n; %loop over sample
    I=mod((i-HlfWdt:i+HlfWdt)-1,n)+1; %local interval with respect to covariate
    GP0(i,:)=gpfit(YExc(I)); %stationary GP fit
end;

%Smoothing spline through parameter estimates for shape
Bet0=mean(GP0(:,1))*ones(size(Pnl,1),1);
[BetShp,LOFShp]=CV(Rgh,@NLLLS,CVnGrp,GP0(:,1),Bss,Pnl,Bet0,1);

%Smoothing spline through parameter estimates for scale
Bet0=mean(GP0(:,2))*ones(size(Pnl,1),1);
[BetScl,LOFScl]=CV(Rgh,@NLLLS,CVnGrp,GP0(:,2),Bss,Pnl,Bet0,1);

%Output argument
Bet0=[BetShp;BetScl];

return;