umcg.genetica.math.stats.FisherExactTest Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of genetica-libraries Show documentation
There is a newer version: 1.0.7
/*
 * FisherExactTest.java
 *
 * Created on June 6, 2007, 10:08 AM
 *
 * To change this template, choose Tools | Template Manager
 * and open the template in the editor.
 */

package umcg.genetica.math.stats;

/**
 *
 * @author ludefranke
 */
public class FisherExactTest {
    private double left;
    private double right;
    private double twotail;
    private double sleft;
    private double sright;
    private double sless;
    private double slarg;
    private int sn11;
    private int sn1_;
    private int sn_1;
    private int sn;
    private double sprob;
    private int n11_;
    private int n12_;
    private int n21_;
    private int n22_;
    
    /** Creates a new instance of FisherExactTest */
    public FisherExactTest() {

        /* 
        int n11 = 10;
        int n12 = 1;
        int n21 = 1;
        int n22 = 10;
        System.out.println(calculateFisherFromMatrix(n11, n12, n21, n22));
        if (1==1) System.exit(0);
        */
    }
    
    public double getFisherPValue(int n11, int n12, int n21, int n22) {
        n11_ = n11;
        n12_ = n12;
        n21_ = n21;
        n22_ = n22;
        if(n11_ < 0)
            n11_ *= -1;
        if(n12_ < 0)
            n12_ *= -1;
        if(n21_ < 0)
            n21_ *= -1;
        if(n22_ < 0)
            n22_ *= -1;
        int n1_ = n11_ + n12_;
        int n_1 = n11_ + n21_;
        int n = n11_ + n12_ + n21_ + n22_;
        exact(n11_, n1_, n_1, n);
        left = sless;
        right = slarg;
        twotail = sleft + sright;
        if(twotail > 1.0D)
            twotail = 1.0D;
        return twotail;
    }
    
    public double getFisherLeftTail() {
        return left;
    }
    
    public double getFisherRightTail() {
        return right;
    }
    
    public double calculateFisherTwoTail() {
        return twotail;
    }
    
    private static double lngamm(int z) {
        double x = 0.0D;
        x += 1.6594701874084621E-07D / (double)(z + 7);
        x += 9.9349371139307475E-06D / (double)(z + 6);
        x -= 0.1385710331296526D / (double)(z + 5);
        x += 12.50734324009056D / (double)(z + 4);
        x -= 176.61502914983859D / (double)(z + 3);
        x += 771.32342877576741D / (double)(z + 2);
        x -= 1259.1392167222889D / (double)(z + 1);
        x += 676.52036812188351D / (double)z;
        x += 0.99999999999951827D;
        return (Math.log(x) - 5.5810614667953278D - (double)z) + ((double)z - 0.5D) * Math.log((double)z + 6.5D);
    }
    
    private static double lnfact(int n) {
        if(n <= 1)
            return 0.0D;
        else
            return lngamm(n + 1);
    }
    
    private static double lnbico(int n, int k) {
        return lnfact(n) - lnfact(k) - lnfact(n - k);
    }
    
    private static double hyper_323(int n11, int n1_, int n_1, int n) {
        return Math.exp((lnbico(n1_, n11) + lnbico(n - n1_, n_1 - n11)) - lnbico(n, n_1));
    }
    
    private double hyper(int n11) {
        return hyper0(n11, 0, 0, 0);
    }
    
    private double hyper0(int n11i, int n1_i, int n_1i, int ni) {
        if((n1_i == 0) & (n_1i == 0) & (ni == 0)) {
            if(n11i % 10 != 0) {
                if(n11i == sn11 + 1) {
                    sprob = sprob * (((double)sn1_ - (double)sn11) / (double)n11i) * (((double)sn_1 - (double)sn11) / (((double)n11i + (double)sn) - (double)sn1_ - (double)sn_1));
                    sn11 = n11i;
                    return sprob;
                }
                if(n11i == sn11 - 1) {
                    sprob = sprob * ((double)sn11 / ((double)sn1_ - (double)n11i)) * ((((double)sn11 + (double)sn) - (double)sn1_ - (double)sn_1) / ((double)sn_1 - (double)n11i));
                    sn11 = n11i;
                    return sprob;
                }
            }
            sn11 = n11i;
        } else {
            sn11 = n11i;
            sn1_ = n1_i;
            sn_1 = n_1i;
            sn = ni;
        }
        sprob = hyper_323(sn11, sn1_, sn_1, sn);
        return sprob;
    }
    
    private double exact(int n11, int n1_, int n_1, int n) {
        int max = n1_;
        if(n_1 < max)
            max = n_1;
        int min = (n1_ + n_1) - n;
        if(min < 0)
            min = 0;
        if(min == max) {
            sless = 1.0D;
            sright = 1.0D;
            sleft = 1.0D;
            slarg = 1.0D;
            return 1.0D;
        }
        double prob = hyper0(n11, n1_, n_1, n);
        sleft = 0.0D;
        double p = hyper(min);
        int i;
        for(i = min + 1; p < 0.99999998999999995D * prob; i++) {
            sleft += p;
            p = hyper(i);
        }
        
        i--;
        if(p < 1.0000000099999999D * prob)
            sleft += p;
        else
            i--;
        sright = 0.0D;
        p = hyper(max);
        int j;
        for(j = max - 1; p < 0.99999998999999995D * prob; j--) {
            sright += p;
            p = hyper(j);
        }
        
        j++;
        if(p < 1.0000000099999999D * prob)
            sright += p;
        else
            j++;
        if(Math.abs(i - n11) < Math.abs(j - n11)) {
            sless = sleft;
            slarg = (1.0D - sleft) + prob;
        } else {
            sless = (1.0D - sright) + prob;
            slarg = sright;
        }
        return prob;
    }
    
    
    
}