Print["Updated: Dec 8, 2000"]; 
(* Author: A.I. McLeod *)
(* Copyright (c), A.I. McLeod, 2000 *)
(* Last update: December 8, 2000 *)
(* bilo *)
(* Summary:
Bivariate Loess
*)

BeginPackage["BILO`"]

Needs["Statistics`ContinuousDistributions`"];

SLPlot::usage=
"SLPlot[fits, res, p:0.5, \[Alpha]:1] plots the improved
Spread-Location plot given vectors of fits and res.
The power transformation applied to the absolute residuals
is specified by p and defaults to 0.5.  The span of the
robust loess local linear smoother is \[Alpha] and defaults
to 1."

RobustLoessPlot::usage=
"RobustLoessPlot[data, \[Alpha]:0.6, \[Lambda]:1, opts]
plots the data and the
robust loess curve to data= {{x\_1, y\_1}, {x\_2, y\_2},  ... , 
{x\_n, y\_n}}.  The\ span \[Alpha]\ is\  >  0. 
The\ degree\ of\ the\ local\ linear\ polynomial\ is\ specified\ 
by\ \[Lambda]\ and\ it\ is\ assumed\ that\ \[Lambda]\  
\[GreaterEqual] \ 
0.  The default span and degree are 0.6 and 1 respectively.
RobustLoessPlot can be given additional standard graphics
options as well as the option MaxIterations.  MaxIterations
puts an upper limit (default is 25) on the number of iterations
in the IRWLS algorithm.  If this upper limit is reached a warning
message will be given."

RobustLoessFit::usage=
"RobustLoessFit[data, \[Alpha]:0.75, \[Lambda]:1, opts] fits a
robust loess curve to data= {{x\_1, y\_1}, {x\_2, y\_2},  ... , 
{x\_n, y\_n}}.  The\ span \[Alpha]\ is\  >  0. 
The\ degree\ of\ the\ local\ linear\ polynomial\ is\ specified\ 
by\ \[Lambda]\ and\ it\ is\ assumed\ that\ \[Lambda]\  
\[GreaterEqual] \ 
0."

LoessFit::usage = 
"LoessFit[x, data, \[Alpha], \[Lambda]] 
fits\ a\ loess\ curve\ at\ 
the\ point\ x\ for\ 
data = {{x\_1, y\_1}, {x\_2, y\_2},  ... , {x\_n, y\_n}} . 
The\ span
\[Alpha]\ is\  >  0. 
The\ degree\ of\ the\ local\ linear\ polynomial\ is\ specified\ 
by\ \[Lambda]\ and\ it\ is\ assumed\ that\ \[Lambda]\  
\[GreaterEqual] \ 
0. The\ point\ x\ may\ be\ a\ vector\ in\ which\ case\ the\ loess\ 
curve\ is\ evaluated\ for\ each\ element\ of\ the\ vector . \ The\ 
default\ values\ of\ \[Alpha]\ and\ \[Lambda]\ are\ 0.75\ and\ 1";

LoessPlot::usage = \
    "LoessPlot[data,\[Alpha],\[Lambda],numvalues] 
plots the data and the loess curve.";

BWPlot::usage=
	"BWPlot[data] box-and-whisker plot of data vector or data matrix.
In the case of a matrix the columns are taken as variables and the
plot is produced horizontally so the distributions of the variables
are most easily compared. The box goes from Q1 to Q3. The median
is shown as a point between Q1 and Q3. The whiskers extend out to
adjacent values which are the most extreme values in (Q1-r,Q3+r),
r=1.5(Q3-Q1). Data values outside this interval are plotted. Under
normality the probability of an outside value is about 0.7%."

(*==========================================================================*)
Begin["`Private`"]

DataRange[x_] := {Min[x], Max[x]}
 
LoessFit[(x_)?NumberQ, data_, \[Alpha]_:0.75, \[Lambda]_:1] := 
  WLSFit[data, LoessWts[x, data, \[Alpha]], \[Lambda], x]
 
LoessFit[(x_)?VectorQ, data_, \[Alpha]_:0.75, \[Lambda]_:1] := 
  Table[LoessFit[x[[i]], data, \[Alpha], \[Lambda]], {i, Length[x]}]
 
SLPlot[fits_, res_, p_:0.5, \[Alpha]_:1] := 
  Module[{a, f, r, data2, s, lines}, 
   data2 = Sort[Transpose[{fits, res}], First[#1] < First[#2] & ]; 
    {f, r} = Transpose[data2]; a = If[p == 0, Log[Abs[r]], Abs[r]^p]; 
    data2 = Transpose[{f, a}]; s = RobustLoessFit[data2, \[Alpha]]; 
    lines = Line[Transpose[{f, s}]]; 
    ListPlot[data2, PlotRange -> All, Axes -> False, Frame -> True, 
     FrameLabel -> {"fit", "Abs[res]^p"}, 
     PlotStyle -> {PointSize[0.02], RGBColor[0, 0, 1]}, 
     Epilog -> {RGBColor[0, 1, 0], Thickness[0.02], lines}]; BWPlot[a - s]]

RobustLoessFit::MaximumReached=
"The maximum number of iterations of the IRWLS algorithm, as 
specified by the option MaxIterations, has been reached and without
convergence of the algorithm.
You could try increasing MaxIterations.";

RobustLoessFit[data_, \[Alpha]_:0.75, \[Lambda]_:1, (opts___)?OptionQ] := 
  Module[{x, y, \[Delta], res, rsum = 0, data2, iter = 0, rprev = 1, 
    r = Table[1, {Length[data]}]}, 
   maxiter = MaxIterations /. {opts} /. Options[RobustLoessFit]; 
    data2 = Sort[data, First[#1] < First[#2] & ]; {x, y} = Transpose[data2]; 
    \[Delta] = Table[LoessWts[x[[i]], data, \[Alpha]], {i, 1, Length[x]}]; 
    While[++iter < maxiter && Abs[rsum - rprev] > 0.001, 
     res = Table[y[[i]] - WLSFit[data, \[Delta][[i]]*r, \[Lambda], x[[i]]], 
        {i, 1, Length[x]}]; r = BiSquare[res/(6*Median[Abs[res]])];
      rprev = rsum; rsum = Abs[Plus @@ r; ]]; 
If[iter==maxiter,Message[RobustLoessFit::MaximumReached]];
y - res]
 
Options[RobustLoessFit] = {MaxIterations -> 25}
 
LoessSummary[data_, \[Alpha]_:0.75, \[Lambda]_:1] := 
  With[{res = LoessResiduals[data, \[Alpha], \[Lambda]]}, 
   {res -> res, \[Sigma] -> Sqrt[Plus @@ (res^2)/Length[res]], 
    \[Mu] -> 1.199999999999999*(\[Lambda] + 1)/\[Alpha]}]
 
LoessResiduals[data_, \[Alpha]_:0.75, \[Lambda]_:1] := 
  Last[Transpose[data]] - LoessFit[First[Transpose[data]], data, \[Alpha], 
    \[Lambda]]
 
RobustLoessPlot[data_, \[Alpha]_:0.6, \[Lambda]_:1, opts___] := 
  Module[{x, y, fits, data2, lines}, 
   data2 = Sort[data, First[#1] < First[#2] & ]; {x, y} = Transpose[data2]; 
    fits = RobustLoessFit[data2, \[Alpha], \[Lambda]]; 
    lines = Line[Transpose[{x, fits}]]; 
    ListPlot[data2, PlotRange -> All, Axes -> False, Frame -> True, 
     PlotStyle -> {PointSize[0.02], RGBColor[0, 0, 1]}, 
     Epilog -> {RGBColor[0, 1, 0], Thickness[0.02], lines},opts]]
 
LoessPlot[data_, \[Alpha]_:0.75, \[Lambda]_:1, numvalues_:30, opts___] := 
  With[{x = EquispaceVector[First[Transpose[data]], numvalues]}, 
   ListPlot[data, opts, PlotStyle -> {PointSize[0.05]}, 
    PlotRange -> ScaleRectangle[data], Frame -> True, Axes -> False, 
    Epilog -> {Thickness[0.02], RGBColor[0, 1, 1], 
      Line[Transpose[{x, LoessFit[x, data, \[Alpha], \[Lambda]]}]]},opts]]
 
WLSFit[data_, wts_, ldegree_:1, x_] := 
  Fit[Transpose[(wts*#1 & ) /@ 
      Join[{Table[1, {Length[data]}]}, Transpose[data]]], 
    Join[{u}, Table[v^i, {i, ldegree}]], {u, v}] /. {u -> 1, v -> x}
 
LoessWts[x_, data_, \[Alpha]_] := 
  Tricube[(x - First[Transpose[data]])/LoessDistance[x, data, \[Alpha]]]
 
Tricube = Compile[{{x, _Real, 1}},If[Abs[#]<1, (1-#^2)^2, 0]& /@x];

Bisquare = Compile[{{x, _Real, 1}},
            If[Abs[#]<1, (1-Abs[#]^3)^3, 0]& /@x];
 
LoessDistance[x_, data_, \[Alpha]_] := 
  Module[{A = Max[1, \[Alpha]], X = First[Transpose[data]], q}, 
   q = Min[Length[X], Ceiling[\[Alpha]*Length[X]]]; A*Sort[Abs[X - x]][[q]]]
 
EquispaceVector[(x_)?VectorQ, numvalues_:30] := 
  Range[Min[x], Max[x], N[(Max[x] - Min[x])/numvalues]]
 
ScaleRectangle[data_] := 
  With[{x = Transpose[data]}, (AddEps[#1] & ) /@ 
    {DataRange[First[x]], DataRange[Last[x]]}]
 
AddEps[{xlo_, xhi_}] := 
  With[{\[Epsilon] = (xhi - xlo)*0.05}, {xlo - \[Epsilon], xhi + \[Epsilon]}]
 
BWPlot[data_] := 
  Module[{datadim, k, datapts, whiskers, box, outsidepts, medpt, dmax, dmin, 
    coldata, Q1, Q3, uplim, dnlim, outside, jitter, drange, epsilon, 
    boxwidth = 0.4}, datadim = Dimensions[data]; 
    k = If[Length[datadim] == 2, datadim[[2]], 1]; 
    datapts = outsidepts = whiskers = box = medpt = dmax = dmin = {}; 
    Do[coldata = If[k == 1, data, Column[data, i]]; 
      datapts = 
       Join[datapts, Transpose[{coldata, Table[i, {Length[coldata]}]}]]; 
      medpt = Join[medpt, {PointSize[0.04], 
         Point[{Quantile[coldata, 0.5], i}]}]; 
      Q1 = Quantile[coldata, 0.25]; Q3 = Quantile[coldata, 0.75]; 
      box = Join[box, {RGBColor[0.690207, 0.7685929999999999, 0.870602], 
         Polygon[{{Q1, i - boxwidth}, {Q1, i + boxwidth}, 
           {Q3, i + boxwidth}, {Q3, i - boxwidth}}]}]; 
      step = 1.5*(Q3 - Q1); uplim = Q3 + step; dnlim = Q1 - step; 
      upadj = Max[Select[coldata, #1 <= uplim & ]]; 
      dnadj = Min[Select[coldata, #1 >= dnlim & ]]; 
      whiskers = 
       Join[whiskers, {Thickness[0.005], Line[{{dnadj, i}, {Q1, i}}], 
         Line[{{Q3, i}, {upadj, i}}]}]; 
      outside = 
       Join[Select[coldata, #1 > uplim & ], Select[coldata, #1 < dnlim & ]]\
; jitter = Table[i + (Random[] - 1/2)/6, {Length[outside]}]; 
      outsidepts = 
       Join[outsidepts, (Circle[#1, Offset[{5, 5}]] & ) /@ 
         Transpose[{outside, jitter}]]; dmax = Max[dmax, Max[coldata]]; 
      (dmin = Min[dmin, Min[coldata]]; ), {i, k}]; 
    epsilon = (dmax - dmin)*0.03; drange = {dmin - epsilon, dmax + epsilon}; 
    ListPlot[datapts, Ticks -> {Automatic, None}, 
     PlotStyle -> {AbsolutePointSize[0]}, PlotRange -> {drange, {0, k + 1}}, 
     Axes -> {True, False}, Epilog -> {outsidepts, box, medpt, whiskers}]]
 
  End[];

EndPackage[]

