File Coverage

blib/lib/Statistics/Descriptive/PDL.pm

Criterion	Covered	Total	%
statement	167	184	90.7
branch	49	68	72.0
condition	16	31	51.6
subroutine	40	41	97.5
pod	14	14	100.0
total	286	338	84.6

package Statistics::Descriptive::PDL;

## no critic (ProhibitExplicitReturnUndef)

use 5.010;
use strict;
use warnings;
use Scalar::Util qw /blessed/;
use POSIX qw /fmod/;

#  avoid loading too much, especially into our name space
use PDL::Lite '2.012';
#  try to keep running if PDL::Stats is not installed
eval 'require PDL::Stats::Basic';
my $has_PDL_stats_basic = $@ ? undef : 1;
#$has_PDL_stats_basic = 0;

#  We could inherit from PDL::Objects, but in this case we want
#  to hide the piddle from the caller to avoid arbitrary changes
#  being applied to it.

our $VERSION = '0.17';

our $Tolerance = 0.0;  #  for compatibility with Stats::Descr, but not used here

my @cache_methods = qw /
  count sum mode median
  mean standard_deviation skewness kurtosis
  geometric_mean harmonic_mean
  max min sample_range
  iqr
  sum_sqr_sample_weights
/;
__PACKAGE__->_make_caching_accessors( \@cache_methods );

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

    ## no critic
    no strict 'refs';
    ## use critic
    foreach my $method (@$methods)
    {
        *{ $pkg . "::" . $method } = do
        {
            my $m = $method;
            sub {
                my $self = shift;
                return $self->{_cache}{$method}
                  if defined $self->{_cache}{$method};

                my $piddle = $self->_get_piddle;
                return undef
                  if !defined $piddle or $piddle->isempty;

                my $call_meth = "_$method";
                my $val = $self->$call_meth;

                if (blessed $val and $val->isa('PDL')) {
                    $val = $val->sclr;
                }
                return $self->{_cache}{$method} = $val;
            };
        };
    }

    return;
}

sub new {
    my $proto = shift;
    my $class = ref($proto) || $proto;

    my $self = {piddle => undef};
    bless $self, $class;

    return $self;
}

sub available_stats {
    my $self = shift;
    my @methods = sort (@cache_methods, qw /percentile variance/);
    return wantarray ? @methods : \@methods;
}

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

    #  have we been passed an ndarray?
    if (blessed $_[0] && $_[0]->isa('PDL')) {
        $data = $_[0];
    }
    else {
        $data
          = ref ($_[0]) eq 'ARRAY'
          ? $_[0]
          : \@_;
        return if !scalar @$data;
    }

    my $piddle;
    my $has_existing_data = $self->count;

    $self->clear_cache;

    # Take care of appending to an existing data set
    if ($has_existing_data) {
        $piddle = $self->_get_piddle;
        $piddle = $piddle->append (PDL->pdl ($data)->flat);
        $self->_set_piddle ($piddle);
    }
    else {
        $self->_set_piddle (PDL->pdl($data)->flat);
    }

    return $self->count;
}

sub get_data {
    my $self = shift;
    my $piddle = $self->_get_piddle;

    my $data = defined $piddle ? $piddle->unpdl : [];

    return wantarray ? @$data : $data;
}

sub get_data_as_hash {
    my $self = shift;

    my $piddle = $self->_get_piddle;
    if (defined $piddle) {
        require Statistics::Descriptive::PDL::SampleWeighted;
        my $wtd_obj = Statistics::Descriptive::PDL::SampleWeighted->new;
        my $wts_piddle = PDL->ones ($piddle->dims);
        $wtd_obj->add_data ($piddle->copy, $wts_piddle);
        return $wtd_obj->get_data_as_hash;
    }

    return wantarray ? () : {};
}

sub values_are_unique {}

#  flatten $data if multidimensional
sub _set_piddle {
    my ($self, $data) = @_;
    $self->{piddle} = PDL->pdl ($data);
}

sub _get_piddle {
    my $self = shift;
    return $self->{piddle};
}

sub clear_cache {
    my $self = shift;
    delete $self->{_cache};
    return;
}

sub _count {
    my $self = shift;
    return $self->_get_piddle->nelem;
}

sub _sum {
    my $self = shift;
    return $self->_get_piddle->sum;
}

sub sum_weights {
    my $self = shift;
    return $self->count;
}

sub sum_sqr_weights {
    my $self = shift;
    return $self->count;
}

sub _sum_sqr_sample_weights {
    my $self = shift;
    return $self->sum_sqr_weights;
}

sub _min {
    my $self = shift;
    return $self->_get_piddle->min;
}

sub _max {
    my $self = shift;
    return $self->_get_piddle->max;
}

sub _mean {
    my $self = shift;
    return $self->_get_piddle->average;
}


sub sd    {return $_[0]->standard_deviation}
sub stdev {return $_[0]->standard_deviation}

sub _standard_deviation {
    my $self = shift;

    my $piddle = $self->_get_piddle;

    my $sd;
    my $n = $piddle->nelem;
    if ($n > 1) {
        if ($has_PDL_stats_basic) {
            $sd = $piddle->stdv_unbiased;
        }
        else {
            my $var = (($piddle ** 2)->sum - $n * $self->mean ** 2);
            $sd = $var > 0 ? sqrt ($var / ($n - 1)) : 0;
        }
    }
    elsif ($n == 1){
        $sd = 0;
    }
    return $sd;
}

sub variance {
    my $self = shift;
    my $sd = $self->standard_deviation;
    return defined $sd ? $sd ** 2 : undef;
}

sub _median {
    my $self = shift;
    return $self->_get_piddle->median;
}


sub _skewness {
    my $self = shift;

    my $piddle = $self->_get_piddle;

    my $n = $piddle->nelem;

    return undef if $n < 3;

    return $piddle->skew_unbiased
      if $has_PDL_stats_basic;

    #  do it ourselves
    my $mean = $self->mean;
    my $sd   = $self->standard_deviation;
    my $sumpow3 = ((($piddle - $mean) / $sd) ** 3)->sum;
    my $correction = $n / ( ($n-1) * ($n-2) );
    my $skew = $correction * $sumpow3;

    return $skew;
}

sub _kurtosis {
    my $self = shift;
    my $piddle = $self->_get_piddle;

    my $n = $piddle->nelem;

    return undef if $n < 4;

    return $piddle->kurt_unbiased
      if $has_PDL_stats_basic;

    #  do it ourselves
    my $mean = $self->mean;
    my $sd   = $self->standard_deviation;
    my $sumpow4 = ((($piddle - $mean) / $sd) ** 4)->sum;

    my $correction1 = ( $n * ($n+1) ) / ( ($n-1) * ($n-2) * ($n-3) );
    my $correction2 = ( 3  * ($n-1) ** 2) / ( ($n-2) * ($n-3) );

    my $kurt = ( $correction1 * $sumpow4 ) - $correction2;
    return $kurt;
}

sub _sample_range {
    my $self = shift;
    my $min = $self->min // return undef;
    my $max = $self->max // return undef;
    return $max - $min;
}


sub _harmonic_mean {
    my $self = shift;
    my $piddle = $self->_get_piddle;

    return undef if $piddle->which->nelem != $piddle->nelem;

    my $hs = (1 / $piddle)->sum;

    return $hs ? $self->count / $hs : undef;
}

sub _geometric_mean {
    my $self = shift;
    my $piddle = $self->_get_piddle;

    my $count = $self->count;

    return undef if !$count;
    return undef if $piddle->where($piddle < 0)->nelem;

    my $exponent = 1 / $self->count();
    my $powered = $piddle ** $exponent;

    my $gm = $powered->dprodover;
}

sub _mode {
    my $self = shift;
    my $piddle = $self->_get_piddle;

    my $count  = $piddle->nelem;
    my $unique = $piddle->uniq;

    return undef if $unique->nelem == $count or $unique->nelem == 1;

    #if (!($count % $unique->nelem)) {
    #    #  might have equal numbers of each value
    #    #  need to check for this, but for now return undef
    #    return undef;
    #}

    my $mode = $piddle->mode;

    #  bodge to handle odd values
    return undef if !$piddle->in($mode)->max;

    return $mode;
}

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

    my $piddle = $self->_get_piddle;

    return
      if !defined $piddle || $piddle->nelem == 0;

    my @vals = map {$self->percentile($_)} @percentiles;

    return @vals;
}


#  caching wrapper
#  need to convert $p to fraction, or perhaps die if it is between 0 and 1
#  hard-coded cache percentiles not ideal
sub percentile {
    my ($self, $p) = @_;
    my $piddle = $self->_get_piddle;

    return undef
      if !defined $piddle || $piddle->nelem == 0;

    return $self->median
      if $p == 50;

    die "Percentile $p outside range 0..100"
      if $p < 0 or $p > 100;

    #  allow for other number formats like '005'
    #  needed for cache
    $p += 0;

    return $self->{_cache}{percentile}{$p}
      if defined $self->{_cache}{percentile}{$p};

    my $pctl = $self->_percentile($p);

    if (blessed $pctl && $pctl->isa('PDL')) {
        $pctl = $pctl->sclr;
    }

    if (int ($p) == $p) {
        $self->{_cache}{percentile}{$p} = $pctl;
    }

    return $pctl;
}

sub _percentile {
    my ($self, $p) = @_;

    return $self->_get_piddle->pct($p / 100);
}



sub _iqr {
    my $self = shift;
    $self->percentile(75) - $self->percentile(25);
}


1;

__END__