File Coverage

blib/lib/Statistics/OLS.pm

Criterion	Covered	Total	%
statement	131	198	66.1
branch	25	60	41.6
condition	1	3	33.3
subroutine	16	18	88.8
pod	8	15	53.3
total	181	294	61.5

#================---------=================================================#
# Statistics::OLS -- Perform Ordinary Least Squares (with statistics) 2-D  #
#                     by Sanford Morton                #
#==========================================================================#

# Revision history for Perl module Statistics::OLS.
#
# 0.01 - 22 March 1998
# 	   - original version
#
# 0.02 - 29 March 1998
# 	   - corrected array bounds check bug in setData
# 	   - included check for divide by zero in standard error of
# 	     coefficients and t-stats
#
# 0.03 - 31 May 1998
# 	   - placed module into standard format using h2xs
#
# 0.04 - 13 July 1998
# 	   - changed the name from Statistics::Ols to Statistics::OLS
#
# 0.05 - 15 Sep 1999
# 	   - corrected error checking bug
# 	   - corrected pod documentation bug
#
# 0.06 - 4 July 2000
# 	   - allowed data in scientific (exponential) notation
#
# 0.07 - 12 October 2000
#          - _sse fix for potential precision problems

package Statistics::OLS;

$Statistics::OLS::VERSION = '0.07';

use strict;

#==================#
#  public methods  #
#==================#

sub new {
    my $class = shift;
    my $self = {};

    bless $self, $class;
    $self->_init (@_);

    return $self;
}


sub setData {
  # check for equal or non-numeric data
  # can receive data either as \@xdata, \@ydata or as @xydata.
  # set refs:  either to $self->{'_xdata'} and $self->{'_ydata'}
  #            or to $self->{'_xydata'}
  # then set $self->{'_flatDataArray'}
  my $self = shift;
  my ($arrayref1, $arrayref2) = @_;
  my ($arrayref, $i);

  if (ref $arrayref2) { # passing data as two data arrays (x0 ...) (y0 ...)

    unless ($#$arrayref1 == $#$arrayref2) { # error checking
      $self->{'_errorMessage'} = "The dataset does not contain an equal number of x and y values. ";
      return 0;
    }

    unless ($#$arrayref1 > 1) { # error checking
      $self->{'_errorMessage'} = "The data set must contain at least three points. ";
      return 0;
    }

    # check whether data are equal and numeric
    for ($i=0; $i<=$#$arrayref1; $i++) {

      unless ($$arrayref1[$i] =~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/) {
	$self->{'_errorMessage'} = "The data element $$arrayref1[$i] is non-numeric. ";
	return 0;
      }
      unless ($$arrayref2[$i] =~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/) {
	$self->{'_errorMessage'} = "The data element $$arrayref2[$i] is non-numeric. ";
	return 0;
      }
    }

    $self->{'_xdata'} = $arrayref1;
    $self->{'_ydata'} = $arrayref2;
    $self->{'_flatDataArray'} = 0; # passed as two data arrays

  } else { # passing data as a single flat data array (x0 y0 ...)

    # check whether array is unbalanced
    if ($#$arrayref1 % 2 == 0) {
      $self->{'_errorMessage'} = "The dataset does not contain an equal number of x and y values.";
      return 0;
    }

    unless ($#$arrayref1 > 4) { # error checking
      $self->{'_errorMessage'} = "The data set must contain at least three points. ";
      return 0;
    }

    # check whether data are numeric
    for ($i=0; $i<=$#$arrayref1; $i++) {
      unless ($$arrayref1[$i] =~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/) {
	$self->{'_errorMessage'} = "The data element $$arrayref1[$i] is non-numeric.";
	return 0;
      }
    }

    $self->{'_xydata'} = $arrayref1;
    $self->{'_flatDataArray'} = 1; # passed as one data array
  }

  $self->{'_dataIsSet'} = 1;
  $self->{'_gotMinMax'} = 0; # recalculate min-max if already calculated
  return 1;
}



sub error {
  # returns the last error message as a string
  my $self = shift;
  return $self->{'_errorMessage'};
}


sub regress {
  my $self = shift;

  unless ($self->{'_dataIsSet'}) {
      $self->{'_errorMessage'} = "No datset has been registered. ";
      return 0;
  }

  my ($sumX, $sumY, $sumXX, $sumYY, $sumXY) = qw (0 0 0 0 0);
  my ($n, $i, $arrayref);

  if ($self->{'_flatDataArray'}) {
    $arrayref = $self->{'_xydata'};
    $n = 1 + $#{ $arrayref };
    for ($i=0; $i<$n; $i+=2) {
      $sumX += $self->{'_xydata'}[$i];
      $sumY += $self->{'_xydata'}[$i+1];
      $sumXX += $self->{'_xydata'}[$i]**2;
      $sumYY += $self->{'_xydata'}[$i+1]**2;
      $sumXY += $self->{'_xydata'}[$i] * $self->{'_xydata'}[$i+1];
    }
    $n /= 2; # number of observations
  } else {
    $arrayref = $self->{'_xdata'};
    $n = $#{ $arrayref };
    $n++; # number of observations
    for ( $i=0; $i<$n; $i++ ) {
      $sumX += $self->{'_xdata'}[$i];
      $sumY += $self->{'_ydata'}[$i];
      $sumXX += $self->{'_xdata'}[$i]**2;
      $sumYY += $self->{'_ydata'}[$i]**2;
      $sumXY += $self->{'_xdata'}[$i] * $self->{'_ydata'}[$i];
    }
  }

  # sum of squared deviations of X and Y
  $self->{'_ssdX'} = $sumXX - $sumX**2/$n;
  $self->{'_ssdY'} = $sumYY - $sumY**2/$n;
  $self->{'_ssdXY'} = $sumXY - $sumX*$sumY/$n;

  # num observations and sample averages
  $self->{'_n'} = $n;
  ($self->{'_avX'}, $self->{'_avY'}) = ($sumX/$n, $sumY/$n);

  # sample var's and cov's (using n-1)
  $self->{'_varX'} = $self->{'_ssdX'} / ($n-1);
  $self->{'_varY'} = $self->{'_ssdY'} / ($n-1);
  $self->{'_covXY'} = $self->{'_ssdXY'} / ($n-1);

  # coefficient estimates
  $self->{'_b2'} = $self->{'_ssdX'} == 0
    ? undef
      : $self->{'_ssdXY'} / $self->{'_ssdX'}; # slope
  $self->{'_b1'} = ($sumY - $self->{'_b2'} * $sumX) / $n; # intercept

  # R-squared
  $self->{'_rsq'} = ($self->{'_ssdX'} == 0 or $self->{'_ssdY'} == 0)
    ? 1.0
      : ($self->{'_ssdXY'} / $self->{'_ssdX'})
	* ($self->{'_ssdXY'} / $self->{'_ssdY'}) ;
#  $self->{'_rsq'} = $self->{'_b2'}**2 * $self->{'_ssdX'} /  $self->{'_ssdY'};

  # error (residual) sum of squares
  $self->{'_sse'} = $self->{'_ssdY'} - $self->{'_ssdX'} * $self->{'_b2'}**2;
  $self->{'_sse'} = 0 if $self->{'_sse'} < 0; # potential precision problems

  # homoscedastic standard deviation of error term
  $self->{'_sigma'} = sqrt ($self->{'_sse'}/($n-2));

  # standard error of coefficients and t-stats
  $self->{'_seB1'} = $self->{'_seB2'} = undef;
  $self->{'_t1'} = $self->{'_t2'} = undef;

  unless ($self->{'_ssdX'} == 0) {
    $self->{'_seB1'} = $self->{'_sigma'} * sqrt ($sumXX / ($n*$self->{'_ssdX'}));
    $self->{'_seB2'} = $self->{'_sigma'} / sqrt $self->{'_ssdX'};
    $self->{'_t2'} = $self->{'_b2'} / $self->{'_seB2'} unless $self->{'_seB2'} == 0;
    $self->{'_t1'} = $self->{'_b1'} / $self->{'_seB1'} unless $self->{'_seB1'} == 0;
  }

  # durbin-watson
  my $sum = 0;
  my ($prevErr, $currentErr);

  if ($self->{'_sse'} == 0) {
    $self->{'_dw'} = undef;
  } else {
    if ($self->{'_flatDataArray'}) {
      $arrayref = $self->{'_xydata'};
      $n = 1+$#{ $arrayref };
      $prevErr = $self->{'_xydata'}[1]
	- $self->{'_b1'} - $self->{'_b2'} * $self->{'_xydata'}[0];
      for ($i=2; $i<$n; $i+=2) {
	$currentErr = $self->{'_xydata'}[$i+1]
	  - $self->{'_b1'} - $self->{'_b2'} * $self->{'_xydata'}[$i];
	$sum += ($currentErr - $prevErr)**2;
	$prevErr = $currentErr;
      }
    } else {
      $arrayref = $self->{'_xdata'};
      $n = 1+$#{ $arrayref };
      $prevErr = $self->{'_ydata'}[0]
	- $self->{'_b1'} - $self->{'_b2'} * $self->{'_xdata'}[0];
      for ( $i=1; $i<$n; $i++ ) {
	$currentErr = $self->{'_ydata'}[$i]
	  - $self->{'_b1'} - $self->{'_b2'} * $self->{'_xdata'}[$i];
	$sum += ($currentErr - $prevErr)**2;
	$prevErr = $currentErr;
      }
    }
    $self->{'_dw'} = $sum / $self->{'_sse'};
  }

  $self->{'_gotMinMax'} = 0; # should recalculate min-max's if already calculated
  return 1;
}

sub minMax {
  my $self = shift;
  $self->_getMinMax() unless $self->{'_gotMinMax'};
  return ($self->{'_xmin'}, $self->{'_xmax'},
	  $self->{'_ymin'}, $self->{'_ymax'});
}

sub coefficients {  my $self = shift;  return ($self->{'_b1'}, $self->{'_b2'}); }

sub rsq {  my $self = shift;  return $self->{'_rsq'}; }

sub tstats {  my $self = shift;  return ($self->{'_t1'}, $self->{'_t2'}); }

sub av {  my $self = shift;  return ($self->{'_avX'}, $self->{'_avY'}); }

sub var {  my $self = shift;  return ($self->{'_varX'}, $self->{'_varY'},
				      $self->{'_covXY'}); }

sub sigma {  my $self = shift;  return $self->{'_sigma'}; }

sub size {   my $self = shift;  return $self->{'_n'}; }

sub dw {   my $self = shift;  return $self->{'_dw'}; }

sub residuals {
  my $self = shift;
  my ($n, $i, $arrayref);
  my @result = ();

  if ($self->{'_flatDataArray'}) { # construct xy data array
    $arrayref = $self->{'_xydata'};
    $n = 1+$#{ $arrayref };
    for ($i=0; $i<$n; $i+=2) {
      $result[$i] = $self->{'_xydata'}[$i];
      $result[$i+1] = $self->{'_xydata'}[$i+1]
                      - $self->{'_b1'} - $self->{'_b2'} * $self->{'_xydata'}[$i];
    }
  } else { # construct y data array
    $arrayref = $self->{'_xdata'};
    $n = 1+$#{ $arrayref };
    for ( $i=0; $i<$n; $i++ ) {
      $result[$i] = $self->{'_ydata'}[$i]
                    - $self->{'_b1'} - $self->{'_b2'} * $self->{'_xdata'}[$i];
    }
  }
  return @result;
}

sub predicted {
  my $self = shift;
  my ($n, $i, $arrayref);
  my @result = ();

  if ($self->{'_flatDataArray'}) {
    $arrayref = $self->{'_xydata'};
    $n = 1+$#{ $arrayref };
    for ($i=0; $i<$n; $i+=2) {
      $result[$i] = $self->{'_xydata'}[$i];
      $result[$i+1] = $self->{'_b1'} + $self->{'_b2'} * $self->{'_xydata'}[$i];
    }
  } else {
    $arrayref = $self->{'_xdata'};
    $n = 1+$#{ $arrayref };
    for ( $i=0; $i<$n; $i++ ) {
      $result[$i] = $self->{'_b1'} + $self->{'_b2'} * $self->{'_xdata'}[$i];
    }
  }
  return @result;
}


#===================#
#  private methods  #
#===================#

# initialization;
# this contains a record of all private data
# this is the place to start if you want to read the code.
sub _init {
  my $self = shift;

# $self->{'_flatDataArray'} = ''; # data passed as one flat or two data arrays?
  $self->{'_dataIsSet'} = 0; # return error if asking to regress
  $self->{'_errorMessage'} = '';

  # will hold references to caller's data array(s)
# $self->{'_xydata'} = $self->{'_xdata'} = $self->{'_ydata'} = '';

# $self->{'_ssdX'} = $self->{'_ssdY'} = $self->{'_ssdXY'} = '';
  $self->{'_n'} = 0; # num observations
# $self->{'_avX'} = $self->{'_avY'} = '';
# $self->{'_varX'} =  $self->{'_vary'} =  $self->{'_covXY'} = '';

  $self->{'_gotMinMax'} = 0; # do not calculate again
# $self->{'_xmin'} = $self->{'_xmax'} = 0;
# $self->{'_ymin'} = $self->{'_ymax'} = 0;

# $self->{'_b1'} = $self->{'_b2'} = '';
# $self->{'_rsq'} = $self->{'_sse'} = $self->{'_sigma'} = '';

# $self->{'_seB1'} = $self->{'_seB2'} = undef;
# $self->{'_t1'} = $self->{'_t2'} = undef;
# $self->{'_dw'} = undef;
}


# sets min and max values of all data (_xmin, _ymin, _xmax, _ymax);
# also sets _xslope, _yslope, _ax and _ay used in _data2pxl;
# usage: $self->_getMinMax
sub _getMinMax {
  my $self = shift;
  my ($i, $n, $arrayref);

  if ($self->{'_flatDataArray'}) {
    $self->{'_xmin'} = $self->{'_xmax'} = $self->{'_xydata'}[0];
    $self->{'_ymin'} = $self->{'_ymax'} = $self->{'_xydata'}[1];
    $arrayref = $self->{'_xydata'};
    $n = 1+$#{ $arrayref };
    for ($i=2; $i<$n; $i+=2) {
      $self->{'_xmin'} = $self->{'_xydata'}[$i]
	if $self->{'_xydata'}[$i] < $self->{'_xmin'};
      $self->{'_xmax'} = $self->{'_xydata'}[$i]
	if $self->{'_xydata'}[$i] > $self->{'_xmax'};
      $self->{'_ymin'} = $self->{'_xydata'}[$i+1]
	if $self->{'_xydata'}[$i+1] < $self->{'_ymin'};
      $self->{'_ymax'} = $self->{'_xydata'}[$i+1]
	if $self->{'_xydata'}[$i+1] > $self->{'_ymax'};
    }
    $n /= 2; # number of observations
  } else {
    $self->{'_xmin'} = $self->{'_xmax'} = $self->{'_xdata'}[0];
    $self->{'_ymin'} = $self->{'_ymax'} = $self->{'_ydata'}[0];
    $arrayref = $self->{'_xdata'};
    $n = 1+$#{ $arrayref };
    for ( $i=1; $i<$n; $i++ ) {
      $self->{'_xmin'} = $self->{'_xdata'}[$i]
	if $self->{'_xdata'}[$i] < $self->{'_xmin'};
      $self->{'_xmax'} = $self->{'_xdata'}[$i]
	if $self->{'_xdata'}[$i] > $self->{'_xmax'};
      $self->{'_ymin'} = $self->{'_ydata'}[$i]
	if $self->{'_ydata'}[$i] < $self->{'_ymin'};
      $self->{'_ymax'} = $self->{'_ydata'}[$i]
	if $self->{'_ydata'}[$i] > $self->{'_ymax'};
    }
  }
  $self->{'_gotMinMax'} = 1;
}

1;

__END__