regress.m package
- To: mathgroup at yoda.physics.unc.edu
- Subject: regress.m package
- From: Gossett <gossett at bethel.edu>
- Date: Wed, 5 Aug 1992 16:11:05 -0500
Jerry B. Keiper of Wolfram Research, Inc. has pointed out some weakness in
the data validation in the regress.m package I recently posted. I have
incorporated his suggestions into the revision that follows.
BeginPackage["Regression`"]
Regress::usage = "Regress[Data,Predictions,X_Plot_Range,Y_Plot_Range] does a linear regression on the data points. Data is either a list with the y values, or a list of {x,y} pairs. Predictions is an optional list of x values at which the least squares function is to be evaluated. X_Plot_Range and Y_Plot_Range are lists with an estimated lower and upper value for the axis limits. For example, Regress[{2,4,5,6,8,9},{3,7},{0,10},{0,10}] or
Regress[{{-1,3},{3,4},{4,4}},{-2,5},{0,6.25}]"
(* This package may be freely used for non-commercial purposes. It has been
designed for educational purposes, but may be useful in other contexts.
Copyright August 1992
Dr. Eric Gossett
Bethel College
3900 Bethel Drive
St. Paul, MN 55112
gossett at bethel.edu
(Version 1.1 fixed problems with vertical and horizontal data sets).
Version 1.2 (Improves the data validation, following suggestions from
Jerry B. Keiper at Wolfram Research, Inc. )
*)
Begin["`Private`"]
Unprotect[Regress]
ValidData[x_] := VectorQ[x, NumberQ] || MatrixQ[x, NumberQ];
Regress[Indata_?ValidData,Predictions_List:{},Domain_List,Range_List] :=
Module[{Data,L,Pmax,Pmin,NumData,Xvals,Yvals,CalcTable,
Xsum,Ysum,XYsum,X2sum,preds},
(* Finalize the data array: get length, add x coordinates if needed *)
NumData = Length[Indata];
If[Length[Dimensions[Indata]]==1,
Data = Table[{i,Indata[[i]]},{i,1,NumData}],
Data = Indata];
(* Check for vertical lines *)
Xvals = Transpose[Data][[1]];
Yvals = Transpose[Data][[2]];
If[Union[Flatten[Outer[(Chop[#1-#2])&,Xvals,Xvals]]] == {0},
Print["The data forms a vertical line"]; Return[]];
(* Compute the least squares line *)
L = Chop[Fit[Data,{1,x},x]];
(* Echo the original data *)
Print[" "];
Print["The Data set is:"];
Print[" "];
Print[TableForm[Data,TableAlignments->{Right,Right}]];
Print[" "];
(* Print the equation of the least squares line *)
If[SameQ[Head[L],Plus],
Print["The equation of the line is y = ",L[[2]][[1]]," x + ",
L[[1]] ],
If[FreeQ[L,x],Print["The equation of the line is y = ",L],
Print["The equation of the line is y = ",L[[1]]," x"]
]
];
Print[" "];
(* Plot the data and the least squares line *)
Pmax = Max[Map[(#[[2]])&,Data],
{L/.x->Domain[[1]],L/.x->Domain[[2]]},Range];
Pmin = Min[Map[(#[[2]])&,Data],
{L/.x->Domain[[1]],L/.x->Domain[[2]]},Range];
Show[ListPlot[Data,PlotRange->{Domain,{Pmin,Pmax}},
PlotStyle->{PointSize[.0125]},
DisplayFunction->Identity],
Plot[L,{x,Domain[[1]],Domain[[2]]},PlotRange->{Domain,{Pmin,Pmax}},
DisplayFunction->Identity],
DisplayFunction->$DisplayFunction];
(* Print any requested estimates using the least squares line *)
If[Length[Predictions] >0,
Print["The predicted values are:"];
Print[" "];
preds=Transpose[{Predictions,Table[L/.x->Predictions[[i]],
{i,1,Length[Predictions]}]}];
preds = Prepend[preds,{"X","Y"}];
Print[TableForm[preds,TableAlignments->{Right,Right}]];
];
(* Print the hand calculation table *)
Print[" "];
Print["The hand calculations are:"];
Print[" "];
CalcTable = Table[{Xvals[[n]],Yvals[[n]],Xvals[[n]] Yvals[[n]],
Xvals[[n]]^2},{n,1,NumData}];
Xsum = Apply[Plus,Xvals];
Ysum = Apply[Plus,Yvals];
XYsum = Apply[Plus,Part[Transpose[CalcTable],3]];
X2sum = Apply[Plus,Part[Transpose[CalcTable],4]];
CalcTable = Prepend[CalcTable,{"X","Y","XY","X^2"}];
CalcTable = Append[CalcTable,{"-----","-----","-----","-----"}];
CalcTable = Append[CalcTable,{Xsum,Ysum,XYsum,X2sum}];
CalcTable = Append[CalcTable,{Xsum^2," "," "," "}];
Print[TableForm[CalcTable,TableAlignments->{Right,Right}]];
] (* end Regress *)
End[]
Protect[Regress]
EndPackage[]