File Coverage

blib/lib/Statistics/Distributions/Bartlett.pm

Criterion	Covered	Total	%
statement	18	60	30.0
branch	0	4	0.0
condition			n/a
subroutine	6	8	75.0
pod	0	2	0.0
total	24	74	32.4

package Statistics::Distributions::Bartlett;

use warnings;
use strict;
use Carp;
use Statistics::Distributions qw/chisqrprob/;
use List::Util qw/sum/;

require Exporter;
our @ISA = qw(Exporter);
our @EXPORT = qw(bartlett);

=head1 NAME

Statistics::Distributions::Bartlett - Bartlett's test for equal sample variances.

=cut

=head1 VERSION

This document describes Statistics::Distributions::Bartlett version 0.0.2

=cut

=head1 SYNOPSIS

    use Statistics::Distributions::Bartlett;

    # Create data to pass as ARRAY/ARRAY references
    my @a = [qw/4534543  543 453 5 543 534 5 4 543 543 6 534 54 3 534  534 54 5/];
    my @b = [qw/ 546 565 65 64  5 434 65 457 67 78 6  34 2 6 57 43  2 556 7 4563/];
    my @c = [qw/ 565 44 535 6678 787 5 5/];

    # Call exported sub routine on ARRAY references of data to print result to STDOUT.
    &bartlett(\@a,\@b,\@c);

    # Call in LIST-context to access results directly:
    my ($K, $p_val, $df, $k, $n_total ) = &bartlett(\@a,\@b,\@c);

=cut

=head1 DESCRIPTION

Bartlett test is used to test if k samples have equal variances. Such homogeneity is often assumed by other statistical
tests and consequently the Bartlett test should be used to verify that assumption. See http://www.itl.nist.gov/div898/handbook/eda/section3/eda357.htm.

=cut

use version; our $VERSION = qv('0.0.2');

sub bartlett {
    my @groups = @_;
    my $k = scalar @groups;
    croak qq{\nThere must be more than one group} if ($k < 2);
    my $vars = &var(\@groups);
    my $n_total = sum map { scalar @{$_} } @groups;
    #my $SS_p = sum map { print qq{\nss $_->[2] and n $_->[0] and }, ($_->[2]-1) * $_->[0] ;($_->[2]-1) * $_->[0]  } @{$SSs};
    my $var_p = sum  map { ($_->[1]-1) * $_->[0]  } @{$vars};
    $var_p /= ($n_total-$k);
    $var_p = log($var_p);
    #my $SS_sum = sum map { print qq{\nss $_->[2] and n $_->[0] and }, log($_->[0]);($_->[2]-1) * log($_->[0])  } @{$SSs};
    my $var_sum = sum map { ($_->[1]-1) * log($_->[0])  } @{$vars};
    my $n_under = sum map { 1 / ($_->[1] - 1) } @{$vars};
    my $bar_k = ( ( ($n_total-$k) * $var_p ) - $var_sum ) / (1 + ( 1 / ( 3 * ($k-1))) * ( $n_under - ( 1 / ($n_total-$k)) ) ) ;
    my $df = $k -1;
    my $pval = &chisqrprob($df,$bar_k);
    if ( !wantarray ) { print qq{\nK = $bar_k\np_val = $pval\ndf = $df\nk = $k\ntotal n = $n_total}; return; }
    else { return ($bar_k, $pval, $df, $k, $n_total) }
    return;
}

sub var {
    my $groups = shift;
    my $result = [];
    for my $a_ref (@{$groups}) {
        # $stat->count()
        my $n = scalar(@{$a_ref});
        # $stat->sum();
        my $sum = sum @{$a_ref};
        # $stat->mean()
        my $mean = ( $sum / $n ) ;
        my $var = sum map { ($_-$mean)**2  } @{$a_ref};
        $var /= ($n-1);
        push @{$result}, [$var, $n];
    }
    return $result;
}

1; # Magic true value required at end of module

__END__