File Coverage

lib/Statistics/Descriptive/PDL/Weighted.pm

Criterion	Covered	Total	%
statement	197	198	99.4
branch	47	62	75.8
condition	8	12	66.6
subroutine	30	30	100.0
pod	6	6	100.0
total	288	308	93.5

package Statistics::Descriptive::PDL::Weighted;

use 5.010;
use strict;
use warnings;

#  avoid loading too much, especially into our name space
use PDL::Lite '2.012';

#  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.

## no critic (ProhibitExplicitReturnUndef)

our $VERSION = '0.17';

use parent 'Statistics::Descriptive::PDL';

my @cache_methods = qw /
  count sum mode median
  mean standard_deviation skewness kurtosis
  geometric_mean harmonic_mean
  sum_weights sum_sqr_weights
/;
__PACKAGE__->_make_caching_accessors( \@cache_methods );


sub new {
    my $proto = shift;
    my $data_type = shift // PDL::double();
    my $class = ref($proto) || $proto;

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

    return $self;
}


sub _wt_type{PDL::double}

sub add_data {
    my ($self, $data, $weights) = @_;

    my ($data_piddle, $weights_piddle);

    my $data_from_hash;

    if (ref $data eq 'HASH') {
        $data_piddle    = PDL->pdl ($self->{data_type}, [keys %$data])->flat;
        $weights_piddle = PDL->pdl ($self->_wt_type, [values %$data])->flat;
        $data_from_hash = 1;
    }
    else {
        $data_piddle    = PDL->pdl ($self->{data_type}, $data)->flat;
        $weights_piddle = PDL->pdl ($self->_wt_type, $weights)->flat;
        die "data and weight vectors not of same length"
          if scalar $data_piddle->nelem != $weights_piddle->nelem;
        die "Cannot pass zero or negative weights"
          if PDL::any($weights_piddle <= 0);
    }

    return if !$data_piddle->nelem;

    my $has_existing_data = $self->count;

    # Take care of appending to an existing data set
    if ($has_existing_data) {
        my $d_piddle = $self->_get_piddle;
        $d_piddle    = $d_piddle->append ($data_piddle);
        $self->_set_piddle ($d_piddle);
        my $w_piddle = $self->_get_weights_piddle;
        $w_piddle    = $w_piddle->append ($weights_piddle);
        $self->_set_weights_piddle ($w_piddle);

        delete $self->{sorted};
    }
    else {
        $self->_set_piddle ($data_piddle);
        $self->_set_weights_piddle ($weights_piddle);
    }

    #  need to clear late because count is cached
    $self->clear_cache;

    #  somewhat awkward but needs to be set after clearing the cache
    if ($data_from_hash && !$has_existing_data) {
        $self->values_are_unique(1);
    }

    return $self->count;
}

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

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

    return wantarray ? ($data, $wts) : [$data, $wts];
}

sub get_data_as_hash {
    my $self = shift;

    my $piddle = $self->_get_piddle;

    return wantarray ? () : {}
      if !defined $piddle;

    my $deduped = $self->_deduplicate;

    my $data = $deduped->_get_piddle->unpdl;
    my $wts  = $deduped->_get_weights_piddle->unpdl;
    my %h;
    @h{@$data} = @$wts;

    return wantarray ? (%h) : \%h;
}


sub values_are_unique {
    my $self = shift;
    if (@_) {
        my $flag = shift;
        $self->{_cache}{deduplicated} = !!$flag;
    }
    return $self->{_cache}{deduplicated};
}

sub _set_weights_piddle {
    my ($self, $data) = @_;
    $self->{weights_piddle} = PDL->pdl ($data);
}

sub _get_weights_piddle {
    my $self = shift;
    return $self->{weights_piddle};
}

sub _count {
    my $self = shift;
    my $piddle = $self->_get_weights_piddle;
    return undef if !defined $piddle;
    return $piddle->sum;
}

sub _sum {
    my $self = shift;
    return ($self->_get_piddle * $self->_get_weights_piddle)->sum;
}

sub _sum_weights {
    my $self = shift;

    return $self->_get_weights_piddle->sum;
}

sub _sum_sqr_weights {
    my $self = shift;
    return $self->_get_weights_piddle->power(2)->sum;
}

#sub min_weight {
#    my $self = shift;
#    return $self->_get_weights_piddle->min;
#}


sub _mean {
    my $self = shift;

    my $data = $self->_get_piddle;
    my $wts  = $self->_get_weights_piddle;
    return ($data * $wts)->sum / $wts->sum;
}


sub _standard_deviation {
    my $self = shift;

    my $data = $self->_get_piddle;
    my $sd;
    my $n = $data->nelem;
    if ($n > 1) {
        #  long winded approach
        my $wts  = $self->_get_weights_piddle;
        my $mean = $self->mean;
        my $sumsqr = ($wts * (($data - $mean) ** 2))->sum;
        my $var = $sumsqr / $self->sum_weights;
        $sd = sqrt $var;
    }
    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;

    my $data = $self->_sort_piddle;
    my $cumsum = $self->_get_cumsum_weight_vector;

    my $target_wt = $self->sum_weights * 0.5;
    #  vsearch should be faster since it uses a binary search
    my $idx = PDL->pdl($target_wt)->vsearch_insert_leftmost($cumsum->reshape);

    return $data->at($idx);
}

sub _sort_piddle {
    my $self = shift;
    my $data = $self->_get_piddle;

    return undef if !defined $data;

    return $data if $self->{sorted};

    my $wts = $self->_get_weights_piddle;
    my $s = $data->qsorti->reshape;
    my $sorted_data = $data->slice($s);
    my $sorted_wts  = $wts->slice($s);

    $self->_set_piddle($sorted_data);
    $self->_set_weights_piddle($sorted_wts);

    $self->{sorted} = 1;
    delete $self->{_cache}{cumsum_weight_vector};

    return $sorted_data;
}

#  de-duplicate if needed, aggregating weights
#  there should be a sumover or which approach that will work better
#  maybe yvals related
sub _deduplicate {
    my ($self, %args) = @_;
    my $piddle = $self->_get_piddle;

    return undef
      if !defined $piddle;

    return $self
      if $self->values_are_unique;

    my $unique = $piddle->uniq;

    if ($unique->nelem == $piddle->nelem) {
        $self->values_are_unique(1);
        return $self
    }

    if (!$self->{sorted}) {
        $unique = $unique->qsort;
    }

    $piddle = $self->_sort_piddle;
    my $wts_piddle = $self->_get_weights_piddle;

    #  could use a map with a hash, but this avoids
    #  stringification and loss of precision
    #  (not that that should cause too many issues for most data)
    #  Should try to reduce looping when there are
    #  not many dups in large data sets
    my $j        = 0;  #  index into deduplicated piddle
    my $last_val = $piddle->at(0);
    my @wts;
    foreach my $i (0..$piddle->nelem-1) {
        my $val = $piddle->at($i);
        if ($val != $last_val) {
            $j++;
            $last_val = $val;
        }
        $wts[$j] += $wts_piddle->at($i);
    }

    if ($args{inplace}) {
        $self->_set_piddle($unique);
        $self->_set_weights_piddle(\@wts);
        $self->clear_cache;
        $self->values_are_unique (1);
        return $self;
    }

    my $new = $self->new($self->{data_type});
    $new->_set_piddle($unique);
    $new->_set_weights_piddle(\@wts);
    $new->values_are_unique (1);

    return $new;
}

sub _get_cumsum_weight_vector {
    my $self = shift;

    return $self->{_cache}{cumsum_weight_vector}
      if defined $self->{_cache}{cumsum_weight_vector};
    return $self->{_cache}{cumsum_weight_vector}
      = $self->_get_weights_piddle->cumusumover->reshape;
}

sub _skewness {
    my $self = shift;

    my $data = $self->_get_piddle;

    #  long winded approach
    my $mean = $self->mean;
    my $sd   = $self->standard_deviation;
    my $wts = $self->_get_weights_piddle;
    my $sumpow3 = ($wts * ((($data - $mean) / $sd) ** 3))->sum;
    my $skew = $sumpow3 / $self->sum_weights;
    return $skew;
}

sub _kurtosis {
    my $self = shift;

    my $data = $self->_get_piddle;

    #  long winded approach
    my $mean = $self->mean;
    my $sd   = $self->standard_deviation;
    my $wts = $self->_get_weights_piddle;
    my $sumpow4 = ($wts * ((($data - $mean) / $sd) ** 4))->sum;
    my $kurt = $sumpow4 / $self->sum_weights - 3;
    return $kurt;
}


sub _harmonic_mean {
    my $self = shift;

    my $data = $self->_get_piddle;

    #  not sure about this...
    return undef if $data->which->nelem != $data->nelem;

    my $wts = $self->_get_weights_piddle;

    my $hs = ((1 / $data) * $wts)->sum;

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

sub _geometric_mean {
    my $self = shift;

    my $data = $self->_get_piddle;

    #  should add a sorted status check, as we can use vsearch in such cases
    return undef if $data->where($data < 0)->nelem;

    my $wts = $self->_get_weights_piddle;

    # formula from https://en.wikipedia.org/wiki/Weighted_geometric_mean
    return exp (($data->log * $wts)->sum / $wts->sum);
}

sub _mode {
    my $self = shift;

    #  de-duplicate and aggregate weights if needed
    my $deduped = $self->_deduplicate;

    my $data = $deduped->_get_piddle;


    my $wts = $deduped->_get_weights_piddle;
    my $mode = $data->at($wts->maximum_ind);
    if ($mode > $data->max) {
        #  PDL returns strange numbers when distributions are flat
        $mode = undef;
    }

    return $mode;
}

#  need to convert $p to fraction, or perhaps die if it is between 0 and 1
sub _percentile {
    my ($self, $p) = @_;

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

    $data = $self->_sort_piddle;

    my $cumsum = $self->_get_cumsum_weight_vector;

    my $target_wt = $self->sum_weights * ($p / 100);

    my $idx = PDL->pdl($target_wt)->vsearch_insert_leftmost($cumsum->reshape);

    return $data->at($idx);
}


1;

__END__