File Coverage

lib/Statistics/Descriptive/PDL/SampleWeighted.pm

Criterion	Covered	Total	%
statement	75	77	97.4
branch	11	16	68.7
condition	4	9	44.4
subroutine	13	13	100.0
pod			n/a
total	103	115	89.5

package Statistics::Descriptive::PDL::SampleWeighted;

use 5.010;
use strict;
use warnings;

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

#  this is otherwise not loaded due to oddities with multiple loading of PDL::Lite
#*pdl = \&PDL::Core::pdl;

#  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::Weighted';

my @cache_methods = qw /
  median
  standard_deviation skewness kurtosis
/;
__PACKAGE__->_make_caching_accessors( \@cache_methods );



sub _wt_type{PDL::long}


sub _standard_deviation {
    my $self = shift;

    my $data = $self->_get_piddle;

    my $wts = $self->_get_weights_piddle;
    my $n   = $wts->sum;

    return 0 if $n == 1;

    my $var = ((($data ** 2) * $wts)->sum - $n * $self->mean ** 2);

    return $var > 0 ? sqrt ($var / ($n - 1)) : 0;
}


sub _median {
    my $self = shift;

    my $piddle = $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);

    #  if the target weight is "on a boundary" between
    #  two values then we need to interpolate
    my $median
      = $target_wt == $cumsum->at($idx)
      ? ($piddle->at($idx) + $piddle->at($idx+1)) / 2
      : $piddle->at($idx);

    return $median;
}


sub _skewness {
    my $self = shift;

    my $n = $self->sum_weights;
    return undef if $n < 3;

    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;
    #  inplace seems not to be faster here.
    #  Possibly PDL is smart enough to do it by default
    #  in such cases
    #my $sumpow3 = ($data - $mean)->inplace->divide($sd, 0)->pow(3)->mult($wts, 0)->sum;
    my $correction = $n / ( ($n-1) * ($n-2) );
    my $skew = $correction * $sumpow3;
    return $skew;
}

sub _kurtosis {
    my $self = shift;

    my $n    = $self->sum_weights;
    return undef if $n <= 3;

    my $data = $self->_get_piddle;
    my $mean = $self->mean;
    my $sd   = $self->standard_deviation;
    my $wts  = $self->_get_weights_piddle;

    my $sumpow4 = ($wts * ((($data - $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) );

    return ( $correction1 * $sumpow4 ) - $correction2;
}

#  crude memoisation - would be nice to use
#  state but it has issues with lists on older perls
my %k_piddle_cache;

#  Uses same basic algorithm as PDL::pctl.
sub _percentile {
    my ($self, $p) = @_;

    my $piddle = $self->_get_piddle;

    return undef
      if !defined $piddle or $piddle->isempty;

    $self->_sort_piddle;

    #  possible slowdown here - users need to dedup before calling to avoid
    # my $deduped = $self->_deduplicate;
    #  there is actually no need to dedup
    my $deduped = $self;
    $piddle = $deduped->_get_piddle;

    # my $wt_piddle = $self->_get_weights_piddle;

    my $cumsum = $deduped->_get_cumsum_weight_vector;
    my $wt_sum = $deduped->sum_weights;

    use POSIX qw /floor/;

    my $target_wt = ($p / 100) * ($wt_sum - 1) + 1;
    my $k = floor $target_wt;
    my $d = $target_wt - $k;

    my $idx = ($k_piddle_cache{$k} //= PDL->pdl(PDL::indx(), [$k]))->vsearch_insert_leftmost($cumsum)->at(0);

    if (scalar keys %k_piddle_cache > 10000) {
        #  Reset if we get too many
        #  - could be more nuanced based on frequency
        #  but then we would have to track it
        %k_piddle_cache = ();
    }

    #  we need to interpolate if our target weight falls between two sets of weights
    #  e.g. target is 1.3, but the cumulative weights are [1,2] or [1,5]
    my $fraction = $target_wt - ($cumsum->at($idx));
    if ($fraction > 0 && $fraction < 1) {
        my $lower_val = $piddle->at($idx);
        my $upper_val = $piddle->at($idx+1);
        my $val = $lower_val + $d * ($upper_val - $lower_val);
        return $val;
    }

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

#  weight for each sample is 1
sub _sum_sqr_sample_weights {
    my $self = shift;
    return $self->sum_weights;
}



1;

__END__