(*******************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which
typically will have the same name as this file except ending in
".nb" instead of ".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent
to using the Evaluate Initialization Cells menu command in the front
end.

DO NOT EDIT THIS FILE.  This entire file is regenerated
automatically each time the parent Notebook file is saved in the
Mathematica front end.  Any changes you make to this file will be
overwritten.
***********************************************************************)





Off[General::spell1]



<<Utilities`FilterOptions`

<<Graphics`Graphics`

<<Statistics`DataManipulation`

<<Statistics`StatisticsPlots`

<<Statistics`ContinuousDistributions`

<<Statistics`DescriptiveStatistics`

<<Statistics`MultiDescriptiveStatistics`





ReadData::notfound="File `1` not found.";

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

ReadData::doc="Document string: `1`";

ReadData::names = "Symbols defined in this file: `1`";

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

ReadData[file_String]:=
  Module[{varSymbols,data,values,preDefined,input=OpenRead[file]},
    If[input===$Failed,
      Message[ReadData::notfound,file]; Return[],
      Message[ReadData::doc,Read[input,String]]];
    data=ReadList[input,Word,RecordLists\[Rule]True];
    Message[ReadData::names,First[ data]];
    varSymbols=(ToHeldExpression/@First[data])/.Hold[x_]\[Rule]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[ReadData::defined,preDefined];];
    Close[input];]



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

WriteData::length="Columns lack same length.";

Attributes[WriteData]={HoldFirst};

WriteData[syms_List,file_String,doc_:" "]:=
  Module[{k,i,lens,types,w,varSymbols=getNames[syms],data=syms,
      output=OpenWrite[file,FormatType\[Rule]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\[Rule]{None,varSymbols},
        TableSpacing\[Rule]{0,2}]];
    Close[output];]





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.";

Options[Bin]={BinPercentage\[Rule].2,BinWidth\[Rule]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 similar 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,i,j,
      n=Length[x]},
    {pct,width}={BinPercentage,BinWidth}/.opts/.Options[Bin];
    sx=x[[in]];
    If[Automatic===width,
      (*constant number ni in each bin*)
      ni=Floor[pct*n];
      Transpose[{Mean/@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)}]
      ]]



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.";

\!\(\(\(Options[KernelDensity] = {KernelFunction \[Rule] Biweight, 
        KernelGrid \[Rule] 25, 
        KernelScale \[Rule] Automatic};\)\(\[IndentingNewLine]\)
  \)\n
  \(KernelScale[Biweight] := 2.77794;\)\n
  \(KernelScale[Gauss] := 1.05922;\)\n
  \(\(KernelScale[f_, a_: - Infinity, b_:  Infinity] := 
    Module[\[IndentingNewLine]{s, t, n, k2, num, g, 
        dg2}, \[IndentingNewLine]g[
          x_] := \((1/
              Sqrt[2\ Pi])\)\ E^\((\(-\((x\^2)\)\)/
                2)\); \[IndentingNewLine]dg2 = 
        D[g[t], {t, 2}]; \[IndentingNewLine]k2 = 
        Integrate[t\^2\ f[t], {t, a, b}]; \[IndentingNewLine]num = 
        Integrate[f[t]\^2, {t, a, b}]; \[IndentingNewLine]den = 
        Integrate[dg2\^2, {t, \(-Infinity\), Infinity}]; \[IndentingNewLine]N[
        PowerExpand[
          k2\^\(\(-2\)/5\)\ \((num\/den)\)\^\(1/5\)]]]\)\(\[IndentingNewLine]\
\)
  \)\n
  \(KernelFunction[Biweight] = 
      If[Abs[#] < 1, N[\((15/16)\) \((1 - #^2)\)^2], 0] &;\)\n
  \(KernelFunction[Gauss] = N[E^\((\(-\((#^2)\)\)/2)\)/Sqrt[2\ Pi]] &;\)\n
  \(\(KernelFunction[f_] := f;\)\(\[IndentingNewLine]\)
  \)\n
  \(KernelDensity::scale = "\<Kernel smoothing constant = `1` s/n^(1/5)\>";\)\
\n
  \(\(KernelDensity::range = "\<Kernel width is `1` on interval [`2`,`3`]\>";\
\)\(\n\)
  \)\[IndentingNewLine]
  KernelDensity[expr_, opts___Rule] := 
    Module[\[IndentingNewLine]{h, scale, kernel, grid, kFunc, min, max, tab, 
        one, n, \[IndentingNewLine]data = expr}, \[IndentingNewLine]n = 
        Length[data]; \[IndentingNewLine]{h, grid, 
          kernel} = \({KernelScale, KernelGrid, KernelFunction} /. {opts}\) /. 
          Options[KernelDensity]; \[IndentingNewLine]kFunc = 
        KernelFunction[kernel]; \[IndentingNewLine]If[h === Automatic, 
        h = KernelScale[kernel]]; \[IndentingNewLine]Message[
        KernelDensity::scale, h]; \[IndentingNewLine]scale = 
        h\ N[Min[StandardDeviation[data], InterquartileRange[data]/1.345]/
              n^\((1/5)\)]; \[IndentingNewLine]If[\[IndentingNewLine]NumberQ[
          grid], \[IndentingNewLine]min = 
          Min[data] - 2\ scale; \[IndentingNewLine]max = 
          Max[data] + 2\ scale; \[IndentingNewLine]grid = 
          Range[min, 
            max, \((max - min)\)/
              grid], \[IndentingNewLine] (*else*) \[IndentingNewLine]min = 
          First[grid]; 
        max = Last[grid];\[IndentingNewLine]]; \[IndentingNewLine]Message[
        KernelDensity::range, scale, min, max]; \[IndentingNewLine]one = 
        Table[1, {n}]; \[IndentingNewLine]tab = 
        Table[\((one . \((kFunc /@ \((\((data - grid[\([i]\)])\)/
                          scale)\))\))\)/scale, {i, 1, Length[grid]}]/
          n; \[IndentingNewLine]Interpolation[
        Transpose[{grid, tab}]]\[IndentingNewLine]]\)



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.";

Options[Boxplot]={BoxVertical\[Rule]False,BoxWidth\[Rule]1,BoxCenter\[Rule]1,
      BoxFill\[Rule]True,AspectRatio\[Rule]0.5,Frame\[Rule]True};

Boxplot[x_List,opts___?OptionQ]:=Module[
    {box,vert,center,width,ar,frame,b,bi,fill,gOpts,
      ticks},{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_}\[RuleDelayed]{b,a},box[[2]]],
        Frame\[Rule]frame,AspectRatio\[Rule]ar,FrameTicks\[Rule]ticks,gOpts,
        PlotRange\[Rule]{{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]]]]}]




\!\(\(\( (*\ \(os\  = \ 
        Sqrt[2]\ Table[
            InverseErf[N[2\ \((i/\((n + 1)\))\) - 1]], {i, 
              n}];\)\ *) \)\(\[IndentingNewLine]\)\( (*\ \(os\  = \ 
        Quantile[NormalDistribution[], Range[n]\/\(n + 1\)];\)\ *) \)\)\)

\!\(QuantilePlot[x_List] := 
    Module[\[IndentingNewLine]{n = Length[x], \ \[IndentingNewLine]m\  = \ 
          Mean[x], \ \[IndentingNewLine]s\  = \ 
          StandardDeviation[x]}, \[IndentingNewLine]os\  = \@2\ Table[
            InverseErf[N[2\ \((i/\((n + 1)\))\) - 1]], {i, 
              n}]; \[IndentingNewLine]QuantilePlot[m\  + \ s\ os, 
        x]\[IndentingNewLine]]\)



ScatterPlot::usage="ScatterPlot[x_, y_] displays a scatter plot of y on x.";

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

ScatterPlotMatrix::usage=
    "ScatterplotMatrix[x_,...] draws a scatter plot matrix ofthe collection \
of input lists.";



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

Options[ScatterPlot]=
    Join[{ScatterPlotFunction\[Rule]None},Options[ListPlot]];

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

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

ScatterPlot[x_List,y_List,opts___?OptionQ]:=
    scat[ListPlot,Transpose[{x,y}],opts];

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

ScatterPlot3D[x_List,y_List,z_List,opts___?OptionQ]:=
  ListPlot3D[Transpose[{x,y,z}],opts]



Attributes[ScatterPlotMatrix] = HoldAll;

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

ScatterPlotMatrix[varSeq__]:=Module[
    {held=Hold[varSeq],
      vars={varSeq},
      i,j,k,c},
    k=Length[held];
    Show[GraphicsArray[
        Table[
          If[j>i,
            
            ListPlot[Transpose[{vars[[j]],vars[[i]]}],AxesLabel\[Rule]None,
              Axes\[Rule]None,DisplayFunction\[Rule]Identity,Frame\[Rule]True,
              FrameTicks\[Rule]None],(*else show correlation or name*)
       
                 If[j\[Equal]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}]]
      ]]



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.";

Options[Smooth]={SmoothSkipPercentage\[Rule].05,
      SmoothBinPercentage\[Rule].40};

Smooth[x_List,y_List,opts___?OptionQ]:=Module[
    {skip,pct,in,
      sx=x,
      sy=y,
      n=Length[x]},
    {skip,
        pct}={SmoothSkipPercentage,SmoothBinPercentage}/.{opts}/.Options[
          Smooth];
    If[Not[OrderedQ[x]],
      in=Ordering[x];sx=x\[LeftDoubleBracket]in\[RightDoubleBracket];
      sy=y\[LeftDoubleBracket]in\[RightDoubleBracket]];
    smthByPct[sx,sy,Ceiling[pct*n],Ceiling[skip*n]]
    ]

fit[x_,y_,at_]:=Module[
    {xBar=Mean[x],
      xDev,b},
    xDev=x-xBar;
    b=First[{xDev.y}/{xDev.xDev}];
    If[ListQ[at],
      Transpose[{at,Mean[y]+b (at-xBar)}],
      (*else*)
      {at,Mean[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+(Mean[x[[Range[binSize]]]]-min) Range[0,half-1]/half;
    xHi=max+(Mean[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,Mean[#1]]&,{Partition[x,binSize,skip],
            Partition[y,binSize,skip]}],
        fit[x[[n+1-Range[binSize]]],y[[n+1-Range[binSize]]],
          xHi]]
      ]]



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.";

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.";



Attributes[getNames] = {HoldFirst};

getNames[n_] :=Table[Extract[Hold[n],{1,i},HoldForm],{i,Length[n]}]



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\[Rule]True,Residuals\[Rule]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=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]




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\[Rule]BiweightWeightFunction,
        RobustConstant\[Rule]4.685,MaxIterations\[Rule]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=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,k,nextB,
        itCount=0,
        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\[Rule]convQ];
      If[itCount\[GreaterEqual]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\[Rule]est,
          Variance\[Rule]var,Scale\[Rule]s,
          WeightFunction\[Rule]weightFunction,FitResiduals\[Rule]e,
          PredictedResponse\[Rule]Y-e,RobustC\[Rule]c,Weights\[Rule]w},
        (*else lean and mean*)
        est]
      ];



On[General::spell1]