File Coverage

blib/lib/Statistics/Descriptive/Weighted.pm

Criterion	Covered	Total	%
statement	234	264	88.6
branch	70	88	79.5
condition	24	40	60.0
subroutine	27	36	75.0
pod			n/a
total	355	428	82.9

package Statistics::Descriptive::Weighted;
$VERSION = '0.7';
use Statistics::Descriptive;
use Data::Dumper;

package Statistics::Descriptive::Weighted::Sparse;
use strict;
use vars qw($AUTOLOAD @ISA %fields);

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

use Carp qw(cluck confess);

##Define a new field to be used as method, to
##augment the ones inherited
%fields = (
           mean                      => 0,
	   weight                    => 0,
	   sum_squares               => 0,
	   weight_homozyg            => 0,
	   biased_variance           => 0,
	   biased_standard_deviation => 0
  );

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

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

sub add_data {
  my $self = shift;  ##Myself
  my $oldmean;
  my $oldweight;
  my ($min,$max);
  my $aref;

  if ( (not ref $_[0] eq 'ARRAY') || (exists $_[1] and (not (ref $_[1] eq 'ARRAY') || @{$_[0]} != @{$_[1]} ) ) ) {
    cluck "WARNING: Expected input are two references to two arrays of equal length; first data, then positive weights. Second array is optional.\n";
    return undef;
  }

  my ($datum,$weight) = @_;

  ##Calculate new mean, pseudo-variance, min and max;
  ##The on-line weighted incremental algorithm for variance is based on West 1979 from Wikipedia
  ##D. H. D. West (1979). Communications of the ACM, 22, 9, 532-535: Updating Mean and Variance Estimates: An Improved Method

  ## NEW in Version 0.4:
  ## I calculate a sample weighted variance based on normalized weights rather than the sample size
  ## correction factor is: 1 / (1 - sum [w_i / (sum w_i) ]^2)

  ## call H = sum [w_i / (sum w_i) ]^2. An online update eq for H is  H_new =  (sum.w_old^2 * H_old) + weight^2) / sum.w^2

  ## correction factor is then 1 / (1 - H_new)

  my $weighterror;
  for (0..$#$datum ) {
    if (not defined $$weight[$_]) {
      $$weight[$_] = 1;
    }
    if ($$weight[$_] <= 0) {
      $weighterror = 1;
      next;
    }
    $oldmean = $self->{mean};
    $oldweight = $self->{weight};
    $self->{weight} += $$weight[$_];
    $self->{weight_homozyg} = ((($oldweight ** 2 * $self->{weight_homozyg}) + $$weight[$_] ** 2) / ( $self->{weight} ** 2 ));
    $self->{count}++;
    $self->{sum} += ($$weight[$_] * $$datum[$_]);
    $self->{mean} += (($$weight[$_] / $self->{weight} ) * ($$datum[$_] - $oldmean));
    $self->{sum_squares} += (($$weight[$_] / $self->{weight} ) * ($$datum[$_] - $oldmean) ** 2) * $oldweight;
    if (not defined $self->{max} or $$datum[$_] > $self->{max}) {
      $self->{max} = $$datum[$_];
    }
    if (not defined $self->{min} or $$datum[$_] < $self->{min}) {
      $self->{min} = $$datum[$_];
    }
  }
  cluck "WARNING: One or more data with nonpositive weights were skipped.\n" if ($weighterror);
    $self->{sample_range} = $self->{max} - $self->{min};

  if ($self->{count} > 1) {
    $self->{variance}     = ($self->{sum_squares} / ((1 - $self->{weight_homozyg}) * $self->{weight}));
    $self->{standard_deviation}  = sqrt( $self->{variance});
    $self->{biased_variance}     = ($self->{sum_squares} / $self->{weight});
    $self->{biased_standard_deviation}  = sqrt( $self->{biased_variance});
  }
  return 1;
}

sub weight {
  my $self = shift;
  if (@_ > 0) {
    cluck "WARNING: Sparse statistics object expects zero arguments to weight function, returns sum of weights.";
  }
  return $self->_weight();
}

## OVERRIDES FOR UNSUPPORTED FUNCTIONS

sub mindex{
  confess "ERROR: Statistics::Descriptive::Weighted does not support this function.";
}

sub maxdex{
  confess "ERROR: Statistics::Descriptive::Weighted does not support this function.";
}

1;

package Statistics::Descriptive::Weighted::Full;

use Carp qw(cluck confess);
use Tree::Treap;
use strict;
use vars qw(@ISA %fields);

@ISA = qw(Statistics::Descriptive::Weighted::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,  ##keys from variate values to a hashref with keys weight, cdf, tail-prob
  did_cdf    => undef,  ##flag to indicate whether CDF/quantile fun has been computed or not
  quantile   => undef,  ##"hash" for quantile function
  percentile => undef,  ##"hash" for percentile function
  maxweight  => 0,
  mode       => undef,
  order      => 1,
  _reserved  => undef,  ##Place holder for this lookup hash
);

##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;
  my $self = $class->SUPER::new();  ##Create my self re SUPER
  $self->{data}           = new Tree::Treap("num");  ## inserts data by numeric comparison
  $self->{did_cdf}        = 0;
  $self->{maxweight}      = 0;
  $self->{quantile}       = new Tree::Treap("num");
  $self->{percentile}     = new Tree::Treap("num");
  $self->{order}          = 1;
  $self->{'_reserved'} = \%fields;
  bless ($self, $class);
  return $self;
}

## The treap gives relatively fast search and good performance on possibly sorted data
## The choice is motivated by heavy intended use for Empirical Distribution Function
## A lot of work is done at insertion for faster computation on search
## THE ACTUAL DATA INSERTION IS DONE AT FUNCTION _addweight
## The data structure loses information. Like a hash keys appear only once.
## The value of a key is its sum of weight for that key, and the cumulative weight
sub add_data {
  my $self = shift;
  my $key;

  if ( (not ref $_[0] eq 'ARRAY') || (exists $_[1] and (not (ref $_[1] eq 'ARRAY') || @{$_[0]} != @{$_[1]} ) ) ) {
    cluck "WARNING: Expected input are two references to two arrays of equal length; first data, then positive weights. Second array is optional.\n";
    return undef;
  }

  my ($datum,$weight) = @_;
  my $filterdatum = [];
  my $filterweight = [];
  my $weighterror;
  my $newweight;
  for (0..$#$datum) {
    if (not defined $$weight[$_]) {
      $$weight[$_] = 1;
    }
    if ($$weight[$_] > 0) {
      push @$filterdatum,$$datum[$_];
      push @$filterweight,$$weight[$_];
      $newweight = $self->_addweight($$datum[$_], $$weight[$_]);
      if ($newweight > $self->{maxweight}) {
	$self->{maxweight} = $newweight;
	$self->{mode} = $$datum[$_];
      }
    }
    else {
      $weighterror = 1;
    }
  }
  cluck "WARNING: One or more data with nonpositive weights were skipped.\n" if ($weighterror);
  $self->SUPER::add_data($filterdatum,$filterweight);  ##Perform base statistics on the data
  ##Clear the did_cdf flag
  $self->{did_cdf} = 0;
  ##Need to delete all cached keys
  foreach $key (keys %{ $self }) { # Check each key in the object
    # If it's a reserved key for this class, keep it
    next if exists $self->{'_reserved'}->{$key};
    # If it comes from the base class, keep it
    next if exists $self->{'_permitted'}->{$key};
    delete $self->{$key};          # Delete the out of date cached key
  }
  return 1;
}

sub count {
  my $self = shift;
  if (@_ == 1) {  ##Inquire
    my $val = $self->{data}->get_val($_[0]);
    return (defined $val ? ${ $val }{'count'} : $val);
  }
  elsif (@_ == 0) {  ##Inquire
    return $self->{count};
  }
  else {
    cluck "WARNING: Only 1 or fewer arguments expected.";
  }
  return 1;
}

sub weight {
  my $self = shift;
  if (@_ == 1) {  ##Inquire
    my $val = $self->{data}->get_val($_[0]);
    return (defined $val ? ${ $val }{'weight'} : $val);
  }
  elsif (@_ == 0) {  ##Inquire
    return $self->{weight};
  }
  else {
    cluck "WARNING: Only 1 or fewer arguments expected.";
  }
  return 1;
}

sub _addweight {
  my $self = shift;
  my $oldweight = ($self->weight($_[0]) || 0);
  my $newweight = $_[1] + $oldweight;
  my $value = $self->{data}->get_val($_[0]);
  my $weights = ($value ? $$value{'weights'} : [] );
  push @$weights, $_[1];
  my $orders  = ($value ? $$value{'order'} : [] );
  push @$orders, $self->{order}++;
  my $newcount  = ($self->count($_[0]) || 0) + 1;
  if (@_ == 2) {  ##Assign
    my $values = {'weight' => $newweight, 'weights' => $weights, 'count' => $newcount, 'order' => $orders, 'cdf' => undef, 'rt_tail_prob' => undef, 'percentile' => undef};
    $self->{data}->insert($_[0],$values);
  }
  else {
    cluck "WARNING: Only two arguments (key, addend) expected.";
  }
  return $newweight;
}

sub _do_cdf {
  my $self = shift;
  my $cumweight = 0;
  foreach my $key ($self->{data}->keys()){
    my $value = $self->{data}->get_val($key);
    my $keyweight = $self->weight($key);

    my $oldcumweight = $cumweight;
    $cumweight += $keyweight;

    my $propcumweight = $cumweight / $self->{weight};
    my $right_tail_prob = (1 - ($oldcumweight /  $self->{weight}));
    my $percentile = ((100 / $self->{weight}) * ($cumweight - ($keyweight / 2)));

    $$value{'cdf'} = $propcumweight;
    $$value{'rt_tail_prob'} = $right_tail_prob;
    $$value{'percentile'} = $percentile;

    $self->{data}->insert($key,$value);
    $self->{quantile}->insert($propcumweight,$key);
    $self->{percentile}->insert($percentile,$key);
  }
  $self->{did_cdf} = 1;
  return 1;
}

sub quantile {
  my $self = shift;
  $self->_do_cdf() unless $self->{did_cdf};
  if (@_ == 1) {  ##Inquire
    my $proportion = shift;
    cluck "WARNING: Expects an argument between 0 and 1 inclusive." if ($proportion < 0 or $proportion > 1);
    my @keys = $self->{quantile}->range_keys($proportion, undef);
    my $key = $keys[0]; ## GET THE SMALLEST QUANTILE g.e. $proportion
    return $self->{quantile}->get_val($keys[0]);
  }
  else {
    cluck "WARNING: Exactly 1 argument expected.";
    return undef;
  }
}

sub percentile {
  my $self = shift;
  $self->_do_cdf() unless $self->{did_cdf};
  if (@_ != 1) {
    cluck "WARNING: Exactly 1 argument expected.";
  }
  my $percent = shift;
  if ($percent < 0 or $percent > 100) {
    cluck "WARNING: Expects an argument between 0 and 100 inclusive.";
  }
  my $percentile;
  if ($percent < $self->{percentile}->minimum()) {
    $percentile = $self->{data}->minimum();
  } elsif ($percent > $self->{percentile}->maximum()) {
    $percentile = $self->{data}->maximum();
  } else {
    my @lekeys =  $self->{percentile}->range_keys(undef,$percent);
    my $lekey = $lekeys[-1];
    my @gekeys =  $self->{percentile}->range_keys($percent, undef);
    my $gekey = $gekeys[0];
    my $leval = $self->{percentile}->get_val($lekey);
    $percentile = $leval;
    if ($gekey != $lekey) {
      my $geval = $self->{percentile}->get_val($gekey);
      $percentile += ($percent - $lekey) / ($gekey - $lekey) * ($geval - $leval);
    }
  }
  return $percentile;
}


sub median {
  my $self = shift;

  ##Cached?
  return $self->{median} if defined $self->{median};
  return $self->{median} = $self->percentile(50);
}

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

sub cdf {
  my $self = shift;
  $self->_do_cdf() unless $self->{did_cdf};
  if (@_ == 1) {  ##Inquire
    my $value = shift;
    return 0 if ($self->{data}->minimum() > $value);
    my @keys = $self->{data}->range_keys(undef, $value);
    my $key = $keys[-1]; ## GET THE LARGEST OBSERVED VALUE l.e. $value
    return ${ $self->{data}->get_val($key) }{'cdf'};
  }
  else {
    cluck "WARNING: Exactly 1 argument expected.";
    return undef;
  }
}

sub survival {
  my $self = shift;
  $self->_do_cdf() unless $self->{did_cdf};
  if (@_ == 1) {  ##Inquire
    my $value = shift;
    return 1 if ($self->{data}->minimum() > $value);
    my @keys = $self->{data}->range_keys(undef, $value);
    my $key = $keys[-1]; ## GET THE LARGEST OBSERVED VALUE l.e. $value
    return 1 - (${ $self->{data}->get_val($key) }{'cdf'});
  }
  else {
    cluck "WARNING: Only 1 argument expected.";
    return undef;
  }
}

sub rtp {
  my $self = shift;
  $self->_do_cdf() unless $self->{did_cdf};
  if (@_ == 1) {  ##Inquire
    my $value = shift;
    return 0 if ($self->{data}->maximum() < $value);
    my @keys = $self->{data}->range_keys($value, undef);
    my $key = $keys[0]; ## GET THE SMALLEST OBSERVED VALUE g.e. $value
    return ${ $self->{data}->get_val($key) }{'rt_tail_prob'};
  }
  else {
    cluck "WARNING: Only 1 argument expected.";
    return undef;
  }
}

sub get_data {
  my $self = shift;
  $self->_do_cdf() unless $self->{did_cdf};
  my ($uniqkeys, $sumweights, $keys, $weights, $counts, $cdfs, $rtps, $percentiles, $order) = ([],[],[],[],[],[],[],[],[]);
  my $key = $self->{'data'}->minimum();
  while (defined $key){
    my $value = $self->{data}->get_val($key);
    push @$uniqkeys, $key;
    push @$sumweights, $$value{'weight'};
    foreach my $weight (@{ $$value{'weights'} } ) {
      push @$keys, $key;
      push @$weights, $weight;
    }
    push @$order, @{ $$value{'order'}};
    push @$counts, $$value{'count'};
    push @$cdfs,   $$value{'cdf'};
    push @$rtps,   $$value{'rt_tail_prob'};
    push @$percentiles,   $$value{'percentile'};
    $key = $self->{'data'}->successor($key);
  }
  return {'uniqvars' => $uniqkeys, 'sumweights' => $sumweights, 'counts' => $counts, 'cdfs' => $cdfs, 'rtps' => $rtps, 'vars' => $keys, 'weights' => $weights, 'percentiles' => $percentiles, 'order' => $order};
}

sub print {
  my $self = shift;
  print Data::Dumper->Dump([$self->get_data()]);
}

## OVERRIDES FOR UNSUPPORTED FUNCTIONS

sub sort_data{
  confess "ERROR: Statistics::Descriptive::Weighted does not support this function.";
}

sub presorted{
  confess "ERROR: Statistics::Descriptive::Weighted does not support this function.";
}

sub harmonic_mean{
  confess "ERROR: Statistics::Descriptive::Weighted does not support this function.";
}

sub geometric_mean{
  confess "ERROR: Statistics::Descriptive::Weighted does not support this function.";
}

sub trimmed_mean{
  confess "ERROR: Statistics::Descriptive::Weighted does not support this function.";
}

sub frequency_distribution{
  confess "ERROR: Statistics::Descriptive::Weighted does not support this function.";
}

sub least_squares_fit{
  confess "ERROR: Statistics::Descriptive::Weighted does not support this function.";
}


1;

package Statistics::Descriptive;

##All modules return true.
1;

__END__