File Coverage

blib/lib/Statistics/Descriptive/Discrete.pm

Criterion	Covered	Total	%
statement	220	226	97.3
branch	62	78	79.4
condition	11	15	73.3
subroutine	17	18	94.4
pod	8	9	88.8
total	318	346	91.9

package Statistics::Descriptive::Discrete;

### This module draws heavily from Statistics::Descriptive

use strict;
use warnings;
use Carp;
use AutoLoader;
use vars qw($VERSION $AUTOLOAD $DEBUG $Tolerance %autosubs);

$VERSION = '0.11';
$DEBUG = 0;

#see Statistics::Descriptive documentation for use of $Tolerance
$Tolerance = 0.0;

#what subs can be autoloaded?
%autosubs = (
  count					=> undef,
  mean					=> undef,
  geometric_mean=> undef,
  harmonic_mean=>undef,
  sum					=> undef,
  mode					=> undef,
  median				=> undef,
  min					=> undef,
  max					=> undef,
  mindex			=> undef,
  maxdex			=> undef,
  standard_deviation	=> undef,
  sample_range			=> undef,
  variance				=> undef,
  text					=> undef,
);


sub new
{
	my $proto = shift;
	my $class = ref($proto) || $proto;
	my $self = {};
	$self->{_permitted} = \%autosubs;
	$self->{data} = ();
	$self->{_dataindex} = (); #index of where each value first seen when adding data
	$self->{dirty} = 1; #is the data dirty?
	$self->{_index} = 0; #current index of number of data items added

	bless ($self,$class);
	print __PACKAGE__,"->new(",join(',',@_),")\n" if $DEBUG;
	return $self;
}

# Clear the stat object & erase all data
# Object will be ready to use as if new was called
# Not sure this is more efficient than just creating a new object but
# maintained for compatability with Statistics::Descriptive
sub clear
{
	my $self = shift;
    my %keys = %{ $self };

	#remove _permitted from the deletion list
    delete $keys{"_permitted"};

    foreach my $key (keys %keys)
	{ # Check each key in the object
		print __PACKAGE__,"->clear, deleting $key\n" if $DEBUG;
        delete $self->{$key};  # Delete any out of date cached key
    }
	$self->{data} = ();
	$self->{_dataindex} = ();
	$self->{dirty} = 1;
	$self->{_index} = 0;
}

sub add_data
{
	#add data but don't compute ANY statistics yet
	my $self = shift;
	print __PACKAGE__,"->add_data(",join(',',@_),")\n" if $DEBUG;

	#get each element and add 0 to force it be a number
	#that way, 0.000 and 0 are treated the same
	my $val = shift;
	while (defined $val)
	{
		$val += 0;
		$self->{data}{$val}++;
		if (not exists $self->{_dataindex}{$val}) {
			$self->{_dataindex}{$val} = $self->{_index};
		}
		$self->{_index}++;
		#set dirty flag so we know cached stats are invalid
		$self->{dirty}++;
		$val = shift; #get next element
	}
}

sub add_data_tuple
{
	#add data but don't compute ANY statistics yet
	#the data are pairs of values and occurrences
	#e.g. 4,2 means 2 occurrences of the value 4
	#thanks to Bill Dueber for suggesting this

	my $self = shift;
	print __PACKAGE__,"->add_data_tuple(",join(',',@_),")\n" if $DEBUG;

	#we want an even number of arguments (tuples in the form (value, count))
	carp "argument list must have even number of elements" if @_ % 2;

	#get each element and add 0 to force it be a number
	#that way, 0.000 and 0 are treated the same
	#if $count is 0, then this will set the dirty flag but have no effect on
	#the statistics
	my $val = shift;
	my $count = shift;
	while (defined $count)
	{
		$val += 0;
		$self->{data}{$val} += $count;
		if (not exists $self->{_dataindex}{$val}) {
			$self->{_dataindex}{$val} = $self->{_index};
		}
		$self->{_index} += $count;
		#set dirty flag so we know cached stats are invalid
		$self->{dirty}++;
		$val = shift; #get next element
		$count = shift;
	}
}

sub _test_for_too_small_val
{
    my $self = shift;
    my $val = shift;

    return (abs($val) <= $Statistics::Descriptive::Discrete::Tolerance);
}

sub _calc_harmonic_mean
{
    my $self = shift;
		my $count = shift;
		my $datakeys = shift; #array ref

    my $hs = 0;

    foreach my $val ( @{$datakeys} )
    {
        ##Guarantee that there are no divide by zeros
        if ($self->_test_for_too_small_val($val))
        {
            return;
        }

				foreach (1..$self->{data}{$val})
				{
        	$hs += 1/$val;
				}
    }

    if ($self->_test_for_too_small_val($hs))
    {
        return;
    }

    return $count/$hs;
}

sub _all_stats
{
	#compute all the stats in one sub to save overhead of sub calls
	#a little wasteful to do this if all we want is count or sum for example but
	#I want to keep add_data as lean as possible since it gets called a lot
	my $self = shift;
	print __PACKAGE__,"->_all_stats(",join(',',@_),")\n" if $DEBUG;

	#if data is empty, set all stats to undef and return
	if (!$self->{data})
	{
		foreach my $key (keys %{$self->{_permitted}})
		{
			$self->{$key} = undef;
		}
		$self->{count} = 0;
		return;
	}

	#count = total number of data values we have
	my $count = 0;
	$count += $_ foreach (values %{$self->{data}});

	my @datakeys = keys %{$self->{data}};

	#initialize min, max, mode to an arbitrary value that's in the hash
	my $default = $datakeys[0];
	my $max  = $default;
	my $min  = $default;
	my $mode = $default;
	my $moden = 0;
	my $sum = 0;

	#find min, max, sum, and mode
	foreach (@datakeys)
	{
		my $n = $self->{data}{$_};
		$sum += $_ * $n;
		$min = $_ if $_ < $min;
		$max = $_ if $_ > $max;

		#only finds one mode but there could be more than one
		#also, there might not be any mode (all the same frequency)
		#todo: need to make this more robust
		if ($n > $moden)
		{
			$mode = $_;
			$moden = $n;
		}
	}
	my $mindex = $self->{_dataindex}{$min};
	my $maxdex = $self->{_dataindex}{$max};

	my $mean = $sum/$count;

	my $stddev = 0;
	my $variance = 0;

	if ($count > 1)
	{
		# Thanks to Peter Dienes for finding and fixing a round-off error
		# in the following variance calculation

		foreach my $val (@datakeys)
		{
			$stddev += $self->{data}{$val} * (($val - $mean) ** 2);
		}
		$variance = $stddev / ($count - 1);
		$stddev = sqrt($variance);
	}
	else {$stddev = undef}

	#find median, and do it without creating a list of the all the data points
	#if n=count is odd and n=2k+1 then median = data(k+1)
	#if n=count is even and n=2k, then median = (data(k) + data(k+1))/2
	my $odd = $count % 2; #odd or even number of points?
	my $even = !$odd;
	my $k = $odd ? ($count-1)/2 : $count/2;
	my $median = undef;
	my $temp = 0;
	MEDIAN: foreach my $val (sort {$a <=> $b} (@datakeys))
	{
		foreach (1..$self->{data}{$val})
		{
			$temp++;
			if (($temp == $k) && $even)
			{
				$median += $val;
			}
			elsif ($temp == $k+1)
			{
				$median += $val;
				$median /= 2 if $even;
				last MEDIAN;
			}
		}
	}

	#compute geometric mean
	my $gm = 1;
	my $exponent = 1/$count;
	foreach my $val (@datakeys)
		{
				if ($val < 0)
				{
						$gm = undef;
						last;
				}
				foreach (1..$self->{data}{$val})
				{
					$gm *= $val**$exponent;
				}
		}

	#compute harmonic mean
	my $harmonic_mean = scalar $self->_calc_harmonic_mean($count, \@datakeys);

	print __PACKAGE__,"count: $count, _index ",$self->{_index},"\n" if $DEBUG;

	$self->{count}  = $count;
	$self->{sum}    = $sum;
	$self->{standard_deviation} = $stddev;
	$self->{variance} = $variance;
	$self->{min}    = $min;
	$self->{max}    = $max;
	$self->{mindex} = $mindex;
	$self->{maxdex} = $maxdex;
	$self->{sample_range} = $max - $min; #todo: does this require any bounds checking
	$self->{mean}    = $mean;
	$self->{geometric_mean} = $gm;
	$self->{harmonic_mean} = $harmonic_mean;
	$self->{median} = $median;
	$self->{mode}   = $mode;

	#clear dirty flag so we don't needlessly recompute the statistics
	$self->{dirty} = 0;
}

sub set_text
{
	my $self = shift;
	$self->{text} = shift;
}

sub get_data
{
	#returns a list of the data in sorted order
	#the list could be very big an this defeat the purpose of using this module
	#use this only if you really need it
	my $self = shift;
	print __PACKAGE__,"->get_data(",join(',',@_),")\n" if $DEBUG;

	my @data;
	foreach my $val (sort {$a <=> $b} (keys %{$self->{data}}))
	{
		push @data, $val foreach (1..$self->{data}{$val});
	}
	return @data;
}

# this is the previous frequency_distribution code
# redid this completely based on current implementation in
# Statistics::Descriptive
# sub frequency_distribution
# {
# 	#Compute frequency distribution (histogram), borrowed heavily from Statistics::Descriptive
# 	#Behavior is slightly different than Statistics::Descriptive
# 	#e.g. if partition is not specified, we use  to set the number of partitions
# 	#     if partition = 0, then we return the data hash WITHOUT binning it into equal bins
# 	#	  I often want to just see how many of each value I saw
# 	#Also, you can manually pass in the bin info (min bin, bin size, and number of partitions)
# 	#I don't cache the frequency data like Statistics::Descriptive does since it's not as expensive to compute
# 	#but I might add that later
# 	#todo: the minbin/binsize stuff is funky and not intuitive -- fix it
# 	my $self = shift;
# 	print __PACKAGE__,"->frequency_distribution(",join(',',@_),")\n" if $DEBUG;

# 	my $partitions = shift; #how many partitions (bins)?
# 	my $minbin = shift; #upper bound of first bin
# 	my $binsize = shift; #how wide is each bin?

# 	#if partition == 0, then return the data hash
# 	if (not defined $partitions || ($partitions == 0))
# 	{
# 		$self->{frequency_partitions} = 0;
# 		%{$self->{frequency}} = %{$self->{data}};
# 		return %{$self->{frequency}};
# 	}

# 	#otherwise, partition better be >= 1
# 	return undef unless $partitions >= 1;

# 	$self->_all_stats() if $self->{dirty}; #recompute stats if dirty, (so we have count)
# 	return undef if $self->{count} < 2; #must have at least 2 values

# 	#set up the bins
# 	my ($interval, $iter, $max);
# 	if (defined $minbin && defined $binsize)
# 	{
# 		$iter = $minbin;
# 		$max = $minbin+$partitions*$binsize - $binsize;
# 		$interval = $binsize;
# 		$iter -= $interval; #so that loop that sets up bins works correctly
# 	}
# 	else
# 	{
# 		$iter = $self->{min};
# 		$max = $self->{max};
# 		$interval = $self->{sample_range}/$partitions;
# 	}
# 	my @k;
# 	my %bins;
# 	while (($iter += $interval) < $max)
# 	{
# 		$bins{$iter} = 0;
# 		push @k, $iter;
# 	}
# 	$bins{$max} = 0;
# 	push @k, $max;

# 	VALUE: foreach my $val (keys %{$self->{data}})
# 	{
# 		foreach my $k (@k)
# 		{
# 			if ($val <= $k)
# 			{
# 				$bins{$k} += $self->{data}{$val};  #how many of this value do we have?
# 				next VALUE;
# 			}
# 		}
# # 	}

# 	%{$self->{frequency}} = %bins;   #save it for later in case I add caching
# 	$self->{frequency_partitions} = $partitions; #in case I add caching in the future
# 	return %{$self->{frequency}};
# }

sub frequency_distribution_ref
{
    my $self = shift;
    my @k = ();

		# If called with no parameters, return the cached hashref
		# if we have one and data is not dirty
		# This is implemented this way because that's how Statistics::Descriptive
		# implements this.  I don't like it.
  	if ((!@_) && (! $self->{dirty}) && (defined $self->{_frequency}))
    {
        return $self->{_frequency};
    }

		$self->_all_stats() if $self->{dirty}; #recompute stats if dirty, (so we have count)

    # Must have at least two elements
    if ($self->count() < 2)
    {
        return undef;
    }

    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 (defined $partitions) && ($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 (keys %{$self->{data}})
    {
        foreach my $limit (@k)
        {
            if ($element <= $limit)
            {
                $bins{$limit} += $self->{data}{$element};
                next ELEMENT;
            }
        }
    }

		$self->{_frequency} = \%bins;
    return $self->{_frequency};
}

sub frequency_distribution {
    my $self = shift;

    my $ret = $self->frequency_distribution_ref(@_);

    if (!defined($ret))
    {
        return undef;
    }
    else
    {
        return %$ret;
    }
}

# return count of unique values in data if called in scalar context
# returns sorted array of unique data values if called in array context
# returns undef if no data
sub uniq
{
	my $self = shift;

	if (!$self->{data})
	{
		return undef;
	}

	my @datakeys = sort {$a <=> $b} keys %{$self->{data}};

	if (wantarray)
	{
		return @datakeys;
	}
	else
	{
		my $uniq = @datakeys;
		return $uniq;
	}
}

sub AUTOLOAD {
	my $self = shift;
	my $type = ref($self)
		or croak "$self is not an object";
	my $name = $AUTOLOAD;
	$name =~ s/.*://;     ##Strip fully qualified-package portion
	return if $name eq "DESTROY";
	unless (exists $self->{_permitted}{$name} ) {
		croak "Can't access `$name' field in class $type";
	}

	print __PACKAGE__,"->AUTOLOAD $name\n" if $DEBUG;

	#compute stats if necessary
	$self->_all_stats() if $self->{dirty};
	return $self->{$name};
}

1;

__END__