Beefy Boxes and Bandwidth Generously Provided by pair Networks Bob
No such thing as a small change
 
PerlMonks  

Fisher's Exact test (probability of 2x2 tables)

by jeroenes (Priest)
on Oct 03, 2001 at 12:19 UTC ( #116378=sourcecode: print w/ replies, xml ) Need Help??

Category: Misc
Author/Contact Info Jeroen Elassaiss-Schaap, just /msg me at jeroenes

Authors of the original code: Kristopher J. Preacher and Nancy E. B

Description: This is a plain translation of this page into perl, kept close to the original code.
package Statistics::Fishersexact;

use strict;
use warnings;

use constant DECIMALS=>8;

sub round{
  my $frac= shift; my $dec = shift;
  sprintf("%".$dec."f",$frac);
}

sub fac {
  my $z = shift;
  if ($z==0) { return 1; }
  if ($z>1) { return $z*fac($z-1); }
  else { return $z; }
}

sub sum {
  my $sum = 0;
  $sum += $_ for @_;
  $sum;
}

sub p_fisher {
  my $f11=shift;
  my $f12=shift;
  my $f21=shift;
  my $f22=shift;
  my $row1 = sum($f11,$f12);
  my $row2 = sum($f21,$f22);
  my $col1 = sum($f11,$f21);
  my $col2 = sum($f12,$f22);
  my $r1=sum($f11,$f12);
  my $r2=sum($f21,$f22);
  my $c1=sum($f11,$f21);
  my $c2=sum($f12,$f22);
  my $total = sum($r1,$r2);
  my ($a,$b,$c,$d,$z,$temp,$p,$ratio,$p1,$aaron,$flag);
  if ($r1<=$r2) { if ($c1<=$c2) { $a=$f11;$b=$f12;$c=$f21;$d=$f22; } }
  if ($r1<=$r2) { if ($c2<=$c1) { $a=$f12,$b=$f22;$c=$f11;$d=$f21; } }
  if ($r2<=$r1) { if ($c1<=$c2) { $a=$f21;$b=$f11;$c=$f22;$d=$f12; } }
  if ($r2<=$r1) { if ($c2<=$c1) { $a=$f22;$b=$f21;$c=$f12;$d=$f11; } }
  if ($b<=$a) { $z=$c;$c=$a;$a=$b;$b=$d;$d=$z; }
  elsif ($c<=$a) { $z=$b;$b=$a;$a=$c;$c=$d;$d=$z; }
  $p1=0;$temp=$a;$ratio=0;$z=0;
  my $i;
  my ($ab,$ac,$bd,$cd,$abcd,$x1,$x2,$x3,$x4,$y1,$y2,$y3,$y4,$ra1,$ra2,
+$ra3,$ra4,$ra5);
  for ($i=0;$i<10000;$i++) {
    $ab=sum($a,$b);$ac=sum($a,$c);$bd=sum($b,$d);
    $cd=sum($c,$d);$abcd=sum($ab,$cd);
    $x1=fac($ab);$x2=fac($ac);$x3=fac($bd);$x4=fac($cd);
    $y1=fac($a); $y2=fac($b); $y3=fac($c); $y4=fac($d);
    for (my $i=0;$i<2;$i++) {
      if ($x2>$x1) { $z=$x1;$x1=$x2;$x2=$z }
      if ($x3>$x2) { $z=$x2;$x2=$x3;$x3=$z }
      if ($x4>$x3) { $z=$x3;$x3=$x4;$x4=$z }
        }
    for (my $i=0;$i<2;$i++) {
      if ($y2>$y1) { $z=$y1;$y1=$y2;$y2=$z }
      if ($y3>$y2) { $z=$y2;$y2=$y3;$y3=$z }
      if ($y4>$y3) { $z=$y3;$y3=$y4;$y4=$z }
    }
    $ra1=$x1/$y1;$ra2=$x2/$y2;$ra3=$x3/$y3; $ra4=$x4/$y4;
    $ra5=1/fac($abcd);
    $ratio=$ra1*$ra2*$ra3*$ra4*$ra5;
    if ($p1==0) {$aaron=$ratio }
    $p1=sum($p1,$ratio);
    $a=sum($a,-1); $b=sum($b,1); $c=sum($c,1);$d=sum($d,-1);
    if ($a<0) { last }
  }
  my $cp1=round($p1,DECIMALS);
  if ($p1<.00000001) { if ($p1>0) { $cp1=".00000001" } }
  if ($p1>.99999999) { if ($p1!=1) { $cp1=".99999999" } }
  if ($r1<=$r2) { if ($c1<=$c2) { $a=$f11;$b=$f12;$c=$f21;$d=$f22 } }
  if ($r1<=$r2) { if ($c2<=$c1) { $a=$f12,$b=$f22;$c=$f11;$d=$f21 } }
  if ($r2<=$r1) { if ($c1<=$c2) { $a=$f21;$b=$f11;$c=$f22;$d=$f12 } }
  if ($r2<=$r1) { if ($c2<=$c1) { $a=$f22;$b=$f21;$c=$f12;$d=$f11 } }
  if ($b<=$a) { $z=$c;$c=$a;$a=$b;$b=$d;$d=$z }
  elsif ($c<=$a) { $z=$b;$b=$a;$a=$c;$c=$d;$d=$z }
  my $p2=0;$ratio=0;$flag=0;
  for (my $i=0;$i<10000;$i++) {
    $a=sum($a,1);$b=sum($b,-1);$c=sum($c,-1);$d=sum($d,1);
    if ($b<0) { last }
    if ($c<0) { last }
    $ab=sum($a,$b);$ac=sum($a,$c);$bd=sum($b,$d);$cd=sum($c,$d);$abcd=
+sum($ab,$cd);
    $x1=fac($ab);$x2=fac($ac);$x3=fac($bd);$x4=fac($cd);
    $y1=fac($a);$y2=fac($b);$y3=fac($c);$y4=fac($d);
    for ($i=0;$i<2;$i++) {
      if ($x2>$x1) { $z=$x1;$x1=$x2;$x2=$z }
      if ($x2>$x1) { $z=$x1;$x1=$x2;$x2=$z }
      if ($x3>$x2) { $z=$x2;$x2=$x3;$x3=$z }
      if ($x4>$x3) { $z=$x3;$x3=$x4;$x4=$z }
    }
    for ($i=0;$i<2;$i++) {
      if ($y2>$y1) { $z=$y1;$y1=$y2;$y2=$z }
      if ($y3>$y2) { $z=$y2;$y2=$y3;$y3=$z }
      if ($y4>$y3) { $z=$y3;$y3=$y4;$y4=$z }
    }
    $ra1=$x1/$y1;$ra2=$x2/$y2;$ra3=$x3/$y3;$ra4=$x4/$y4;$ra5=1/fac($ab
+cd);
    $ratio=$ra1*$ra2*$ra3*$ra4*$ra5;
    if ($ratio<=$aaron) { $flag=1 }
    if ($flag==1) { $p2=sum($p2,$ratio) }
  }
  my $cp2=round(sum($p1,$p2),DECIMALS);
  if ($p2<.00000001) { if ($p2>0) { $cp2=".00000001" } }
  if ($p2>.99999999) { if ($p2!=1) { $cp2=".99999999" } }
  return ($cp1,$cp2);
}
 
 
=head1 DESCRIPTION
 
This is perlization of
http://quantrm2.psy.ohio-state.edu/kris/fisher/fisher.htm. It calculat
+es
fisher's Exact test, one and two tailed, from a 2x2 table.
 
($p1, $p2) = p_fisher( $cell11, $cell12, $cell21, $cell22);
 
=cut
 
 
1;

Comment on Fisher's Exact test (probability of 2x2 tables)
Download Code
Re: Fisher's Exact test (probability of 2x2 tables)
by Zaxo (Archbishop) on Oct 03, 2001 at 12:57 UTC

    I'd like to recommend a non-recursive &fishersexact::fac, for both speed and memory usage. Something like

    fac { my ($arg, $res) = ($_[0], 1); die "Negative argument to factorial\n" if $arg < 0; $res *= $_ for 1..$arg; return $res; }

    Additionally factorials soon overflow perl int(). An implementation in terms of Math::BigInt, Math::GMP or Math::Pari might be in order. Overflow might be abated by rewriting the formulas in terms of binomial coefficients.

    After Compline,
    Zaxo

      If I look at the code, there is a very lot more to improve. It seems to neglect an indecent number of software design rules. Right now, I'm just happy to have kept to my promise in Re: Fisher exact - one tailed

      Jeroen
      "We are not alone"(FZ)

Re: Fisher's Exact test (probability of 2x2 tables)
by merlyn (Sage) on Oct 03, 2001 at 18:30 UTC
    The package name should follow the rules for module package names. You need an uppercase character (lowercase are reserved for pragmas), and you should probably get a name for this module from the module naming czars before distributing it (even posting on perlmonks is a form of "distribution").

    -- Randal L. Schwartz, Perl hacker

      I obviously didn't go that official. Even worse, I have never named any module that I wrote in the official way. Will do from now one, of course.

      Thanx for the tip.

      Jeroen
      "We are not alone"(FZ)

Back to Code Catacombs

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: sourcecode [id://116378]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others making s'mores by the fire in the courtyard of the Monastery: (7)
As of 2014-04-19 03:05 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    April first is:







    Results (475 votes), past polls