Statistics::Descriptive - Module of basic descriptive statistical functions.

  use Statistics::Descriptive;
  $stat = Statistics::Descriptive::Full->new();
  $stat->add_data(1,2,3,4); $mean = $stat->mean();
  $var  = $stat->variance();
  $tm   = $stat->trimmed_mean(.25);
  $Statistics::Descriptive::Tolerance = 1e-10;

package Statistics::Descriptive;

use strict;
use warnings;

##This module draws heavily from perltoot v0.4 from Tom Christiansen.

require 5.00404;  ##Yes, this is underhanded, but makes support for me easier
		  ##Not only that, but it's the latest "safe" version of
		  ##Perl5.  01-03 weren't bug free.
use vars (qw($VERSION $Tolerance));

$VERSION = '3.0201';

$Tolerance = 0.0;

package Statistics::Descriptive::Sparse;

use vars qw($VERSION);

$VERSION = '3.0201';

use vars qw(%fields);
use Carp;

sub _make_accessors
{
    my ($pkg, $methods) = @_;

    no strict 'refs';
    foreach my $method (@$methods)
    {
        *{$pkg."::".$method} =
            do {
                my $m = $method;
                sub {
                    my $self = shift;

                    if (@_)
                    {
                        $self->{$m} = shift;
                    }
                    return $self->{$m};
                };
            };
    }

    return;
}

sub _make_private_accessors
{
    my ($pkg, $methods) = @_;

    no strict 'refs';
    foreach my $method (@$methods)
    {
        *{$pkg."::_".$method} =
            do {
                my $m = $method;
                sub {
                    my $self = shift;

                    if (@_)
                    {
                        $self->{$m} = shift;
                    }
                    return $self->{$m};
                };
            };
    }

    return;
}

##Define the fields to be used as methods
%fields = (
  count			=> 0,
  mean			=> 0,
  sum			=> 0,
  sumsq			=> 0,
  min			=> undef,
  max			=> undef,
  mindex		=> undef,
  maxdex		=> undef,
  sample_range		=> undef,
  variance => undef,
  );

__PACKAGE__->_make_accessors( [ grep { $_ ne "variance" } keys(%fields) ] );
__PACKAGE__->_make_accessors( ["_permitted"] );
__PACKAGE__->_make_private_accessors(["variance"]);

sub new {
  my $proto = shift;
  my $class = ref($proto) || $proto;
  my $self = {
    %fields,
  };
  bless ($self, $class);
  $self->_permitted(\%fields);
  return $self;
}

sub _is_permitted
{
    my $self = shift;
    my $key = shift;

    return exists($self->_permitted()->{$key});
}

sub add_data {
  my $self = shift;  ##Myself
  my $oldmean;
  my ($min,$mindex,$max,$maxdex,$sum,$sumsq,$count);
  my $aref;

  if (ref $_[0] eq 'ARRAY') {
    $aref = $_[0];
  }
  else {
    $aref = \@_;
  }

  ##If we were given no data, we do nothing.
  return 1 if (!@{ $aref });

  ##Take care of appending to an existing data set
  
  if (!defined($min = $self->min()))
  {
      $min = $aref->[$mindex = 0];
  }
  else
  {
      $mindex = $self->mindex();
  }

  if (!defined($max = $self->max()))
  {
      $max = $aref->[$maxdex = 0];
  }
  else
  {
      $maxdex = $self->maxdex();
  }

  $sum = $self->sum();
  $sumsq = $self->sumsq();
  $count = $self->count();

  ##Calculate new mean, sumsq, min and max;
  foreach ( @{ $aref } ) {
    $sum += $_;
    $sumsq += $_**2;
    $count++;
    if ($_ >= $max) {
      $max = $_;
      $maxdex = $count-1;
    }
    if ($_ <= $min) {
      $min = $_;
      $mindex = $count-1;
    }
  }

  $self->min($min);
  $self->mindex($mindex);
  $self->max($max);
  $self->maxdex($maxdex);
  $self->sample_range($max - $min);
  $self->sum($sum);
  $self->sumsq($sumsq);
  $self->mean($sum / $count);
  $self->count($count);
  ##indicator the value is not cached.  Variance isn't commonly enough
  ##used to recompute every single data add.
  $self->_variance(undef());
  return 1;
}

sub standard_deviation {
  my $self = shift;  ##Myself
  return undef if (!$self->count());
  return sqrt($self->variance());
}

##Return variance; if needed, compute and cache it.
sub variance {
  my $self = shift;  ##Myself
  my $div = @_ ? 0 : 1;
  my $count = $self->count();
  if ($count < 1 + $div) {
      return 0;
  }

  if (!defined($self->_variance())) {
    my $variance = ($self->sumsq()- $count * $self->mean()**2);

    # Sometimes due to rounding errors we get a number below 0.
    # This makes sure this is handled as gracefully as possible.
    #
    # See:
    #
    # https://rt.cpan.org/Public/Bug/Display.html?id=46026
    if ($variance < 0)
    {
        $variance = 0;
    }

    $variance /= $count - $div;

    $self->_variance($variance);
  }
  return $self->_variance();
}

##Clear a stat.  More efficient than destroying an object and calling
##new.
sub clear {
  my $self = shift;  ##Myself
  my $key;

  return if (!$self->count());
  while (my($field, $value) = each %fields) {
    $self->{$field} = $value;
  }
}

1;

package Statistics::Descriptive::Full;

use vars qw($VERSION);

$VERSION = '3.0201';

use Carp;

use POSIX ();

use vars qw(@ISA $a $b %fields);

@ISA = qw(Statistics::Descriptive::Sparse);

##Create a list of fields not to remove when data is updated
%fields = (
  _permitted => undef,  ##Place holder for the inherited key hash
  data       => undef,  ##Our data
  presorted  => undef,  ##Flag to indicate the data is already sorted
  _reserved  => undef,  ##Place holder for this lookup hash
);

__PACKAGE__->_make_private_accessors(
    [qw(data frequency geometric_mean harmonic_mean 
        least_squares_fit median mode
        skewness kurtosis
       )
    ]
);
__PACKAGE__->_make_accessors([qw(presorted _reserved _trimmed_mean_cache)]);

sub _clear_fields
{
    my $self = shift;

    # Empty array ref for holding data later!
    $self->_data([]);
    $self->_reserved(\%fields);
    $self->presorted(0);
    $self->_trimmed_mean_cache(+{});

    return;
}

##Have to override the base method to add the data to the object
##The proxy method from above is still valid
sub new {
  my $proto = shift;
  my $class = ref($proto) || $proto;
  # Create my self re SUPER
  my $self = $class->SUPER::new();  
  bless ($self, $class);  #Re-anneal the object
  $self->_clear_fields();
  return $self;
}

sub _is_reserved
{
    my $self = shift;
    my $field = shift;

    return exists($self->_reserved->{$field});
}

sub _delete_all_cached_keys
{
    my $self = shift;

    KEYS_LOOP:
    foreach my $key (keys %{ $self }) { # Check each key in the object
        # If it's a reserved key for this class, keep it
        if ($self->_is_reserved($key) || $self->_is_permitted($key))
        {
            next KEYS_LOOP;
        }
        delete $self->{$key};          # Delete the out of date cached key
    }
    return;
}

##Clear a stat.  More efficient than destroying an object and calling
##new.
sub clear {
    my $self = shift;  ##Myself
    my $key;

    if (!$self->count())
    {
        return;
    }

    $self->_delete_all_cached_keys();
    $self->SUPER::clear();
    $self->_clear_fields();
}

sub add_data {
  my $self = shift;
  my $aref;

  if (ref $_[0] eq 'ARRAY') {
    $aref = $_[0];
  }
  else {
    $aref = \@_;
  }
  $self->SUPER::add_data($aref);  ##Perform base statistics on the data
  push @{ $self->_data() }, @{ $aref };
  ##Clear the presorted flag
  $self->presorted(0);

  $self->_delete_all_cached_keys();

  return 1;
}

sub get_data {
  my $self = shift;
  return @{ $self->_data() };
}

sub sort_data {
  my $self = shift;

  if (! $self->presorted())
  {
      ##Sort the data in descending order
      $self->_data([ sort {$a <=> $b} @{$self->_data()} ]);
      $self->presorted(1);
      ##Fix the maxima and minima indices
      $self->mindex(0);
      $self->maxdex($#{$self->_data()});
  }

  return 1;
}

sub percentile {
  my $self = shift;
  my $percentile = shift || 0;
  ##Since we're returning a single value there's no real need
  ##to cache this.

  ##If the requested percentile is less than the "percentile bin
  ##size" then return undef.  Check description of RFC 2330 in the
  ##POD below.
  my $count = $self->count();
  return undef if $percentile < 100 / $count;

  $self->sort_data();
  my $num = $count*$percentile/100;
  my $index = &POSIX::ceil($num) - 1;
  my $val = $self->_data->[$index];
  return wantarray
    ? ($val, $index)
    : $val
    ;
}

sub _calc_new_median
{
    my $self = shift;
    my $count = $self->count();

    ##Even or odd
    if ($count % 2)
    {   
        return $self->_data->[($count-1)/2];
    }
    else
    {
        return
        (
            ($self->_data->[($count)/2] + $self->_data->[($count-2)/2] ) / 2
        );
    }
}

sub median {
    my $self = shift;

    ##Cached?
    if (! defined($self->_median()))
    {
        $self->sort_data();
        $self->_median($self->_calc_new_median());
    }
    return $self->_median();
}

sub quantile {
    my ( $self, $QuantileNumber ) = @_;

    unless ( defined $QuantileNumber and $QuantileNumber =~ m/^0|1|2|3|4$/ ) {
       carp("Bad quartile type, must be 0, 1, 2, 3 or 4\n");
       return;
    }
    
    $self->sort_data();

    return $self->_data->[0] if ( $QuantileNumber == 0 );

    my $count = $self->count();

    return $self->_data->[ $count - 1 ] if ( $QuantileNumber == 4 );

    my $K_quantile = ( ( $QuantileNumber / 4 ) * ( $count - 1 ) + 1 );
    my $F_quantile = $K_quantile - POSIX::floor($K_quantile);
    $K_quantile = POSIX::floor($K_quantile);

    # interpolation
    my $aK_quantile     = $self->_data->[ $K_quantile - 1 ];
    return $aK_quantile if ( $F_quantile == 0 );
    my $aKPlus_quantile = $self->_data->[$K_quantile];
    
    # Calcul quantile
    my $quantile = $aK_quantile
      + ( $F_quantile * ( $aKPlus_quantile - $aK_quantile ) );

    return $quantile;
}

sub _real_calc_trimmed_mean
{
    my $self = shift;
    my $lower = shift;
    my $upper = shift;

    my $lower_trim = int ($self->count()*$lower); 
    my $upper_trim = int ($self->count()*$upper); 
    my ($val,$oldmean) = (0,0);
    my ($tm_count,$tm_mean,$index) = (0,0,$lower_trim);

    $self->sort_data();
    while ($index <= $self->count() - $upper_trim -1)
    {
        $val = $self->_data()->[$index];
        $oldmean = $tm_mean;
        $index++;
        $tm_count++;
        $tm_mean += ($val - $oldmean) / $tm_count;
    }

    return $tm_mean;
}

sub trimmed_mean
{
    my $self = shift;
    my ($lower,$upper);
    #upper bound is in arg list or is same as lower
    if (@_ == 1)
    {
        ($lower,$upper) = ($_[0],$_[0]);
    }
    else
    {
        ($lower,$upper) = ($_[0],$_[1]);
    }

    ##Cache
    my $thistm = join ':',$lower,$upper;
    my $cache = $self->_trimmed_mean_cache();
    if (!exists($cache->{$thistm}))
    {
        $cache->{$thistm} = $self->_real_calc_trimmed_mean($lower, $upper);
    }

    return $cache->{$thistm};
}

sub _test_for_too_small_val
{
    my $self = shift;
    my $val = shift;

    return (abs($val) <= $Statistics::Descriptive::Tolerance);
}

sub _calc_harmonic_mean
{
    my $self = shift;

    my $hs = 0;

    foreach my $item ( @{$self->_data()} )
    {
        ##Guarantee that there are no divide by zeros
        if ($self->_test_for_too_small_val($item))
        {
            return;
        }

        $hs += 1/$item;
    }

    if ($self->_test_for_too_small_val($hs))
    {
        return;
    }

    return $self->count()/$hs;
}

sub harmonic_mean
{
    my $self = shift;

    if (!defined($self->_harmonic_mean()))
    {
        $self->_harmonic_mean(scalar($self->_calc_harmonic_mean()));
    }

    return $self->_harmonic_mean();
}

sub mode
{
    my $self = shift;

    if (!defined ($self->_mode()))
    {
        my $mode = 0;
        my $occurances= 0;
        my $flag = 1;

        my %count;

        foreach my $item (@{ $self->_data() })
        {
            $count{$item}++;
            $flag = 0 if ($count{$item} > 1);
        }

        #Distribution is flat - no mode exists
        if ($flag)
        {
            return undef;
        }

        foreach my $val (keys %count)
        {
            if ($count{$val} > $occurances)
            {
                $occurances = $count{$val};
                $mode = $val;
            }
        }

        $self->_mode($mode);
    }

    return $self->_mode();
}

sub geometric_mean {
    my $self = shift;

    if (!defined($self->_geometric_mean()))
    {
        my $gm = 1;
        my $exponent = 1/$self->count();

        for my $val (@{ $self->_data() })
        {
            if ($val < 0)
            {
                return undef;
            }
            $gm *= $val**$exponent;
        }

        $self->_geometric_mean($gm);
    }

    return $self->_geometric_mean();
}

sub skewness {
    my $self = shift;

    if (!defined($self->_skewness()))
    {
        my $n    = $self->count();
        my $sd   = $self->standard_deviation();
        
        my $skew;
        
        #  skip if insufficient records
        if ( $sd && $n > 2) {
            
            my $mean = $self->mean();
            
            my $sum_pow3;
            
            foreach my $rec ( $self->get_data ) {
                my $value  = (($rec - $mean) / $sd);
                $sum_pow3 +=  $value ** 3;
            }
            
            my $correction = $n / ( ($n-1) * ($n-2) );
            
            $skew = $correction * $sum_pow3;
        }

        $self->_skewness($skew);
    }

    return $self->_skewness();
}

sub kurtosis {
    my $self = shift;

    if (!defined($self->_kurtosis()))
    {
        my $kurt;
        
        my $n  = $self->count();
        my $sd   = $self->standard_deviation();
        
        if ( $sd && $n > 3) {

            my $mean = $self->mean();
            
            my $sum_pow4;
            foreach my $rec ( $self->get_data ) {
                $sum_pow4 += ( ($rec - $mean ) / $sd ) ** 4;
            }
            
            my $correction1 = ( $n * ($n+1) ) / ( ($n-1) * ($n-2) * ($n-3) );
            my $correction2 = ( 3  * ($n-1) ** 2) / ( ($n-2) * ($n-3) );
            
            $kurt = ( $correction1 * $sum_pow4 ) - $correction2;
        }
        
        $self->_kurtosis($kurt);
    }

    return $self->_kurtosis();
}


sub frequency_distribution_ref
{
    my $self = shift;
    my @k = ();
    # Must have at least two elements
    if ($self->count() < 2)
    {
        return undef;
    }

    if ((!@_) && (defined $self->_frequency()))
    {
        return $self->_frequency()
    }

    my %bins;
    my $partitions = shift;

    if (ref($partitions) eq 'ARRAY')
    {
        @k = @{ $partitions };
        return undef unless @k;  ##Empty array
        if (@k > 1) {
            ##Check for monotonicity
            my $element = $k[0];
            for my $next_elem (@k[1..$#k]) {
                if ($element > $next_elem) {
                    carp "Non monotonic array cannot be used as frequency bins!\n";
                    return undef;
                }
                $element = $next_elem;
            }
        }
        %bins = map { $_ => 0 } @k;
    }
    else
    {
        return undef unless $partitions >= 1;
        my $interval = $self->sample_range() / $partitions;
        foreach my $idx (1 .. ($partitions-1))
        {
            push @k, ($self->min() + $idx * $interval);
        }

        $bins{$self->max()} = 0;

        push @k, $self->max();
    }

    ELEMENT:
    foreach my $element (@{$self->_data()})
    {
        foreach my $limit (@k)
        {
            if ($element <= $limit)
            {
                $bins{$limit}++;
                next ELEMENT;
            }
        }
    }

    return $self->_frequency(\%bins);
}

sub frequency_distribution {
    my $self = shift;

    my $ret = $self->frequency_distribution_ref(@_);

    if (!defined($ret))
    {
        return undef;
    }
    else
    {
        return %$ret;
    }
}

sub least_squares_fit {
  my $self = shift;
  return () if $self->count() < 2;

  ##Sigma sums
  my ($sigmaxy, $sigmax, $sigmaxx, $sigmayy, $sigmay) = (0,0,0,0,$self->sum);
  my ($xvar, $yvar, $err);

  ##Work variables
  my ($iter,$y,$x,$denom) = (0,0,0,0);
  my $count = $self->count();
  my @x;

  ##Outputs
  my ($m, $q, $r, $rms);

  if (!defined $_[1]) {
    @x = 1..$self->count();
  }
  else {
    @x = @_;
    if ( $self->count() != scalar @x) {
      carp "Range and domain are of unequal length.";
      return ();
    }
  }
  foreach $x (@x) {
    $y = $self->_data->[$iter];
    $sigmayy += $y * $y;
    $sigmaxx += $x * $x;
    $sigmaxy += $x * $y;
    $sigmax  += $x;
    $iter++;
  }
  $denom = $count * $sigmaxx - $sigmax*$sigmax;
  return ()
    unless abs( $denom ) > $Statistics::Descriptive::Tolerance;

  $m = ($count*$sigmaxy - $sigmax*$sigmay) / $denom;
  $q = ($sigmaxx*$sigmay - $sigmax*$sigmaxy ) / $denom;

  $xvar = $sigmaxx - $sigmax*$sigmax / $count;
  $yvar = $sigmayy - $sigmay*$sigmay / $count;

  $denom = sqrt( $xvar * $yvar );
  return () unless (abs( $denom ) > $Statistics::Descriptive::Tolerance);
  $r = ($sigmaxy - $sigmax*$sigmay / $count )/ $denom;

  $iter = 0;
  $rms = 0.0;
  foreach (@x) {
    ##Error = Real y - calculated y
    $err = $self->_data->[$iter] - ( $m * $_ + $q );
    $rms += $err*$err;
    $iter++;
  }

  $rms = sqrt($rms / $count);
  
  $self->_least_squares_fit([$q, $m, $r, $rms]);

  return @{ $self->_least_squares_fit() };
}

1;

package Statistics::Descriptive;

##All modules return true.
1;

__END__