% ### EXresampleLoess2.m ###    11.16.14    {C. Bergevin}
% % Resampling approach to demonstrate means to estimate uncertainty for
% non-parametric locally-weighted regression (loess) trends
% --> this version allows user to specify whether to use replacement or not
% [e.g., for bootstrapping, set range= 1 and replace= 1]
% [if so, then this would be 'bootstrapping']
% [see also EXloess1.m; requires loess.m and boundedline.m]

clear;
% -----------------------------
N= 100;     % # of points for fit (over interval [min(x) max(x)])
xR= [0 1];  % range of x-values
scaleN= 0.3;    % scale factor for noise
alpha= 0.17;     % loess fit parameter (between 0 and 1)
order= 1;        % order of polynomial for fit
M= 20;  % # of resamples
range= 1; % range to vary span over [0,1]
replace= 1; % use replacement? [0-no, 1-yes]
% -------------------------------------------------
x= linspace(min(xR),max(xR),N); % determine fit x-values
polyS= ceil(10*rand(1));   % randomly determine polynomial order
Pv= 2*rand(polyS,1)-1;  % polynomial coefficients
y= polyval(Pv,x)+ rand(1)*sin(2*pi*x)+ scaleN*randn(N,1)'; % pseudo-random 'data'
N= numel(x);   % total number of points
% ---
% loess fit to ENTIRE data set
xFit= linspace(min(xR),max(xR),N); % can 'resample' if desired (i.e., need not have xFit=x)
yFitW= loess(x,y,xFit,alpha,order,[],1,0,0.1);   % loess fit via external function
% ------------
% Resampling: regression on randomly chosen subsets of the data
for nn=1:M
    % create randomresampling index
    if replace==0
        temp= randperm(N);  % each entry is unique
    elseif replace==1
        temp= floor(N*rand(1,N))+1;  % allow for possible repeated values
    end
    tINDX= temp(1:ceil(range*N));   % trim to subset of specified size
    temp2(nn,:)= loess(x(tINDX),y(tINDX),xFit,alpha,order,[],1,0,0.1);
end
yFitBS= mean(temp2);
errBS= std(temp2);
% ---
% visualize
figure(1); clf; hold on; grid on;
boundedline(xFit,yFitBS,2*errBS,'b');   % trend from bootstrapping
%plot(xFit,yFitBS,'b-','LineWidth',2);    % without boundedline.m, uncomment these lines
%plot(xFit,yFitBS-2*errBS,'b-'); plot(xFit,yFitBS+2*errBS,'b-'); % 95% CIs for bootstrapped trend
plot(xFit,yFitW,'r-','LineWidth',2);    % trend for whole data set
plot(x,y,'ko');
xlabel('x'); ylabel('y');
legend('resampled CI','mean resampled trend','trend to whole set','data','Location','SouthWest');