File Coverage

blib/lib/Statistics/LineFit.pm

Criterion	Covered	Total	%
statement	211	215	98.1
branch	119	138	86.2
condition	7	12	58.3
subroutine	22	22	100.0
pod	12	17	70.5
total	371	404	91.8

package Statistics::LineFit;
use strict;
use Carp qw(carp);
BEGIN {
        use Exporter ();
        use vars qw ($AUTHOR $VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
        $AUTHOR      = 'Richard Anderson ';
        @EXPORT      = @EXPORT_OK = qw();
        %EXPORT_TAGS = ();
        @ISA         = qw(Exporter);
        $VERSION     = 0.06;
}

sub new {
#
# Purpose: Create a new Statistics::LineFit object
#
    my ($caller, $validate, $hush) = @_;
    my $self = { doneRegress  => 0,
                 gotData      => 0,
                 hush         => defined $hush ? $hush : 0,
                 validate     => defined $validate ? $validate : 0,
               };
    bless $self, ref($caller) || $caller;
    return $self;
}

sub coefficients {
#
# Purpose: Return the slope and intercept from least squares line fit
#
    my $self = shift;
    unless (defined $self->{intercept} and defined $self->{slope}) {
        $self->regress() or return;
    }
    return ($self->{intercept}, $self->{slope});
}

sub computeSums {
#
# Purpose: Compute sum of x, y, x**2, y**2 and x*y (private method)
#
    my $self = shift;
    my ($sumX, $sumY, $sumXX, $sumYY, $sumXY) = (0, 0, 0, 0, 0);
    if (defined $self->{weight}) {
        for (my $i = 0; $i < $self->{numXY}; ++$i) {
            $sumX += $self->{weight}[$i] * $self->{x}[$i];
            $sumY += $self->{weight}[$i] * $self->{y}[$i];
            $sumXX += $self->{weight}[$i] * $self->{x}[$i] ** 2;
            $sumYY += $self->{weight}[$i] * $self->{y}[$i] ** 2;
            $sumXY += $self->{weight}[$i] * $self->{x}[$i]
                * $self->{y}[$i];
        }
    } else {
        for (my $i = 0; $i < $self->{numXY}; ++$i) {
            $sumX += $self->{x}[$i];
            $sumY += $self->{y}[$i];
            $sumXX += $self->{x}[$i] ** 2;
            $sumYY += $self->{y}[$i] ** 2;
            $sumXY += $self->{x}[$i] * $self->{y}[$i];
        }
    }
    return ($sumX, $sumY, $sumXX, $sumYY, $sumXY);
}

sub durbinWatson {
#
# Purpose: Return the Durbin-Watson statistic
#
    my $self = shift;
    unless (defined $self->{durbinWatson}) {
        $self->regress() or return;
        my $sumErrDiff = 0;
        my $errorTMinus1 = $self->{y}[0] - ($self->{intercept} + $self->{slope}
            * $self->{x}[0]);
        for (my $i = 1; $i < $self->{numXY}; ++$i) {
            my $error = $self->{y}[$i] - ($self->{intercept} + $self->{slope}
                * $self->{x}[$i]);
            $sumErrDiff += ($error - $errorTMinus1) ** 2;
            $errorTMinus1 = $error;
        }
        $self->{durbinWatson} = $self->sumSqErrors() > 0 ?
            $sumErrDiff / $self->sumSqErrors() : 0;
    }
    return $self->{durbinWatson};
}

sub meanSqError {
#
# Purpose: Return the mean squared error
#
    my $self = shift;
    unless (defined $self->{meanSqError}) {
        $self->regress() or return;
        $self->{meanSqError} = $self->sumSqErrors() / $self->{numXY};
    }
    return $self->{meanSqError};
}

sub predictedYs {
#
# Purpose: Return the predicted y values
#
    my $self = shift;
    unless (defined $self->{predictedYs}) {
        $self->regress() or return;
        $self->{predictedYs} = [];
        for (my $i = 0; $i < $self->{numXY}; ++$i) {
            $self->{predictedYs}[$i] = $self->{intercept}
                + $self->{slope} * $self->{x}[$i];
        }
    }
    return @{$self->{predictedYs}};
}

sub regress {
#
# Purpose: Do weighted or unweighted least squares 2-D line fit (if needed)
#
# Description:
# The equations below apply to both the weighted and unweighted fit: the
# weights are normalized in setWeights(), so the sum of the weights is
# equal to numXY.
#
    my $self = shift;
    return $self->{regressOK} if $self->{doneRegress};
    unless ($self->{gotData}) {
        carp "No valid data input - can't do regression" unless $self->{hush};
        return 0;
    }
    my ($sumX, $sumY, $sumYY, $sumXY);
    ($sumX, $sumY, $self->{sumXX}, $sumYY, $sumXY) = $self->computeSums();
    $self->{sumSqDevX} = $self->{sumXX} - $sumX ** 2 / $self->{numXY};
    if ($self->{sumSqDevX} != 0) {
        $self->{sumSqDevY} = $sumYY - $sumY ** 2 / $self->{numXY};
        $self->{sumSqDevXY} = $sumXY - $sumX * $sumY / $self->{numXY};
        $self->{slope} = $self->{sumSqDevXY} / $self->{sumSqDevX};
        $self->{intercept} = ($sumY - $self->{slope} * $sumX) / $self->{numXY};
        $self->{regressOK} = 1;
    } else {
        carp "Can't fit line when x values are all equal" unless $self->{hush};
        $self->{sumXX} = $self->{sumSqDevX} = undef;
        $self->{regressOK} = 0;
    }
    $self->{doneRegress} = 1;
    return $self->{regressOK};
}

sub residuals {
#
# Purpose: Return the predicted Y values minus the observed Y values
#
    my $self = shift;
    unless (defined $self->{residuals}) {
        $self->regress() or return;
        $self->{residuals} = [];
        for (my $i = 0; $i < $self->{numXY}; ++$i) {
            $self->{residuals}[$i] = $self->{y}[$i] - ($self->{intercept}
                + $self->{slope} * $self->{x}[$i]);
        }
    }
    return @{$self->{residuals}};
}

sub rSquared {
#
# Purpose: Return the correlation coefficient
#
    my $self = shift;
    unless (defined $self->{rSquared}) {
        $self->regress() or return;
        my $denom = $self->{sumSqDevX} * $self->{sumSqDevY};
        $self->{rSquared} = $denom != 0 ? $self->{sumSqDevXY} ** 2 / $denom : 1;
    }
    return $self->{rSquared};
}

sub setData {
#
# Purpose: Initialize (x,y) values and optional weights
#
    my ($self, $x, $y, $weights) = @_;
    $self->{doneRegress} = 0;
    $self->{x} = $self->{y} = $self->{numXY} = $self->{weight}
        = $self->{intercept} = $self->{slope} = $self->{rSquared}
        = $self->{sigma} = $self->{durbinWatson} = $self->{meanSqError}
        = $self->{sumSqErrors} = $self->{tStatInt} = $self->{tStatSlope}
        = $self->{predictedYs} = $self->{residuals} = $self->{sumXX}
        = $self->{sumSqDevX} = $self->{sumSqDevY} = $self->{sumSqDevXY}
        = undef;
    if (@$x < 2) {
        carp "Must input more than one data point!" unless $self->{hush};
        return 0;
    }
    $self->{numXY} = @$x;
    if (ref $x->[0]) {
        $self->setWeights($y) or return 0;
        $self->{x} = [ ];
        $self->{y} = [ ];
        foreach my $xy (@$x) {
            push @{$self->{x}}, $xy->[0];
            push @{$self->{y}}, $xy->[1];
        }
    } else {
        if (@$x != @$y) {
            carp "Length of x and y arrays must be equal!" unless $self->{hush};
            return 0;
        }
        $self->setWeights($weights) or return 0;
        $self->{x} = [ @$x ];
        $self->{y} = [ @$y ];
    }
    if ($self->{validate}) {
        unless ($self->validData()) {
            $self->{x} = $self->{y} = $self->{weights} = $self->{numXY} = undef;
            return 0;
        }
    }
    $self->{gotData} = 1;
    return 1;
}

sub setWeights {
#
# Purpose: Normalize and initialize line fit weighting factors (private method)
#
    my ($self, $weights) = @_;
    return 1 unless defined $weights;
    if (@$weights != $self->{numXY}) {
        carp "Length of weight array must equal length of data array!"
            unless $self->{hush};
        return 0;
    }
    if ($self->{validate}) { $self->validWeights($weights) or return 0 }
    my $sumW = my $numNonZero = 0;
    foreach my $weight (@$weights) {
        if ($weight < 0) {
            carp "Weights must be non-negative numbers!" unless $self->{hush};
            return 0;
        }
        $sumW += $weight;
        if ($weight != 0) { ++$numNonZero }
    }
    if ($numNonZero < 2) {
        carp "At least two weights must be nonzero!" unless $self->{hush};
        return 0;
    }
    my $factor = @$weights / $sumW;
    foreach my $weight (@$weights) { $weight *= $factor }
    $self->{weight} = [ @$weights ];
    return 1;
}

sub sigma {
#
# Purpose: Return the estimated homoscedastic standard deviation of the
#          error term
#
    my $self = shift;
    unless (defined $self->{sigma}) {
        $self->regress() or return;
        $self->{sigma} = $self->{numXY} > 2 ?
            sqrt($self->sumSqErrors() / ($self->{numXY} - 2)) : 0;
    }
    return $self->{sigma};
}

sub sumSqErrors {
#
# Purpose: Return the sum of the squared errors (private method)
#
    my $self = shift;
    unless (defined $self->{sumSqErrors}) {
        $self->regress() or return;
        $self->{sumSqErrors} = $self->{sumSqDevY} - $self->{sumSqDevX}
            * $self->{slope} ** 2;
        if ($self->{sumSqErrors} < 0) { $self->{sumSqErrors} = 0 }
    }
    return $self->{sumSqErrors};
}

sub tStatistics {
#
# Purpose: Return the T statistics
#
    my $self = shift;
    unless (defined $self->{tStatInt} and defined $self->{tStatSlope}) {
        $self->regress() or return;
        my $biasEstimateInt = $self->sigma() * sqrt($self->{sumXX}
            / ($self->{sumSqDevX} * $self->{numXY}));
        $self->{tStatInt} = $biasEstimateInt != 0 ?
            $self->{intercept} / $biasEstimateInt : 0;
        my $biasEstimateSlope = $self->sigma() / sqrt($self->{sumSqDevX});
        $self->{tStatSlope} = $biasEstimateSlope != 0 ?
            $self->{slope} / $biasEstimateSlope : 0;
    }
    return ($self->{tStatInt}, $self->{tStatSlope});
}

sub validData {
#
# Purpose: Verify that the input x-y data are numeric (private method)
#
    my $self = shift;
    for (my $i = 0; $i < $self->{numXY}; ++$i) {
        if (not defined $self->{x}[$i]) {
            carp "Input x[$i] is not defined" unless $self->{hush};
            return 0;
        }
        if ($self->{x}[$i] !~
            /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/)
        {
            carp "Input x[$i] is not a number: $self->{x}[$i]"
                unless $self->{hush};
            return 0;
        }
        if (not defined $self->{y}[$i]) {
            carp "Input y[$i] is not defined" unless $self->{hush};
            return 0;
        }
        if ($self->{y}[$i] !~
            /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/)
        {
            carp "Input y[$i] is not a number: $self->{y}[$i]"
                unless $self->{hush};
            return 0;
        }
    }
    return 1;
}

sub validWeights {
#
# Purpose: Verify that the input weights are numeric (private method)
#
    my ($self, $weights) = @_;
    for (my $i = 0; $i < @$weights; ++$i) {
        if (not defined $weights->[$i]) {
            carp "Input weights[$i] is not defined" unless $self->{hush};
            return 0;
        }
        if ($weights->[$i]
            !~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/)
        {
            carp "Input weights[$i] is not a number: $weights->[$i]"
                unless $self->{hush};
            return 0;
        }
    }
    return 1;
}

sub varianceOfEstimates {
#
# Purpose: Return the variances in the estimates of the intercept and slope
#
    my $self = shift;
    unless (defined $self->{intercept} and defined $self->{slope}) {
        $self->regress() or return;
    }
    my @predictedYs = $self->predictedYs();
    my ($s, $sx, $sxx) = (0, 0, 0);
    if (defined $self->{weight}) {
        for (my $i = 0; $i < $self->{numXY}; ++$i) {
            my $variance = ($predictedYs[$i] - $self->{y}[$i]) ** 2;
            next if 0 == $variance;
            $s += 1.0 / $variance;
	    $sx += $self->{weight}[$i] * $self->{x}[$i] / $variance;
	    $sxx += $self->{weight}[$i] * $self->{x}[$i] ** 2 / $variance;
        }
    } else {
        for (my $i = 0; $i < $self->{numXY}; ++$i) {
            my $variance = ($predictedYs[$i] - $self->{y}[$i]) ** 2;
            next if 0 == $variance;
            $s += 1.0 / $variance;
	    $sx += $self->{x}[$i] / $variance;
	    $sxx += $self->{x}[$i] ** 2 / $variance;
        }
    }
    my $denominator = ($s * $sxx - $sx ** 2);
    if (0 == $denominator) {
        return;
    } else {
        return ($sxx / $denominator, $s / $denominator);
    }
}

1;

__END__