#!/usr/bin/perl -w

use Statistics::Table::F;

$designs = [[18, 22, 17, 10, 34, 15, 12, 20, 21],
            [21, 34, 18, 18, 20, 22, 17, 19, 14, 10, 21],
            [21, 25, 28, 27, 30, 18, 26, 25, 25, 29],
            [25, 17, 19, 22, 18, 18, 22, 30]];

if (($F = anova($designs)) >=
        F(@$designs-1, count_elements($designs) - @$designs, 0.05)) {
    print "F is $F; the difference between designs is significant.\n";
} else {
    print "F is $F; the data are not sufficient for significance.\n";
}

sub count_elements {
    my ($arraysref) = shift;
    my $result;
    foreach (@$arraysref) { $result += @$_ }
    return $result;
}

__END__
# This code is already in Statistics::Table::F, but we include it here 
# for convenience.
#

sub mean {
    my ($arrayref) = @_;
    my $result;
    foreach (@$arrayref) { $result += $_ }
    return $result / @$arrayref;
}

sub estimate_variance {
    my ($arrayref) = @_;
    my ($mean) = mean($arrayref);
    my ($result);
    foreach (@$arrayref) {
        $result += ($_ - $mean) ** 2;
    }
    return $result / $#{$arrayref};
}

sub square_sum {
    my ($arraysref) = shift;
    my (@arrays) = @$arraysref;
    my ($result, $arrayref);
    foreach $arrayref (@arrays) {
        foreach (@$arrayref) { $result += $_ ** 2 }
    }
    return $result;
}

sub sum {
    my ($arraysref) = shift;
    my (@arrays) = @$arraysref;
    my ($result, $arrayref);
    foreach $arrayref (@arrays) {
        foreach (@$arrayref) { $result += $_ }
    }
    return $result;
}

sub square_groups {
    my ($arraysref) = shift;
    my (@arrays) = @$arraysref;
    my ($result, $arrayref);
    foreach $arrayref (@arrays) {
        my $sum = 0;
        foreach (@$arrayref) { $sum += $_ }
        $result += ($sum ** 2) / @$arrayref;
    }
    return $result;
}

# Performs a one-way analysis of variance, returning the F-ratio.
sub anova {
    my ($all) = shift;
    my $num_of_elements = count_elements($all);
    my $square_of_everything = square_sum($all);
    my $sum_of_everything = sum($all);
    my $sum_of_groups = square_groups($all);
    my $degrees_of_freedom_within  = $num_of_elements - @$all;
    my $degrees_of_freedom_between = @$all - 1;
    my $sum_of_squares_within = $square_of_everything - $sum_of_groups;
    my $mean_of_squares_within = $sum_of_squares_within /
        $degrees_of_freedom_within;
    my $sum_of_squares_between = $sum_of_groups -
        ($sum_of_everything ** 2)/$num_of_elements;
    my $mean_of_squares_between = $sum_of_squares_between /
        $degrees_of_freedom_between;
    return $mean_of_squares_between / $mean_of_squares_within;
}