BeginPackage["DataTools`",
	"Utilities`FilterOptions`",
	"Graphics`Graphics`",
	"Statistics`DataManipulation`",
	"Statistics`InverseStatisticalFunctions`"]

Average::usage = "Average[x_] computes average of list x.";

Bin::usage = "Bin[x_List, y_List] returns a nested list of the values of 
y which share a common value of x. The first element in each item of 
the result is the binning value.";

BinPercentage::usage = "BinPercentage controls the proportion of data 
binned by Bin.";

BinWidth::usage = "BinWidth sets an absolute window size for binning data.";

Boxplot::usage = "Boxplot[x_, opts_] draws either horizontal or vertical 
boxplots with optional shading.";

BoxVertical::usage = "BoxVertical is an option to boxplot which determines 
vertical or horzontal layout.";

BoxWidth::usage = "BoxWidth sets the non-semantic width of boxplot.";

BoxCenter::usage = "BoxCenter locates the box on the non-semantic axis.";

BoxFill::usage = "BoxFill determines whether a boxplot is shaded."

Correlation::usage = "Correlation[x_,y_] computes the correlation between 
lists x and y.";

Histogram::usage = "Histogram[x_List, opts] displays a histogram of 
the list x.";

AsDensity::usage = "AsDensity forces a histogram to be normalized to 
integrate to 1."

InterquartileRange::usage = "InterQuartileRange[x_] computes 
interquartile range of list x.";

KernelDensity::usage = "KernelDensity[x_List, options] creates an 
interpolating function which is the kernel density estimate.";

KernelFunction::usage = "KernelFunction option defines kernel used in 
density estimation.";

KernelGrid::usage = "KernelGrid determines number of points used in 
kernel estimation.";

KernelScale::usage = "KernelScale is the width factor used in kernel 
density estimation.";
	
BiweightKernel::usage = "BiweightKernel is a possible value for 
KernelFunction.";

GaussianKernel::usage="GaussianKernel is a possible value for KernelFunction.";

MAD::usage = "MAD[x_] computes median absolute deviation about the 
median of list x.";

Median::usage = "Median[x_] computes median of list x.";

NormalSample::usage = "NormalSample[n_?NumberQ] returns a list of 
n standard N[0,1] r.v.  NormalSample[n_, s2_] returns a list of n 
with variance s2.";

OLS::usage="OLS[x_,y_] returns the fitted values and residuals from 
a least squares regression of y on x.";

Residuals::usage = "Residuals determines if regression should compute 
residuals or return the fitted parameters.";

ReadDataFile::usage = "ReadDataFile[nameString] reads a file of data 
and defines the associated columns.";

RobustRegression::usage = "RobustRegression[x_List, y_, opts___] builds 
a robust regression model and prints a summary of the fit. By default, 
the program uses a biweight influence function.";

RobustConstant::usage = "RobustConstant sets returns the robustness 
constant used in estimating the coefficients in a robust regression.";

RobustTolerance::usage = "RobustTolerance is maximum percentage change 
permitted for the robust regression iterations to converge.";

RobustWeightFunction::usage = "RobustWeightFunction is an option of 
Robust regression used to set the weight function in robust regression."

HuberWeightFunction::usage = "HuberWeightFunction is the Huber weight 
function for a robust regression.";

BiweightWeightFunction::usage = "BiweightWeightFunction is the biweight 
weight function for a robust regression.";
	
Scatterplot::usage = "Scatterplot[x_, y_] displays a scatterplot of y on x.";

ScatterplotFunction::usage = "ScatterplotFunction defines a function 
to be superimposed on a scatterplot.";

ScatterplotMatrix::usage = "ScatterplotMatrix[x_,...] draws a scatterplot 
matrix ofthe collection of input lists.";

Smooth::usage = "Smooth[x_,y_,options_] generates a piecewise linear 
scatterplot smoother.";

SmoothSkipPercentage::usage = "SmoothSkipPercentage determines proportion 
of data skiped between evaluations.";

SmoothBinPercentage::usage = "SmoothBinPercentage sets the proportion 
of data used for linear fit.";

StandardDeviation::usage = "StandardDeviation[x_] computes standard 
deviation of list x.";

Quantile::usage = "Quantile[x_List,q_] returns qth quantile from data 
in list x.";

QQPlot::usage = "QQPlot[x_List] generates a normal quantile-quantile 
plot with reference line.";

Variance::usage = "Variance[x_List] computes the sample variance of input 
data x.";

WriteDataFile::usage = "WriteDataFile[syms_List, file_String, doc_] writes 
the lists identified by the symbols as a data file with an optional 
documentation string.";

Begin["`Private`"]

MakeLabel[x_] := StringDrop[StringDrop[ ToString[x] ,-1],5]

Attributes[getNames] = {HoldAll, Listable};
getNames[n_] := HoldForm[n];

(* filtering data...  *)

okQ[x_] := NumberQ[x] || StringQ[x];
okListQ[x:{_?okQ...}] := True;
okListQ[_] := False;

dataFilter::found = "`1` missing values removed.";

dataFilter[x_?MatrixQ] := With	\
[	{n = Length[x[[1]]], 
	 y = Transpose[ Select[Transpose[x],okListQ] ]
	},
	If[Length[y[[1]]]==n,
		x,
		Message[dataFilter::found,n-Length[y[[1]]] ]; y
	]
]

(* input already as list of rows *)
dataFilter[x_?MatrixQ, asRows] := With	\
[	{n = Length[x], 
	 y = Select[x,okListQ]
	},
	If[Length[y]==n,
		x,
		Message[dataFilter::found,n-Length[y] ]; y
	]
]

dataFilter[x_List] := With	\
[ 	{n=Length[x],y=Select[x,okQ]
	},
	If[Length[y]==n,
		x,
		Message[dataFilter::found,n-Length[y]]; y
	]
]

(* Reading data files.............  15 Jun 95  ... 24 Jun 95 *)

ReadDataFile::defined = "Variable symbol(s) `1` already defined. 
Clear and read again.";

ReadDataFile[file_String] := Module	\
[	{varSymbols, data, values, preDefined, input = OpenRead[file] },
	If[input===$Failed, Return[]];
	Print["File documentation: ", Read[input, String]];
	data = ReadList[input, Word, RecordLists->True];
	varSymbols = (ToHeldExpression/@First[data]) /. Hold[x_]->HoldForm[x];
	preDefined = Select[varSymbols, Head[#[[1]] ]=!=Symbol&];
	If[{}===preDefined,
		varSymbols = #[[1]]& /@ varSymbols; (* remove Hold *)
		Print["Variables to be read: ", varSymbols];
		values = Transpose[(ToExpression /@ #)& /@ Drop[data,1]];
		Print[Length[First[values]], " cases were read."];
		MapThread[Set[#1,#2]&, {varSymbols, values}],
	(* else *)
		Message[ReadDataFile::defined, preDefined];
	];
	Close[input];
]

WriteDataFile::length = "Variables are not columns of equal length.";

Attributes[WriteDataFile] = {HoldFirst};

WriteDataFile[syms_List, file_String, doc_:" "] := Module	\
[	{k,i, lens, types, w,
	 varSymbols = getNames[syms],
	 data = syms,
	 output = OpenWrite[file, FormatType->OutputForm]
	},
	If[output===$Failed, Return[]];
	types = Union[Head /@ data];
	lens = Union[Length /@ data];
	If[ Length[lens]>1 || Length[types]>1 || First[types]=!=List,
		Message[WriteDataFile::length]; Close[output]; Return[],
		lens = First[lens]];
	k = Length[varSymbols];
	WriteString[output, "\"",doc,"\"\n"];
	(* Do[WriteString[output,varSymbols[[i]]," "], {i,k}];
	   WriteString[output,"\n"];
	   w[{x__}] := Write[output,x];
	   w /@ Transpose[data];
	*)
	Write[output, TableForm[Transpose[data], 
		TableHeadings->{None,varSymbols}, TableSpacing->{0,2}]];
	Close[output];
]

(* Random Number Generation.......................... 10 Nov 94 *)

NormalSample[n_Integer] := Table[nr[1],{n}]

NormalSample[n_Integer,s2_?NumberQ /; s2>0] := Table[nr[s2],{n}]

	pi2 = N[2 Pi];
	nr = Compile[{s2},
		Sqrt[-s2 2 Log[ Random[] ] ]
			Cos[ pi2 Random[] ]
		];
				
Average[x_List] := Module	\
[	{n, y=dataFilter[x],
	},
	n=Length[y];
	(Table[1,{n}].y)/n
]

Average[x_List,None] := With[{n=Length[x]},(Table[1,{n}].x)/n]

Median[x_List] := Module	\
[	{n, sx = Sort[dataFilter[x]]},
	n = Length[sx];
	(Plus @@ sx[[Floor[(n + {1,2})/2]]])/2
]
	
MAD[x_List] := With \
[	{y=dataFilter[x]},
	Median[Abs[y-Median[y]]]
]



Variance[x_List] := Module	\
[	{y = dataFilter[x], dev},
	dev = y - Average[y,None];
	(dev.dev) / (Length[dev]-1)
]

StandardDeviation[x:_List] := Sqrt[Variance[x]];

(* Quantile assumes sorted *)
		
Quantile::input = "Quantile expects a sorted numeric list with second 
input between 0 and 1.";

Quantile[sx:{_?NumberQ...}, q_?NumberQ/;(0<q)&&(q<1)] := With 	\
[	{n = Length[sx], f=Mod[q*Length[sx],1]},
	{sx[[ Max[1,Floor[q*(n+1)]] ]],sx[[ Min[n,Ceiling[q*(n+1)]] ]]} . {1-f,f}
]

Quantile[x_List,q_] := Message[Quantile::input];
			
InterquartileRange[x:{_?NumberQ...}] := Quantile[x,.75] - Quantile[x,.25];

Correlation[x_List, y_List] := Module 	\
[	{xy = dataFilter[{x,y}],
     dx,dy
    },
    dx = xy[[1]] - Average[xy[[1]],None];
    dy = xy[[2]] - Average[xy[[2]],None];
	(dx . dy) / Sqrt[(dx.dx)(dy.dy)]
]

Options[Bin] = {BinPercentage->.2, BinWidth->Automatic};

Bin::TooMany = "Appear to have too many (`1`) bins.";

(* collect y values with common x *)
Bin[ x_List, y_List ] := With	\
[	{ux = Union[x]},
	If[ (Length[ux] > .9 Length[x]) && (Length[ux] > 10),
		Message[Bin::TooMany, Length[ux]];
		{x,y},
	(* else *)
		{#,y[[ Flatten[Position[x,#]] ]]}& /@ ux
	]
]

(* collect values of y with close values of x *)
Bin[x_List,y_List, opts___?OptionQ] := Module	\
[	{in = If[OrderedQ[x],
			Range[Length[x]],
			Sort[Range[Length[x]], (x[[#1]]<x[[#2]])&]],
	 pct, width, ni, sx, binLim, n=Length[x],i, j
	},
	{pct,width} = {BinPercentage, BinWidth} /. opts /. Options[Bin];
	sx = x[[in]];
	If[ Automatic===width, (* constant number ni in each bin *)
		ni = Floor[pct*n];
		Transpose[{Average /@ Partition[sx, ni],
					Partition[y[[in]],ni]}],
	(* else constant width *)
		ni = BinCounts[x, {sx[[1]], sx[[n]]+.001, width}];
		Transpose[{sx[[1]] + (.5 * width) + width * Range[0,Length[ni]-1],
				   y[[#]]& /@ (FoldList[Plus,0,Drop[ni,-1]]+Range/@ni)}]
		]
	]

Attributes[Histogram] = {HoldFirst};

Options[Histogram] = {BinWidth->Automatic, AsDensity->True};

Histogram[expr_, opts___?OptionQ] := Module		\
[	{ held = Hold[expr],
	  x = dataFilter[expr],
	  width, asDen, min, max, eps
	},
	{width, asDen} = {BinWidth,AsDensity} /. {opts} /. Options[Histogram];
	min = Min[x]; max=Max[x]; eps = 0.1 * (max-min);
	If[width===Automatic, 
		width = (max-min)/N[Ceiling[5 Sqrt[Length[x]]]]];
	hist[x, width, asDen, MakeLabel[held], min-eps, max+eps, opts]
]

hist[data_, dx_, scaled_, title_, min_, max_, opts___?OptionQ] := Module	\
[	{ n = Length[data], nBins=Ceiling[(max-min)/dx],
	  left, bin, p
	},
	left = min + dx * Range[0,nBins-1];
	bin = BinCounts[data, {min,max, dx}];
	If[Length[bin]!=nBins,	bin = Drop[bin, -1]]; (* force length *)
	(* Print["Bin (",nBins,") counts: ", N[bin]]; *)
	If[scaled, bin = bin / (n * dx)];
	p = GeneralizedBarChart[
		Transpose[{left, bin, Table[dx, {nBins}]}],
		opts,
		BarStyle->{GrayLevel[.9]},
		PlotLabel->title, PlotRange->All, AxesOrigin->{min,0}  ]
]

Options[KernelDensity] = {KernelFunction->Biweight,
						  KernelGrid->25,
						  KernelScale->Automatic};


KernelScale[Biweight] := 2.77794;
KernelScale[Gauss] := 1.05922;
KernelScale[f_, a_:-Infinity, b_:Infinity] := Module	\
[	{s, t, n, k2, num, g, dg2},
	g[x_] := (1/Sqrt[2 Pi]) E^(- (x^2)/2);
	dg2 = D[g[t],{t,2}];
	k2 = Integrate[ t^2 f[t],{t,a,b}];
	num = Integrate[ f[t]^2, {t,a,b}];
	den = Integrate[ dg2^2, {t,-Infinity, Infinity}];
	N[ PowerExpand[ k2^(-2/5) (num / den)^(1/5)] ]
]

KernelFunction[Biweight] = If[Abs[#]<1, N[(15/16)(1-#^2)^2], 0]& ;
KernelFunction[Gauss]   = N[E^(-(#^2)/2) / Sqrt[2 Pi]]&;
KernelFunction[f_] := f;

KernelDensity::scale = "Kernel smoothing constant = `1` s/n^(1/5)";
KernelDensity::range = "Kernel width is `1` on interval [`2`,`3`]";

KernelDensity[expr_, opts___Rule] := Module	\
[	{h, scale, kernel, grid, kFunc, min, max, tab,one, n,
	 data = dataFilter[expr]
	},
	n = Length[data];
	{h, grid, kernel} = { KernelScale, KernelGrid, KernelFunction } /. 
				{opts} /. Options[KernelDensity];
	kFunc = KernelFunction[kernel];
	If[h === Automatic, h = KernelScale[kernel]];
	Message[KernelDensity::scale, h];
	scale = h N[Min[ StandardDeviation[data],
	               InterquartileRange[data]/1.345] / n^(1/5)];
	If[NumberQ[grid], 
		min = Min[data] - 2 scale;
		max = Max[data] + 2 scale;
		grid= Range[min, max, (max-min)/grid],
	(* else *)
		min = First[grid];
		max = Last[grid];
	];
	Message[KernelDensity::range,scale,min,max];
	one = Table[1,{n}];
	tab = Table[
			  (one . (kFunc /@ ((data - grid[[i]])/scale) ))/scale,
			  {i,1,Length[grid]}]
		  / n;
	Interpolation[Transpose[{grid,tab}]]
]

Attributes[QQPlot] = {HoldFirst};

QQPlot[expr_] := Module \
[	{data = Sort[dataFilter[Release[expr]]],
	 label = MakeLabel[Hold[expr]],
	 m,s,n,os,i},
    m = Average[data];
    s = StandardDeviation[data];
    n = Length[data];
    os = Sqrt[2] Table[InverseErf[N[2 (i/(n+1))-1]],{i,Length[data]}];
    ListPlot[Transpose[{os,data}],
			AxesLabel->{"Normal", label},
			Epilog->Line[{{os[[1]], m + s os[[1]]},
						  {os[[n]], m + s os[[n]]}}]  ]
]

Options[Boxplot] = 
	{BoxVertical->False, BoxWidth->1, BoxCenter->1,
	 BoxFill->True, AspectRatio->0.5, Frame->True};
					 
Boxplot[inx_List, opts___?OptionQ ] := Module	\
[	{box, vert, center, width, ar, frame, b, bi, fill, gOpts, ticks,
	 x = dataFilter[inx]},
	{fill, vert,width, center, ar, frame} = {BoxFill, BoxVertical, BoxWidth, BoxCenter, AspectRatio, Frame}
	 /. {opts} /. Options[Boxplot];
	If[MatrixQ[x],
		b = Length[x]; bi = Range[b];
		box = MapThread[boxplot[#1,width,center+#2, fill] &, {x, 2*width*Range[0,b-1]}];
		box = {{Min[ box[[bi,1,1]] ], Max[ box[[bi,1,2]] ]}, Join@@box[[bi,2]]},
	(* else *)
		b = 1;
		box = boxplot[x, width, center, fill] ];
	If[vert, 
		ar=1/ar; ticks={None, Automatic},
		ticks={Automatic, None}];
	gOpts = FilterOptions[Plot, opts];
	Show[
	  Graphics[If[vert, box[[2]] /. {a_?NumberQ ,b_} :> {b,a}, box[[2]] ],
	  			Frame->frame, AspectRatio->ar, FrameTicks->ticks,
				gOpts,
	  			PlotRange->{{center-2*width, center+2*b*width},box[[1]]}[[ If[vert,{1,2},{2,1}] ]]
	  ]]
]

(* assumes sorted *)
fiveNum[ sx_List ] := With	\
[	{n=Length[sx]},
	Join[{sx[[1]]}, Quantile[sx,#]& /@ {.25,.5,.75}, {sx[[n]]}]
]

(* Returns {min,max} and graphics list *)
boxplot[x_List, width_, center_, fill_] := Module		\
[	{sx = Sort[x], w=width/2, n=Length[x],
		lo, hi, fn, i, iqr, lowAt, hiAt},
	fn = fiveNum[sx];
	iqr = fn[[4]]-fn[[2]];
	lo = fn[[2]] - iqr * 1.5;
	hi = fn[[4]] + iqr * 1.5;
	i = 1;
	While[ sx[[i]] < lo, ++i]; lowAt = i; lo = sx[[lowAt]];
	i = n;
	While[ sx[[i]] > hi, --i]; hiAt = i;  hi = sx[[hiAt]];
	{{fn[[1]]-.1*iqr, fn[[5]]+.1*iqr},
	Join[If[fill,{	GrayLevel[0.9],
		  			Polygon[{{fn[[2]],center-w}, {fn[[4]],center-w},
		  			 		 {fn[[4]],center+w}, {fn[[2]],center+w}}],
		  			GrayLevel[0.0]}, {}],
		 {
			Line[{{fn[[2]],center},{lo,center}}],	(* whiskers *)
			Line[{{fn[[4]],center},{hi,center}}],
			Line[{{fn[[3]],center-w},{fn[[3]],center+w}}],	(* median line *)
			Line[{{fn[[2]],center-w}, {fn[[4]],center-w},	(* border *)
		  			 {fn[[4]],center+w}, {fn[[2]],center+w},
		  			 {fn[[2]],center-w}}],
		  	AbsolutePointSize[4]
		 },
	 		Point[{#,center}]& /@ sx[[ Range[1, lowAt-1] ]],
			Point[{#,center}]& /@ sx[[ Range[hiAt+1, n] ]]
		]}
]

Options[Scatterplot] = Join[{ScatterplotFunction->None}, Options[ListPlot]];

Scatterplot[{x_, y_}, opts___?OptionQ] := Scatterplot[x,y,opts]

scat[plotFunction_, data_, opts___] := Module \
[	{pf, lp, x, 
     dispFunc = DisplayFunction /.{opts} /. Options[Scatterplot],
     f = ScatterplotFunction /. {opts} /. {ScatterplotFunction->None}
    },
    lp = plotFunction[data, 
    	PlotRange->All, DisplayFunction->Identity,
    	FilterOptions[plotFunction, opts]
	];
	If[f === None,
		Show[lp, DisplayFunction->dispFunc],
	(* else *)
		x = First/@data;
		pf = Plot[f[z],{z,Min[x], Max[x]}, DisplayFunction->Identity];
		Show[lp,pf,DisplayFunction->dispFunc]
	]
]

Scatterplot[x_List, opts___?OptionQ] := \
	scat[ListPlot, dataFilter[Transpose[{Range[Length[x]],x}]], opts]

Scatterplot[x_List, y_List, opts___?OptionQ] := 
	scat[ListPlot, dataFilter[Transpose[{x,y}], asRows], opts];
	
Scatterplot[x_List, y_List, labels_List, opts___?OptionQ] := 
	scat[LabeledListPlot, dataFilter[Transpose[{x,y,labels}], asRows], opts];


Scatterplot3D[x_List, y_List, z_List, opts___?OptionQ] := 
	ListPlot3D[
		dataFilter[Transpose[{x,y,z}], asRows],
		opts]
		
Attributes[ScatterplotMatrix] = {HoldAll};


ScatterplotMatrix[varSeq__] := Module	\
[	{i,j,held = Hold[varSeq], vars={varSeq}, k, c},
	k=Length[held];
	Show[
	  GraphicsArray
		[	Table[
				If[j>i,
					Scatterplot[vars[[j]],vars[[i]],
						AxesLabel->None, Axes->None,
						DisplayFunction->Identity,
						Frame->True, FrameTicks->None],
				(* else show correlation or name *)
					If[j==i,
						Graphics[ Text[MakeLabel[HeldPart[held,i]],{0,0}] ],
						c = Correlation[vars[[i]], vars[[j]]];
						Graphics[Text[N[c,2], {0,0}]]
					]
				],
				{i,1,k},{j,1,k}
			]
		]
	]]

(* Smoothing...................................................... 13 Jun 95
   Piecewise linear smoothing using either of two types of binning:  Take a 
   contant percentage of the data or by a constant distance (a la kernel).
   Reference to faster algorithms from JCGS.  Alternative to Fourier based 
   methods. 
*)

Options[Smooth] = {SmoothSkipPercentage->.05, SmoothBinPercentage->.40};

Smooth[x_List, y_List, opts___?OptionQ] := Module	\
[	{skip, pct, in, sx, sy, n=Length[x]},
	{skip,pct} = {SmoothSkipPercentage, SmoothBinPercentage} /. {opts} /. Options[Smooth];
	{sx,sy} = dataFilter[{x,y}];
	If[Not[OrderedQ[x]],
		in = Sort[Range[n], (sx[[#1]] <= sx[[#2]])&];		(* sort indices *)
		sx = sx[[in]]; sy = sy[[in]];
	];
	smthByPct[sx, sy, Ceiling[pct*n], Ceiling[skip*n]]
]


fit[x_, y_, at_] := Module	\
[	{xBar = Average[x], xDev, b},
	xDev = x-xBar;
	b = First[{xDev.y}/{xDev.xDev}];
	If[ListQ[at],
		Transpose[{at,Average[y] + b (at - xBar)}],
	(* else *)
		{at,Average[y] + b (at - xBar)}
	]
]

(* assumes sorted on x *)
smthByPct[x_,y_, binSize_, skip_] := Module	\
[	{n = Length[x],
	 half = Floor[binSize/(2 skip)],   (* number at each end *)
	 xLo, xHi, max, min
	},
	min=x[[1]]; max = x[[n]];
	xLo = min + (Average[x[[ Range[binSize] ]]]- min)
	 				 Range[0,half-1]/half;
	xHi = max + (Average[x[[ n+1-Range[binSize] ]]] - max) 
						Reverse[Range[0,half-1]]/half;
	(* Print["binSize=", binSize, " skip=", skip]; *)
	Interpolation[
		Join[fit[x[[ Range[binSize] ]],y[[ Range[binSize] ]], xLo],
			 MapThread[fit[#1,#2,Average[#1]]&,
	 			{Partition[x,binSize,skip],
	 			 Partition[y,binSize,skip]}
	 		 ],
	 		 fit[x[[ n+1-Range[binSize] ]], y[[ n+1-Range[binSize] ]], xHi]
	 	]
	]
]

ols[x_, y_] := Block \
			[{q,r},
			 {q,r} = QRDecomposition[x];
			 Inverse[r] . q . y
			]
	
	wls[x_,y_,w_] := With[ {sr=Sqrt[w]},
							ols[sr x, sr y]
						 ]
	
	
	diag[m_] := Transpose[m,{1,1}]
	
	Options[OLS] := {Constant->True, Residuals->True};
	
    Attributes[OLS] = {HoldFirst};
    
    OLS[xNames_?MatrixQ, iny_List, opts___?OptionQ] := Module	\
	[	{q, r, ri, x = xNames, xtxi, xMat, b, z, se, s2, const, resids,
		 k, yx, y, n
		},
		yx = dataFilter[Join[{iny},x]];
		y = yx[[1]]; x = Drop[yx,1];
		n = Length[y];
		{const,resids} = {Constant, Residuals} /. {opts} /. Options[OLS];
		If[const, xMat = Prepend[x,Table[1,{n}]],
				  xMat = x];
		{q,r} = QRDecomposition[Transpose[xMat]];
		ri = Inverse[r];
		b = ri . q . y;
		k = Length[b];
		z = y - b . xMat;
		s2 = (z.z)/(n-k);
		se = N[Sqrt[s2 diag[ri . Transpose[ri]]]];
		Print[TableForm[Transpose[  {
			Join[{"Factor"}    , If[const,{1},{}]  , getNames[xNames] ],
			Join[{"Coefficient"} , b],
			Join[{"Std. Err."} , se],
			Join[{"t-Stat"}    , b/se] }  ]]];
		Print[""];
		Print["s = ",Sqrt[s2],
			"  R",Superscript[2]," = ",1-(z.z)/((n-1)Variance[y])];
		If[resids, {y-z,z}, b]
	]
	
	OLS[x_,y_, opts___?OptionQ] := OLS[{x},y, opts]

(* Robust Regression ................ 26 June 95 *)

(* Some weight functions *)

	Attributes[HuberWeightFunction] = {Listable};
	HuberWeightFunction[x_/; (Abs[x]<1)]	:= 1
	HuberWeightFunction[x_]        			:= 1 / Abs[x];
	
	Attributes[BiweightWeightFunction] = {Listable};
	BiweightWeightFunction[x_/; (Abs[x]<1)] := (1 - x^2)^2
	BiweightWeightFunction[x_] = 0;

(* Suggested scaling constants *)

	RobustConstant[HuberWeightFunction]    = 1.365;
	RobustConstant[BiweightWeightFunction] = 4.685;
	
(* Convergence test *)
	
	RobustTolerance = 0.001;
	
	convQ[x_,y_] := Max[ Abs[x-y]/x ] < RobustTolerance


Options[RobustRegression] = Join[Options[OLS],
	{RobustWeightFunction->BiweightWeightFunction, 
 	 RobustConstant->4.685,
 	 MaxIterations->10}];
							 	
Attributes[RobustRegression] = {HoldFirst};
			 
RobustRegression::iter  = "Iteration `1` @ `2`";
RobustRegression::itLim = "Number of iterations exceeds limit.";

RobustRegression[xNames_?MatrixQ, iny_, opts___?OptionQ] := Module	\
[	{rr, W, c, se, est, const, resids, itLim, yx,y,
	 x=xNames
	},
	yx = dataFilter[Join[{iny},x]];
	y = yx[[1]]; x = Drop[yx,1];
	{W, c, const, resids,itLim} =
			{RobustWeightFunction, RobustConstant, Constant, Residuals, MaxIterations}
			/. {opts} /. Options[RobustRegression];
	rr = robRegr[ If[const, Prepend[ x, Table[1,{Length[y]}]],x],
				  y, W, c, itLim, resids];
	est = Estimates /. rr;
	se = Sqrt[diag[Variance /. rr]];
	Print[""];
	Print[TableForm[Transpose[  {
		 Join[{"Factor"},If[const,{1},{}]  , getNames[xNames] ],
		 Join[{"Coefficient"}, est],
		 Join[{"Asym. S.E."} , se],
		 Join[{"Coef/S.E." } , est/se]}
	]]];
	Print[""]; Print["Robust scale = ", Scale/.rr];			
	rr
]

RobustRegression[x_, y_, opts___?OptionQ] := RobustRegression[{x},y,opts]

robRegr[Xt_, Y_, weightFunction_, c_, itLim_, full_:True] := Module \
[	{e, est, s, w, h, one, b, Mi, Q, itCount=0,k,nextB,
	 X = Transpose[Xt],
     n = Length[Y]
    },
    k = Length[First[X]];
    Si = Inverse[N[Xt . X]];
	h = N[(X  (X . Si)) . Table[1,{k}]];
	inf = Sqrt[1-h];
	nextB[b_] := 
	(	e = N[Y - X . b];				(* w, e are global to robRegr *)
		s = Median[ Abs[e] ]/.6745;
		w = weightFunction[ e / (inf (c s)) ];
		est = Inverse[N[Xt . (w X)]] . N[Xt . (w Y)];
		itCount += 1;
		Message[RobustRegression::iter, itCount, est];
		est
	);
	b = Si . Xt . Y;  
	Message[RobustRegression::iter, 0, b];
    b = FixedPoint[nextB, b, itLim, SameTest->convQ];
	If[ itCount>=itLim, Message[RobustRegression::itLim]];
	If[full,					(* lots of extra output *)
		Mi = Inverse[ (Xt . (w e e X)) ];
		Q  = Xt . (w w e e X);
		var= (s^4)  Mi . Q . Mi;	(* Note the s^4 *)
		{Estimates -> est, Variance->var, Scale->s, WeightFunction->weightFunction,
			FitResiduals-> e, PredictedResponse-> Y-e, RobustC->c,
			Weights->w},
	(* else lean and mean *)
		est
	]
];

(* Regression Trailer
   Closing Formalities
*)

End[];  (* Private *)
Protect[ Average, 
		 Bin, BinPercentage, BinWidth,
		 Boxplot, BoxWidth, BoxCenter, BoxFill, BoxVertical,
		 Correlation,
		 Histogram, AsDensity, 
		 InterquartileRange,
		 KernelDensity, KernelScale, KernelGrid, KernelFunction, 
		 	BiweightKernel, GaussianKernel,
		 MAD, Median, 
		 NormalSample,
		 OLS, Residuals,
		 ReadDataFile, 
		 RobustRegression, RobustConstant, RobustTolerance, RobustWeightFunction,
		 	BiweightWeightFunction, HuberWeightFunction,
		 Scatterplot, ScatterplotFunction, ScatterplotMatrix,
		 Smooth, SmoothSkipPercentage, SmoothBinPercentage,
		 StandardDeviation, 
		 QQPlot, Quantile, Variance]
EndPackage[];