# rightFisher.pm version 0.1

# Statistical library package to calculate the Fishers exact test (right 
# sided). This should be used with statistic.pl and rank.pl
 
# A right sided test is calculated by adding the probabilities of all
# the possible two by two contingency tables formed by fixing the
# marginal totals and changing the value of n11 to greater than or
# equal to the given value. A right sided Fisher's Exact Test tells us
# how likely it is to randomly sample a table where n11 is greater
# than observed. In other words, it tells us how likely it is to sample
# an observation where the two words are more dependent than currently
# observed.

# Copyright (C)	2001-2003,
# Nitin Varma, University of Minnesota, Duluth
# varm0003@d.umn.edu
# Ted Pedersen,	University of Minnesota, Duluth
# tpederse@d.umn.edu
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation;	either version 2
# of the License, or (at your option) any later	version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along	with this program; if not, write to the	Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.


package	rightFisher;
require	Exporter;
@ISA = qw ( Exporter );
@EXPORT	= qw ( getStatisticName	calculateStatistic );

sub getStatisticName
{
    return "right sided Fisher's Exact test ";
}

sub calculateStatistic
{
    ($freq, $freq1, $freq2, $totalBigrams) = @_;

    my $n11 = $freq;
    my $npp = $totalBigrams;
    my $n1p = $freq1;
    my $np1 = $freq2;
    my $n2p = $npp - $n1p;
    my $np2 = $npp - $np1;

    # we shall have two	arrays one for the numerator and one for the
    # denominator. the arrays will contain the factorial upper limits. we
    # shall be arrange these two arrays	in descending order. while doing the
    # actual calculation, we shall take	a numerator/denominator	pair, and
    # go from the lower	value to the higher value, in effect doing a
    # "cancellation" of	sorts.

    # first create the numerator
    my @numerator = sort { $b <=> $a } ($n1p, $np1, $n2p, $np2);

    # now to the real calculation!!!
    my $probability = 0;
    my $i;
    my $j;

    # we shall calculate for n11 = 0. thereafter we shall just multiply	and
    # divide the result	for 0 with correct numbers to obtain result for	i,
    # i>0, i<=n11!! :o)
    $lower_limit = $n11;
    if($freq1 > $freq2)
    {
	$upper_limit = $freq2;
    }
    else
    {
	$upper_limit = $freq1;
    }

    my $n12 = $n1p-$n11;
    my $n21 = $np1-$n11;
    my $n22 = $n2p - $n21;
    while($n22 < 0 or $n12 < 0 or $n21 < 0)
    {
	$n11++;
        $n12 = $n1p - $n11;
        $n21 = $np1 - $n11;
        $n22 = $n2p - $n21;
    }

    my @denominator = sort { $b	<=> $a } ($npp,	$n22, $n12, $n21, $n11);

    # now that we have our two arrays all nicely sorted	and in place,
    # lets do the calculations!

    my @dLimits	= ();
    my @nLimits	= ();
    my $dIndex = 0;
    my $nIndex = 0;

    for	( $j = 0; $j < 4; $j ++	)
    {
	if ( $numerator[$j] > $denominator[$j] )
	{
	    $nLimits[$nIndex] =	$denominator[$j] + 1;
	    $nLimits[$nIndex+1]	= $numerator[$j];
	    $nIndex += 2;
	}

	elsif (	$denominator[$j] > $numerator[$j] )
	{
	    $dLimits[$dIndex] =	$numerator[$j] + 1;
	    $dLimits[$dIndex+1]	= $denominator[$j];
	    $dIndex += 2;
	}
    }
    $dLimits[$dIndex] =	1;
    $dLimits[$dIndex+1]	= $denominator[4];

    my $product	= 1;
    while ( defined ( $nLimits[0] ) )
    {
	while (	( $product < 10000 ) &&	( defined ( $nLimits[0]	) ) )
	{
	    $product *=	$nLimits[0];
	    $nLimits[0]++;
	    if ( $nLimits[0] > $nLimits[1] )
	    {
		shift @nLimits;
		shift @nLimits;
	    }
	}

	while (	$product > 1 )
	{
	    $product /=	$dLimits[0];
	    $dLimits[0]++;
	    if ( $dLimits[0] > $dLimits[1] )
	    {
		shift @dLimits;
		shift @dLimits;
	    }
	}
    }

    while ( defined ( $dLimits[0] ) )
    {
	$product /= $dLimits[0];
	$dLimits[0]++;
	if ( $dLimits[0] > $dLimits[1] )
	{
	    shift @dLimits;
	    shift @dLimits;
	}
    }

    # $product now has the hypergeometric probability for n11 =	0. add it to
    # the cumulative probability
    $probability += $product;

    # now for the rest of n11's	!!

    for	( $i = $n11+1; $i <=  $upper_limit; $i++ )
    {
	$product *= $n12;
	$n22++;
	$product /= $n22;
	$product *= $n21;
	$n12--;
	$n21--;
	$product /= $i;

	# thats	our new	probability for	n11 = i! :o)) cool eh? ;o))
	# add it to the	main probability! :o))

	$probability +=	$product; # !! :o)
    }

    # thats abt	it!
    return (1-$probability);
}

1;
