Back to library index.

Package regress (in regress.i) -

Index of documented functions or symbols:


DOCUMENT model = regress(y, [x1,x2,x3,...])
       or model = regress(y, [x1,x2,x3,...], axes, vals, chi2)

   Performs linear regression analysis, that is, a linear least squares
   fit of the parametric model
      y = model(1)*x1 + model(2)*x2 + model(3)*x3 + ...
        = model(+) * [x1,x2,x3,...](..,+)
   to the supplied data Y, [X1,X2,X3,...].  Y is an array, and the
   Xi arrays should have the same dimensions as Y, or more accurately
   [X1,X2,X3,...] must have one trailing dimension beyond the dimensions
   of Y; the length of the trailing dimension is the number of parameters
   in the model.

   Use the sigy= keyword to pass in the standard deviations of the Y
   values; sigy must be conformable with Y.  The returned model has
   minimum chi2, where
     chi2 = sum( (y - model(+) * [x1,x2,x3,...](..,+))^2 / sigy^2 )
   The default sigy=1, that is, all points have equal weight in chi2.

   The optional AXES, VALS, and CHI2 arguments are outputs.  AXES are
   the axes of the error ellipsoid of the fitted parameters, and VALS
   are the corresponding singular values.  AXES(i,) is the vector in
   the space of model() corresponding to VALS(i).  The VALS are
   non-negative and arranged in descending order.  CHI2 is the
   chi2 value for the returned model, divided by the number of
   degrees of freedom in the fit, numberof(y)-sum(vals>rcond*vals(1)).
   AXES(i,) are mutually orthogonal unit vectors; you should arrange
   the magnitudes of X and Y so this makes sense - their scale factors
   may matter.  Usually this means arranging that the magnitudes of
   the elements of the returned model not differ too wildly.

   Often and unexpectedly, the data do not permit a definitive choice
   of the model(); there may be entire subspaces of the model space
   which produce indistinguishably good fits to the data.  The signature
   of this problem is that some of the VALS may be zero or very small.
   The rcond= keyword permits you to specify how small a singular
   value, relative to the largest, VALS(1), you are willing to consider.
   The AXES(i,) corresponding to VALS(i) smaller than rcond*VALS(1)
   will be given zero contribution to the returned model().  The default
   value is rcond=1.e-9.  If any of the returned VALS is less than
   rcond*VALS(1) you should either change the scales of Y or some X,
   or remove some of the X (the ones which aren't contributing) to
   eliminate the problem.  The regress function reverses the sign of
   any VALS which are less than rcond*VALS(1), so you can quickly
   identify them.

     ab = regress(y, [1,x]);
       ab(1)+ab(2)*x   is best fit line to y(x)
     ab = regress(y, [1,x,x^2,x^3]);
       poly(x,ab(1),ab(2),ab(3))   is best fit cubic to y(x)
     ab = regress(y, [cos(x),sin(x)]);
       ab(1)*cos(x)+ab(2)*sin(x)   is best fit period 2*pi sine wave to y(x)

SEE ALSO: regress_cov


DOCUMENT cov = regress_cov(axes, vals)
       or cov = chi2 * regress_cov(axes, vals)

   Return the covariance matrix for the model returned by regress,
   where AXES and VALS are the values returned by regress.  This
   is a symmetric matrix representing
   that is, the variances (on the diagonal i=j) and correlations
   of the model parameters.

   If you did not specify the sigy= keyword in the original call to
   regress, then the absolute magnitudes of the covariance matrix
   elements do not mean much.  In case you wish to use the quality
   of the fit itself as an estimate of the errors sigy in the
   original Y values, assuming they are all equal, you can multiply
   by the CHI2 returned by regress, as in the second form above,
   in order to produce an estimated error in the fit.  See Numerical
   Recipes for caveats, but this is what commercial statistics
   software generally does.

   cov(1::numberof(vals)+1) are the variances of the model parameters

SEE ALSO: regress