C:/Documents and Settings/Stephen/Desktop/Biostatistics Library and Program/Biostatistics/FisherExactArray.cs

Go to the documentation of this file.

using System;

namespace StephenAshley.Biostatistics
{
   using NUMBER = Decimal;
   public class FisherExactArray : ChiSquaredArray
   {
      public FisherExactArray() : base(2, 2) { }

      public override NUMBER p()
      {
         return FisherProbability((UInt16)outcomesArray[0, 0].iObserved,
               (UInt16)outcomesArray[0, 1].iObserved, (UInt16)outcomesArray[1, 0].iObserved,
               (UInt16)outcomesArray[1, 1].iObserved);
      }

      private NUMBER FisherProbability(UInt16 uA1, UInt16 uB1, UInt16 uA2, UInt16 uB2)
      {
         Int32[,] f = new Int32[2, 2];
         Int32[,] a = new Int32[2, 2];
         Int32[] iRowSums = new Int32[2];
         Int32[] iColumnSums = new Int32[2];

         int ntotal, 
            i, j, 
            ismall, jsmall;
         NUMBER ptable, poriginal_table, ptotal;
         NUMBER numerator;
         UInt16[,] data = new UInt16[2, 2];

         poriginal_table = 0.0m;
         data[0, 0] = uA1;
         data[0, 1] = uB1;
         data[1, 0] = uA2;
         data[1, 1] = uB2;


         /* copy data into integer table */
         ntotal = 0;
         for (i = 0; i < 2; i++)
            for (j = 0; j < 2; j++)
            {
               f[i, j] = data[i, j];
               ntotal = ntotal + f[i, j];
            }

         /* find value in table to be put in 1,1 position */
         if (f[0, 0] * f[1, 1] <= f[0, 1] * f[1, 0])
         {
            if (f[0, 0] < f[1, 1])
            {
               ismall = 1;
               jsmall = 1;
            }
            else
            {
               ismall = 2;
               jsmall = 2;
            }
         }
         else
         {
            if (f[0, 1] < f[1, 0])
            {
               ismall = 1;
               jsmall = 2;
            }
            else
            {
               ismall = 2;
               jsmall = 1;
            }
         }

         /* reflect the table to get the smallest element in the 1,1 position */
         if ((ismall == 1) && (jsmall == 1))
         {
            a[0, 0] = f[0, 0];
            a[0, 1] = f[0, 1];
            a[1, 0] = f[1, 0];
            a[1, 1] = f[1, 1];
         }
         else if ((ismall == 1) && (jsmall == 2))
         {
            a[0, 0] = f[0, 1];
            a[0, 1] = f[0, 0];
            a[1, 0] = f[1, 1];
            a[1, 1] = f[1, 0];
         }
         else if ((ismall == 2) && (jsmall == 1))
         {
            a[0, 0] = f[1, 0];
            a[0, 1] = f[1, 1];
            a[1, 0] = f[0, 0];
            a[1, 1] = f[0, 1];
         }
         else if ((ismall == 2) && (jsmall == 2))
         {
            a[0, 0] = f[1, 1];
            a[0, 1] = f[1, 0];
            a[1, 0] = f[0, 1];
            a[1, 1] = f[0, 0];
         }

         for (i = 0; i < 2; i++)
         {
            iRowSums[i] = a[i, 0] + a[i, 1];
            iColumnSums[i] = a[0, i] + a[1, i];
         }

         /* find probability of first tail */
         numerator = DecMath.Exp(LogFactorial(iRowSums[0]) + LogFactorial(iRowSums[1]) +
               LogFactorial(iColumnSums[0]) + LogFactorial(iColumnSums[1]) - 
               LogFactorial(ntotal));

         ptotal = 0.0m;
         ptable = 0.0m;
         for (i = 0; i <= a[0, 0]; i++)
         {
            ptable = numerator / DecMath.Exp(LogFactorial(a[0, 0] - i) + 
               LogFactorial(a[0, 1] + i) + LogFactorial(a[1, 0] + i)
                  + LogFactorial(a[1, 1] - i));

            if (i == 0)
               poriginal_table = ptable;

            ptotal = ptotal + ptable;
         }

         /* find probability of second tail */
         if (a[0, 1] > a[1, 0])
            i = a[1, 0];
         else
            i = a[0, 1];

         while ((i != 0) && (ptable <= poriginal_table))
         {
            ptable = numerator / DecMath.Exp(LogFactorial(a[0, 0] + i) + 
               LogFactorial(a[0, 1] - i) + LogFactorial(a[1, 0] - i) + 
               LogFactorial(a[1, 1] + i));

            if (ptable <= poriginal_table)
               ptotal = ptotal + ptable;

            i = i - 1;
         }

         return ptotal;
      }
      /* FISHER */

      /* Computes the natural log of f! */
      private NUMBER LogFactorial(Int32 f)
      {
         NUMBER rtempd;

         rtempd = 0.0m;
         for (int i = 1; i <= f; i++)
            rtempd += DecMath.Log((NUMBER)i);

         return rtempd;
      }
   }
}

---