Loess local regression

Acknowledgements

This code was jointly developed by Luci Ellis, RBA,  and David Burmaster, Alceon Corp, 1997–98.  See "Visualising Data" for more information about the procedure.  The procedure requires the simpleLinRegress[] function developed by Luci Ellis, from the econometric_functions.nb notebook.

Startup

The cell below just loads in some Mathematica add-ons.

Needs["Statistics`NormalDistribution`"];
Needs["Statistics`ContinuousDistributions`"];
Needs["Statistics`DataManipulation`"];
Needs["Statistics`DescriptiveStatistics`"];
Needs["Statistics`MultiDescriptiveStatistics`"];
Needs["Statistics`LinearRegression`"];
Needs["Statistics`DataSmoothing`"];
Needs["Graphics`Graphics`"];

EGloess::usage  =  
      "This function implements the fixed-window (variable sample size) \
version of the loess local regression procedure where there is a single \
explanator.  Loess is useful when there are few explanatory variables, a \
fairly large sample size, but the functional form is unknown and potentially \
complicated."

EGloessX::usage  =  
      "This function implements the fixed-sample-size (nearest neighbour) \
version of the loess local regression procedure where there is a single \
explanator.  Loess is useful when there are few explanatory variables, a \
fairly large sample size, but the functional form is unknown and potentially \
complicated."

EGloess[	xypairs:{..{_Real,_Real}},
			nwtfn_			/; 1 <= nwtfn,
			npoly_Integer	/; 1  == npoly || 2 == npoly,
			ndiv_Integer	/; 10 <= ndiv,
			swin_			/; 0  <  swin		] :=
			(*the conditions for the function to be valid*)
Module[{xvec, yvec, xmin, xmax, wtfn, yloesspt,
		xloessvec, yloessvec, wtvecvec,xypairssorted,regressdata},
		(* these are the names for the local variables *)
		
xypairssorted = Sort[xypairs];
{xvec, yvec}	=	Transpose[xypairssorted];
{xmin, xmax}	=	{First[xvec], Last[xvec]};
(* all pretty straight-forward *)

xloessvec =	
(N[xmin + # ((xmax-xmin)/(ndiv-1))]&) /@ Range[0,ndiv-1];
(* these are the points at which the fit is evaluated *)
wtvecvec  =	
	Table[N[Exp[-Abs[((xvec[[i]]-xloessvec[[j]])/swin)^nwtfn]]],
	{j, 1, ndiv},{i, 1, Length[xypairs]}];
(* a matrix where the rows are the vector of weights for each evaluation of the WLS fit *)

regressdata = If[npoly==1,Reverse /@ xypairssorted,Transpose[Join[{yvec},{xvec},{xvec}^2]]];
(* puts the data in a form appropriate for use in the simpleLinRegress[] auxiliary function *)

yloessvec=	Table[(Table[q^k,{k,0,npoly}]. 	
		Evaluate[simpleLinRegress[ regressdata,   	
		GLSWeights->wtvecvec[[jj]] ]][[1]]) /. q->xloessvec[[jj]],
				 {jj, 1, ndiv}];
(* finds the WLS-fitted values for y at gridpoints of x *)

Transpose[{xloessvec, yloessvec}]
(* returns vector of ndiv {x,y} pairs for loess points *)
]

EGloessX[	xypairs:{..{_Real,_Real}},
			nwtfn_			/; 1 <= nwtfn,
			npoly_Integer	/; 1  == npoly || 2 == npoly,
			ndiv_Integer    /; 10 <= ndiv,
			samplesize_Integer	/; 20 <= samplesize,
			swin_			/; 0  <  swin		] :=
			(*the conditions for the function to be valid*)
Module[{xvec, yvec, xmin, xmax, wtfn, yloesspt,
		xloessvec, yloessvec, wtvecvec,xypairssorted,baseregressdata,regressdata},
		(* these are the names for the local variables *)
		
xypairssorted = Sort[xypairs];
{xvec, yvec}	=	Transpose[xypairssorted];
{xmin, xmax}	=	{First[xvec], Last[xvec]};
(* all pretty straight-forward *)

xloessvec =	
(N[xmin + # ((xmax-xmin)/(ndiv-1))]&) /@ Range[0,ndiv-1];
(* these are the points at which the fit is evaluated *)
wtvecvec  =	
	Table[N[Exp[-Abs[((xvec[[i]]-xloessvec[[j]])/swin)^nwtfn]]],
	{j, 1, ndiv},{i, 1, Length[xypairs]}];
(* a matrix where the rows are the vector of weights for each evaluation of the WLS fit *)

baseregressdata = If[npoly==1,Reverse /@ xypairssorted, Transpose[Join[{yvec},{xvec},{xvec}^2]]];
(* puts the data in a form appropriate for use in the simpleLinRegress[] auxiliary function *)

regressdata = If[npoly==1, 
(Take[#,samplesize]&)/@  (RotateLeft[baseregressdata, Max[0,Floor[Position[xvec, Floor[Select[xvec,Min[Abs[xvec-#]] ]] ]-(samplesize/2)]] ] &/@xloessvec);
(* Partitions the data into multiple samples of the desired size, adjusting for end-sample constraints  *)

yloessvec=	Table[(Table[q^k,{k,0,npoly}]. 	
		Evaluate[simpleLinRegress[ regressdata,   	
		GLSWeights->wtvecvec[[jj]] ]][[1]]) /. q->xloessvec[[jj]],
				 {jj, 1, ndiv}];
(* finds the WLS-fitted values for y at gridpoints of x *)

Transpose[{xloessvec, yloessvec}]
(* returns vector of ndiv {x,y} pairs for loess points *)
]