File Coverage

blib/lib/Statistics/Descriptive/Full.pm

Criterion	Covered	Total	%
statement	395	414	95.4
branch	97	118	82.2
condition	18	23	78.2
subroutine	43	45	95.5
pod	26	26	100.0
total	579	626	92.4

package Statistics::Descriptive::Full;
$Statistics::Descriptive::Full::VERSION = '3.0702';
use strict;
use warnings;

use Carp;
use POSIX ();
use Statistics::Descriptive::Smoother;

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

use parent qw(Statistics::Descriptive::Sparse);

use List::MoreUtils ();
use List::Util ();

##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
  samples    => undef,  ##Number of samples for each value of the data set
  presorted  => undef,  ##Flag to indicate the data is already sorted
  _reserved  => undef,  ##Place holder for this lookup hash
);

__PACKAGE__->_make_private_accessors(
    [qw(data samples frequency geometric_mean harmonic_mean
        least_squares_fit median mode
        skewness kurtosis median_absolute_deviation
       )
    ]
);
__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->_samples([]);
    $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;

    my %keys = %{ $self };

    # Remove reserved keys for this class from the deletion list
    delete @keys{keys %{$self->_reserved}};
    delete @keys{keys %{$self->_permitted}};
    delete $keys{_trimmed_mean_cache};

    KEYS_LOOP:
    foreach my $key (keys %keys) { # Check each key in the object
        delete $self->{$key};  # Delete any out of date cached key
    }
    $self->{_trimmed_mean_cache} = {};  #  just reset this one
    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;  ##Myself

    my $aref;

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

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

    my $oldmean;
    my ($min, $max, $sum, $sumsq);
    my $count = $self->count;

    #  $count is modified lower down, but we need this flag after that
    my $has_existing_data = $count;

    # Take care of appending to an existing data set
    if ($has_existing_data) {
        $min   = $self->min();
        $max   = $self->max();
        $sum   = $self->sum();
        $sumsq = $self->sumsq();
    }
    else {
        $min   = $aref->[0];
        $max   = $aref->[0];
        $sum   = 0;
        $sumsq = 0;
    }

    #  need to allow for already having data
    $sum    += List::Util::sum (@$aref);
    $sumsq  += List::Util::sum (map {$_ ** 2} @$aref);
    $max    =  List::Util::max ($max, @$aref);
    $min    =  List::Util::min ($min, @$aref);
    $count  +=  scalar @$aref;
    my $mean = $sum / $count;

    $self->min($min);
    $self->max($max);
    $self->sample_range($max - $min);
    $self->sum($sum);
    $self->sumsq($sumsq);
    $self->mean($mean);
    $self->count($count);

    ##Variance isn't commonly enough
    ##used to recompute every single data add, so just clear its cache.
    $self->_variance(undef);

    push @{ $self->_data() }, @{ $aref };

    #  no need to clear keys if we are a newly populated object,
    #  and profiling shows it takes a long time when creating
    #  and populating many stats objects
    if ($has_existing_data) {
        ##Clear the presorted flag
        $self->presorted(0);
        $self->_delete_all_cached_keys();
    }

    return 1;
}


sub add_data_with_samples {
    my ($self,$aref_values) = @_;

    return 1 if (!@{ $aref_values });

    my $aref_data = [map { keys %$_ } @{ $aref_values }];
    my $aref_samples = [map { values %$_ } @{ $aref_values }];

    $self->add_data($aref_data);
    push @{ $self->_samples() }, @{ $aref_samples };

    return 1;
}


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

sub get_data_without_outliers {
    my $self = shift;

    if ($self->count() < $Statistics::Descriptive::Min_samples_number) {
        carp("Need at least $Statistics::Descriptive::Min_samples_number samples\n");
        return;
    }

    if (!defined $self->{_outlier_filter}) {
        carp("Outliers filter not defined\n");
        return;
    }

    my $outlier_candidate_index = $self->_outlier_candidate_index;
    my $possible_outlier = ($self->_data())->[$outlier_candidate_index];
    my $is_outlier = $self->{_outlier_filter}->($self, $possible_outlier);

    return $self->get_data unless $is_outlier;
    # Removing the outlier from the dataset
    my @good_indexes = grep { $_ != $outlier_candidate_index } (0 .. $self->count() - 1);

    my @data = $self->get_data;
    my @filtered_data = @data[@good_indexes];
    return @filtered_data;
}

sub set_outlier_filter {
    my ($self, $code_ref) = @_;

    if (!$code_ref || ref($code_ref) ne "CODE") {
        carp("Need to pass a code reference");
        return;
    }

    $self->{_outlier_filter} = $code_ref;
    return 1;
}

sub _outlier_candidate_index {
    my $self = shift;

    my $mean = $self->mean();
    my $outlier_candidate_index = 0;
    my $max_std_deviation = abs(($self->_data())->[0] - $mean);
    foreach my $idx (1 .. ($self->count() - 1) ) {
        my $curr_value = ($self->_data())->[$idx];
        if ($max_std_deviation  <  abs($curr_value - $mean) ) {
            $outlier_candidate_index = $idx;
            $max_std_deviation = abs($curr_value - $mean);
        }
    }
    return $outlier_candidate_index;
}

sub set_smoother {
    my ($self, $args) = @_;

    $args->{data}    = $self->_data();
    $args->{samples} = $self->_samples();

    $self->{_smoother} = Statistics::Descriptive::Smoother->instantiate($args);
}

sub get_smoothed_data {
    my ($self, $args) = @_;

    if (!defined $self->{_smoother}) {
        carp("Smoother object not defined\n");
        return;
    }
    $self->{_smoother}->get_smoothed_data();
}

sub maxdex {
    my $self = shift;

    return undef if !$self->count;
    my $maxdex;

    if ($self->presorted) {
        $maxdex = $self->count - 1;
    }
    else {
        my $max = $self->max;
        $maxdex =  List::MoreUtils::first_index {$_ == $max} $self->get_data;
    }

    $self->{maxdex} = $maxdex;

    return $maxdex;
}

sub mindex {
    my $self = shift;

    return undef if !$self->count;
    #my $maxdex = $self->{maxdex};
    #return $maxdex if defined $maxdex;
    my $mindex;

    if ($self->presorted) {
        $mindex = 0;
    }
    else {
        my $min = $self->min;
        $mindex = List::MoreUtils::first_index {$_ == $min} $self->get_data;
    }

    $self->{mindex} = $mindex;

    return $mindex;
}

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 - no, this is unnecessary now we have methods
      #$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();

    if ((! $count) || ($percentile < 100 / $count))
    {
        return;  #  allow for both scalar and list context
    }

    $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;

    return undef if !$self->count;

    ##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;
    }

    #  check data count after the args are checked - should help debugging
    return undef if !$self->count;

    $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]);
    }

    #  check data count after the args
    return undef if !$self->count;

    ##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 %count;

        foreach my $item (@{ $self->_data() })
        {
            my $count = ++$count{$item};
            if ($count > $occurances)
            {
                $mode = $item;
                $occurances = $count;
            }
        }

        $self->_mode(
            ($occurances > 1)
            ? {exists => 1, mode => $mode}
            : {exists => 0,}
        );
    }

    my $m = $self->_mode;

    return $m->{'exists'} ? $m->{mode} : undef;
}

sub geometric_mean {
    my $self = shift;

    return undef if !$self->count;

    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 ) {
                $sum_pow3 += (($rec - $mean) / $sd) ** 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() };
}

sub median_absolute_deviation {
    my ($self) = @_;

    if (!defined($self->_median_absolute_deviation()))
    {
        my $stat = $self->new;
        $stat->add_data(map { abs($_ - $self->median) } $self->get_data);
        $self->_median_absolute_deviation($stat->median);
    }

    return $self->_median_absolute_deviation();
}

sub summary {
    my ($self) = @_;

    my $FMT = '%.5e';

    return sprintf("Min: $FMT\nMax: $FMT\nMean: $FMT\nMedian: $FMT\n" .
        "1st quantile: $FMT\n3rd quantile: $FMT\n",
        $self->min,
        $self->max,
        $self->mean,
        $self->median,
        $self->quantile(1),
        $self->quantile(3),
    );

}
1;

__END__