% ### 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');