home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / STATISTI.PAK / LINEARRE.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  18.4 KB  |  508 lines
  1.  
  2.  
  3. (*:Version: Mathematica 2.0 *)
  4.  
  5. (*:Name: Statistics`LinearRegression *)
  6.  
  7. (*:Title: Linear Regression Analysis *)
  8.  
  9.  
  10. (*:Legal:
  11.   Copyright (c) 1990, Wolfram Research, Inc.
  12. *)
  13.  
  14. (*:Summary:
  15.   Least squares or weighted least squares linear regression
  16.   for data whose errors are assumed to be normally and 
  17.   independently distributed.
  18.   
  19. *)
  20.  
  21. (*:Keywords: linear regression, design matrix  *)
  22.  
  23. (*:Requirements: No special system requirements. *)
  24.  
  25. (*:Warning: None. *)
  26.  
  27. (*:Sources: Basic statistics texts. *)
  28.  
  29.  
  30. BeginPackage["Statistics`LinearRegression`",
  31.          "Statistics`Common`DistributionsCommon`",
  32.              "Statistics`NormalDistribution`",
  33.              "Statistics`ConfidenceIntervals`",
  34.              "Statistics`DescriptiveStatistics`"]
  35.  
  36. Regress::usage =
  37. "Regress[data, funs, vars] finds a least-squares fit to a list of data as
  38. a linear combination of the functions funs of variables vars. The data
  39. can have the form {{x1, y1, ..., f1}, {x2, y2, ..., f2}, ...}, where the
  40. number of coordinates x, y, ... is equal to the number of variables in
  41. the list vars.  The data can also be of the form {f1, f2, ...}, with a
  42. single coordinate assumed to take values 1, 2, .... The argument funs can
  43. be any list of functions that depend only on the objects vars."
  44.  
  45. Options[Regress] = {OutputList -> Null, OutputControl->Automatic,
  46.    IncludeConstant -> True, BasisNames -> Automatic, 
  47.    Weights -> Equal}  ~Join~   Options[SingularValues] ~Join~ 
  48.    Options[StudentTCI]
  49.  
  50. DesignedRegress::usage =
  51. "DesignedRegress[designmatrix,rdata] finds a least-squares fit given a 
  52. design matrix and a vector of the response data. "
  53.  
  54. Options[DesignedRegress] = { OutputList -> Null, Weights -> Equal,
  55.     BasisNames -> Automatic, OutputControl->Automatic } ~Join~
  56.     Options[SingularValues] ~Join~ Options[StudentTCI]
  57.  
  58.  
  59. DesignMatrix::usage = 
  60. "DesignMatrix[data,funs,vars] gives the design matrix for modeling the
  61. data as a linear combination of the functions. The data can have the 
  62. form {{x1, y1, ..., f1}, {x2, y2, ..., f2}, ...}, where the
  63. number of coordinates x, y, ... is equal to the number of variables in
  64. the list vars.  The data can also be of the form {f1, f2, ...}, with a
  65. single coordinate assumed to take values 1, 2, .... The argument funs can
  66. be any list of functions that depend only on the objects vars."
  67.  
  68.  
  69.         (*  options  *)
  70.  
  71. BasisNames::usage = 
  72. "BasisNames is an option to regression functions and is used to 
  73. include headings for each of the basis functions in the output."
  74.  
  75. IncludeConstant::usage =
  76. "IncludeConstant is an option to the Regress function and is
  77. used to indicate whether a constant term is automatically 
  78. included in the models."
  79.  
  80. OutputList::usage =
  81. "OutputList is an option to regression functions and is used to 
  82. specify a list of statistics to be reported about the fit."
  83.  
  84. OutputControl::usage = 
  85. "OutputControl is an option to regression functions and is used to 
  86. specify whether standard information is automatically included
  87. in the output."
  88.  
  89. Weights::usage = 
  90. "Weights is an option to regression functions and is used to define a
  91. vector of weights for the variances.  When weights are specified, the
  92. weighted least-squares fit is calculated."
  93.  
  94. NoPrint::usage =
  95. "NoPrint specifies that the OutputControl option should force only the
  96. information specified by the OutputList option to be returned."
  97.  
  98.     (*  Possible list elements for OutputList  *)
  99.  
  100.  
  101. BestFit::usage =
  102. "BestFit is used in the output of regression functions to
  103. identify the best fit."
  104.  
  105. BestFitCoefficients::usage =
  106. "BestFitCoefficients is used in the output of regression
  107. functions to identify the list of values of the fitted parameters
  108. in the best fit."
  109.  
  110. ANOVATable::usage =
  111. "ANOVATable is used in the output of regression functions to
  112. identify the analysis of variance table.  This table includes F-ratio
  113. tests for comparing the models specified in the function, as well as
  114. the degrees of freedom and sums of squares used to obtain these
  115. F ratios."
  116.  
  117. ConfidenceIntervalTable::usage =
  118. "ConfidenceIntervalTable is used in the output of regression
  119. functions to identify a table of confidence intervals for the best
  120. fit parameters."
  121.  
  122. CovarianceMatrix::usage =
  123. "CovarianceMatrix is used in the output of regression functions
  124. to identify the estimated covariance matrix of the fit coefficients."
  125.  
  126. CorrelationMatrix::usage =
  127. "CorrelationMatrix is used in the output of regression functions
  128. to identify the estimated correlation matrix of the fit coefficients."
  129.  
  130. EstimatedVariance::usage =
  131. "EstimatedVariance is used in the output of regression
  132. functions to identify the variance estimate, or the residual 
  133. mean square."
  134.  
  135. FitResiduals::usage =
  136. "FitResiduals is used in the output of regression functions
  137. to identify the list of differences between the response data and
  138. the best fit evaluated at the same abscissa points."
  139.  
  140. ParameterTable::usage =
  141. "ParameterTable is used in the output of regression functions 
  142. to identify a table of information about the best fit parameters."
  143.  
  144. PredictedResponse::usage =
  145. "PredictedResponse is used in the output of regression functions
  146. to identify the best fit evaluated at the data points."
  147.  
  148. RSquared::usage =
  149. "RSquared is used in the output of regression functions to
  150. identify the square of the multiple correlation coefficient."
  151.  
  152. AdjustedRSquared::usage =
  153. "AdjustedRSquared is used in the output of regression functions
  154. to identify the multiple correlation coefficient adjusted for the number
  155. of degrees of freedom in the fit."
  156.  
  157.  
  158.  
  159.         Begin["`private`"]
  160.  
  161. Regress[arguments___] :=
  162.     Block[{answer = vRegress[arguments]},
  163.         answer /; answer =!= Fail
  164.     ]
  165.  
  166. vRegress[data_, basis_?VectorQ, vars_, options___] :=
  167.     Module[
  168.      {varlist, rdata, cbasis, cpresent,report,localDesignMatrix,
  169.     name,optionlist,old},
  170.     optionlist = {options};
  171.         If[Head[vars] === List, varlist = vars, varlist = {vars}];
  172.     {idata, rdata} = GetData[data, Length[varlist]];
  173.         If[idata === Fail, Return[Fail]];
  174.            cbasis = basis;
  175.     cpresent = Or @@ Apply[And, Outer[FreeQ, basis, varlist], {1}];
  176.     If[And[!cpresent, IncludeConstant /.optionlist/. Options[Regress]],
  177.             PrependTo[cbasis, 1];
  178.             cpresent = True
  179.         ];
  180.     localDesignMatrix = DesignMatrix[idata, cbasis, varlist];
  181.     If[Length[rdata] <= Length[cbasis],
  182.             Message[Regress::mindata]
  183.         ];
  184.     name = BasisNames/.optionlist/.Options[Regress];
  185.     If[name == Automatic,
  186.         optionlist=ChangeValue[optionlist,BasisNames,cbasis]];
  187.         report = DesignedRegress[localDesignMatrix, rdata, optionlist];
  188.     If[!FreeQ[report, BestFit],
  189.       old = GetValue[report,BestFit];
  190.           report = ChangeValue[report,BestFit,old.cbasis]];
  191.     report
  192.     ]
  193.  
  194. vRegress[badargs___]:=
  195.     (If[badargs =!= Fail,Message[Regress::regf]];Fail)
  196.  
  197. DesignMatrix[arguments___] :=
  198.     Block[{answer = vDesignMatrix[arguments]},
  199.         answer /; answer =!= Fail
  200.     ]
  201.  
  202. vDesignMatrix[data_, basis_?VectorQ, vars_] :=
  203.     Module[{PureFunction,varlist,ndata,func,slots,rules,n},
  204.     If[Head[vars] === List, varlist = vars, varlist = {vars}];
  205.     If[Length[data[[1]]]==0,
  206.         ndata = Transpose[{Range[Length[data]],data}],
  207.         ndata = data];
  208.     If[!ListQ[data],Message[DesignMatrix::desd]; Return[Fail]];
  209.     PureFunction = Map[Release,PrepareFunction[basis,varlist]];
  210.     If[!FreeQ[basis,Slot],
  211.            Message[Regress::slot];
  212.        slots = Table[Slot[n], {n, 1, Length[varlist]}];
  213.        rules = Thread[slots -> varlist];
  214.        PureFunction = PureFunction /.rules ]  ;      
  215.     Apply[PureFunction,ndata, 1]
  216.     ]
  217.  
  218. vDesignMatrix[badargs___]:=
  219.     (If[badargs =!= Fail,Message[DesignMatrix::desf]];Fail)
  220.  
  221.  
  222. PrepareFunction[basis_, vars_] := 
  223.      If[VectorQ[vars],Return[Function[vars, basis]] ]
  224.  
  225. DesignedRegress[arguments___] :=
  226.     Block[{answer = vDesignedRegress[arguments]},
  227.         answer /; answer =!= Fail
  228.     ]
  229.  
  230. vDesignedRegress[designmatrix_List, rdata_?VectorQ, options___] :=
  231.     (* designmatrix::design matrix for regression                       *)
  232.     (* rdata::vector of response values                                 *)
  233.     (* optionlist::list of options                                      *)
  234.     Module[{names, optionlist, weights, sqrtweights, beta,
  235.       ValidWeights = False,
  236.           localRData = rdata, localDesignMatrix = designmatrix,
  237.           tol, SVD, u, w, v, (* input and output of SingularValues *)
  238.           OrthoResponse, (* orthogonal projection of localRData *)
  239.           fitlist, (* best fit coefficients *)
  240.           localPredictedResponse, localFitResiduals,
  241.           ModelDOF, TotalDOF, ErrorDOF,
  242.           ModelSS, TotalSS, ErrorSS,
  243.           TotalResiduals, (* residuals from fitting a constant *)
  244.           localRSquared, 
  245.           ResponseVariance, (* estimated variance from best fit *)
  246.           localPValue, localFRatio, CM, rootCM, 
  247.           report, (* output accumulation *)
  248.           SEList, (* list of standard errors *)
  249.           CILevel, CIList, ctable,
  250.           localCorrelationMatrix,
  251.           TRatioList, tpvalues, ttable        },
  252.  
  253.      (* READ IN OPTIONS *)
  254.  
  255.            optionlist = {options};
  256.     If[Length[optionlist] == 1 && !VectorQ[optionlist],
  257.         {optionlist} = optionlist];
  258.         tol = Tolerance /. optionlist /. Options[DesignedRegress];
  259.     weights = Weights /. optionlist /. Options[DesignedRegress];
  260.     CILevel = ConfidenceLevel /. optionlist /. Options[DesignedRegress];
  261.         ValidWeights = VectorQ[weights] && Length[weights]==Length[rdata] &&
  262.         Apply[And, Map[NumberQ, N[weights]]];
  263.     If[ValidWeights,
  264.             sqrtweights = Sqrt[weights];
  265.         localDesignMatrix = sqrtweights localDesignMatrix;
  266.         localRData = sqrtweights localRData,
  267.         If[weights =!= Equal, Message[Regress::regw, Length[rdata]]]
  268.     ];
  269.  
  270.      (* NUMERICALIZE DESIGN MATRIX & RESPONSE DATA *)
  271.         {localDesignMatrix, localRData} =
  272.                 MakeNumerical[localDesignMatrix, localRData];
  273.         If[localDesignMatrix === Fail, Return[Fail]];
  274.  
  275.     (*  GET SVD AND COMPUTE COEFFICIENTS *) 
  276.  
  277.         SVD = SingularValues[localDesignMatrix, Tolerance -> tol];
  278.         If[Head[SVD] === List && Length[SVD] === 3,
  279.             {u, w, v} = SVD,
  280.         (* else *)
  281.             Message[Regress::noSVD];
  282.             Return[Fail]
  283.         ];
  284.         If[Length[w] < Length[localDesignMatrix[[1]]], 
  285.         Message[Regress::notdep]];  
  286.  
  287.         OrthoResponse = u . localRData;
  288.  
  289.         fitlist = Transpose[v] . (OrthoResponse / w);   
  290.  
  291.         localPredictedResponse = localDesignMatrix . fitlist;
  292.         localFitResiduals = localRData - localPredictedResponse;
  293.         ErrorDOF = Length[rdata] - Length[Transpose[localDesignMatrix]];
  294.         ErrorSS = localFitResiduals . localFitResiduals;
  295.  
  296.      (*  COMPARING MODELS  *)
  297.  
  298.         If[IncludeConstant /.optionlist/. Options[Regress],
  299.       beta = If[ValidWeights,
  300.             (* numericalize sqrtweights *)
  301.             If[Precision[sqrtweights] == Infinity,
  302.                 sqrtweights = N[sqrtweights]];
  303.             sqrtweights (localRData.sqrtweights)/Apply[Plus,weights],
  304.             Mean[localRData]];
  305.       TotalResiduals = localRData -  beta;
  306.           TotalDOF = Length[rdata] - 1,
  307.         (* else *)
  308.           TotalResiduals = localRData;
  309.           TotalDOF = Length[rdata],
  310.         (* error *)
  311.           Message[Regress::uncon];
  312.           Return[Fail]
  313.         ];
  314.         TotalSS = TotalResiduals . TotalResiduals;
  315.         ModelDOF = TotalDOF - ErrorDOF;
  316.         ModelSS = TotalSS - ErrorSS;
  317.  
  318.         ResponseVariance = ErrorSS/ErrorDOF;
  319.     localRSquared = ModelSS/TotalSS;
  320.  
  321.  
  322.     (*  REPORT GENERATION *)
  323.  
  324.         report = GetOutputList[optionlist];
  325.      If[MemberQ[report, ANOVATable],
  326.             localFRatio = (ModelSS/ModelDOF)/ResponseVariance;
  327.             localPValue = Chop[1-CDF[FRatioDistribution[ModelDOF, 
  328.         ErrorDOF],localFRatio],.00001];
  329.             report = report /. ANOVATable ->
  330.                 (ANOVATable ->
  331.                 TableForm[{{ModelDOF, ModelSS, ModelSS/ModelDOF,
  332.                                         localFRatio, localPValue},
  333.                        {ErrorDOF, ErrorSS, ErrorSS/ErrorDOF},
  334.                        {TotalDOF, TotalSS}},
  335.             TableHeadings->{{"Model", "Error", "Total"},
  336.                {"DoF", "SoS", "MeanSS", "FRatio", "PValue"}}
  337.                           ]  );
  338.         ];
  339.   
  340.         (* compute things involving the covariance matrix *)
  341.         If[Or @@ Map[MemberQ[report, #]&, {CovarianceMatrix,
  342.                 CorrelationMatrix, ParameterTable,
  343.                 ConfidenceIntervalTable}],
  344.         names = BasisNames/.optionlist/.Options[Regress];
  345.             CM = Transpose[v] . (v / w^2); 
  346.             rootCM = Sqrt[DiagonalElements[CM]]; 
  347.             SEList = Sqrt[ResponseVariance] rootCM;  
  348.             CIList := Apply[N[StudentTCI[#1, #2, ErrorDOF,
  349.                                    ConfidenceLevel -> CILevel]]&,
  350.                            Transpose[{fitlist, SEList}], 1];        
  351.             localCorrelationMatrix := Transpose[Transpose[CM/rootCM]/rootCM];
  352.             TRatioList = fitlist / SEList;
  353.             tpvalues = Map[
  354.                      CDF[StudentTDistribution[ErrorDOF],#]&,
  355.                     TRatioList];
  356.         tpvalues =Chop[2 Map[If[#>.5, 1-#, #]&,tpvalues],.00001];
  357.             ttable := Transpose[{fitlist, SEList, TRatioList, tpvalues}];
  358.         ctable := Transpose[{fitlist, CIList}];
  359.         report = report /. CovarianceMatrix ->
  360.                               (CovarianceMatrix -> MatrixForm[CM]);
  361.          report = report /. CorrelationMatrix ->
  362.                               (CorrelationMatrix ->
  363.                                    MatrixForm[localCorrelationMatrix]);
  364.         report = report /. ParameterTable ->
  365.              (ParameterTable -> TableForm[ttable,
  366.               TableHeadings -> {names, {"Estimate", "SE", "TStat", "PValue"}}]);
  367.             report = report /. ConfidenceIntervalTable ->
  368.              (ConfidenceIntervalTable -> TableForm[ctable, TableDepth -> 2,
  369.           TableHeadings -> {names, {"Estimate", "CI"}}])
  370.         ];
  371.         report = report /. RSquared -> (RSquared -> localRSquared);
  372.         report = report /. AdjustedRSquared ->
  373.             (AdjustedRSquared -> 1-(1-localRSquared)(TotalDOF/ErrorDOF));
  374.     report = report /. EstimatedVariance ->
  375.                 (EstimatedVariance -> ResponseVariance);
  376.     report = report /. BestFit -> (BestFit->fitlist);
  377.         report = report /. BestFitCoefficients ->
  378.                 (BestFitCoefficients -> fitlist);
  379.           (* correct the predicted response and residuals if weighted.  *)
  380.     If[ValidWeights,
  381.          localPredictedResponse /= sqrtweights;
  382.          localFitResiduals /= sqrtweights ];
  383.         report = report /. PredictedResponse :>
  384.                 (PredictedResponse -> localPredictedResponse);
  385.         report = report /. FitResiduals :>
  386.                 (FitResiduals -> localFitResiduals);
  387.         report
  388.     ] (* end of DesignedRegress *)
  389.  
  390. vDesignedRegress[badargs___]:=
  391.     (If[badargs =!= Fail,Message[Regress::regd]];Fail)
  392.  
  393. DiagonalElements[matrix_] :=
  394.     Module[{k, size = Length[matrix]},
  395.         Table[matrix[[k, k]], {k, 1, size}]
  396.     ]
  397.  
  398.  
  399. GetValue[list_,name_]:=Module[{rule},
  400.     rule = Select[list,!FreeQ[#,name]&];
  401.     name /. rule]
  402.  
  403. ChangeValue[list_,name_,newval_]:=
  404.     If[FreeQ[list,name],Prepend[list,(name->newval)],
  405.     Map[If[!FreeQ[#,name],
  406.     Replace[name->_,(name->_)->(name->newval)],#]&,list]]
  407.  
  408.  
  409. GetData[data_?VectorQ, varcount_] :=
  410.     If[varcount === 1,
  411.         {Transpose[{Range[Length[data]],data}], data},
  412.         Message[Regress::fitc, 1, varcount]; {Fail, Fail}
  413.     ]
  414.  
  415. GetData[data_?MatrixQ, varcount_] :=
  416.     If[Length[data[[1]]] > varcount,
  417.         {data, Map[Last, data]},
  418.         Message[Regress::fitc, Length[data[[1]]]-1, varcount];
  419.         {Fail, Fail}
  420.     ]
  421.  
  422. GetData[___] := (Message[Regress::regd]; {Fail, Fail})
  423.  
  424. MakeNumerical[DesignMatrix_, rdata_] :=
  425.     (* The purpose of MakeNumerical is to determine if DesignMatrix     *)
  426.     (* is in a form appropriate for SingularValues.                  *)
  427.   Module[{dmPrecision = Precision[{DesignMatrix, rdata}],
  428.           nDesignMatrix, nrdata},
  429.     {nrdata, nDesignMatrix} =
  430.        If[dmPrecision < Infinity,
  431.           N[{rdata, DesignMatrix}, dmPrecision],
  432.           N[{rdata, DesignMatrix}]
  433.        ];
  434.        If[MatrixQ[nDesignMatrix],
  435.           {nDesignMatrix, nrdata},
  436.         (* else *)
  437.           Message[Regress::notnum];
  438.           {Fail, Fail}
  439.        ]
  440.   ]
  441.  
  442. GetOutputList[options_] :=
  443.     Module[{localReport,control,outlist},
  444.           control = OutputControl /. options /.Options[Regress];
  445.     outlist = OutputList /. options /.Options[Regress];
  446.      If[outlist === Null, localReport = {}, localReport = outlist];
  447.      Which[control === Automatic,
  448.       localReport = Join[localReport,{ParameterTable,RSquared,
  449.         AdjustedRSquared,EstimatedVariance,ANOVATable}],
  450.        control =!= NoPrint,
  451.         Message[Regress::optx, control, OutputControl -> control]
  452.     ];
  453.        localReport           ]
  454.  
  455.  
  456. DesignMatrix::desf = "The second argument to DesignMatrix is not a list 
  457. of functions."
  458.  
  459. DesignMatrix::desd = "The first argument to DesignMatrix is not a list 
  460. or rectangular array of data."
  461.  
  462. Regress::regf = "The second argument to Regress is not a list of 
  463. functions."
  464.  
  465. Regress::optx = "Unknown option `1` in `2`."
  466.  
  467. Regress::uncon = "Regress was unable to determine that a constant term is 
  468. present in the model. (This should not happen under normal conditions.)"
  469.  
  470. Regress::notnum = "The basis functions yielded non-numerical results
  471. when evaluated at the data points."
  472.  
  473. Regress::regc = "Number of independent coordinates (`1`) does not
  474. equal the number of parameters (`2`)."
  475.  
  476. Regress::slot = "Warning: The basis functions contain Slot specifiers
  477. that will be interpreted as the corresponding variables from the
  478. variable list."
  479.  
  480. Regress::regd = "The first argument to Regress is not a list or 
  481. rectangular array of data."
  482.  
  483. Regress::regw = "Warning: the Weights option must be set to a vector of numbers
  484. having length ``, the length of the first argument to Regress.  Weights ->
  485. Equal assumed."
  486.  
  487. Regress::notdep = "The fit is numerically independent of a linear
  488. combination of the basis functions."
  489.  
  490. Regress::noSVD = "Regress was unable to obtain the singular value
  491. decomposition for the design matrix of this problem."
  492.  
  493. Regress::mindata = "There are as many or more basis functions than data
  494. points.  Subsequent results may be misleading."
  495.  
  496.         End[]
  497.  
  498.  
  499. SetAttributes[ Regress ,ReadProtected];
  500. SetAttributes[ DesignedRegress ,ReadProtected];
  501. SetAttributes[ DesignMatrix ,ReadProtected];
  502.  
  503. Protect[ Regress, DesignedRegress, DesignMatrix];
  504.  
  505.  
  506.         EndPackage[]
  507.  
  508.