AssociationScoreTools.java

package edu.umd.hooka.alignment;

public class AssociationScoreTools {

	/**
	 * Computes the log-likelihood ratio 
	 * @param countEF
	 * @param countE
	 * @param countF
	 * @param N
	 */
	public static double computeLLR(int countEF, int countE, int countF, int N) {
		double cEF = countEF;
		double cE = countE;
		double cF = countF;
		double numS = N;
		double countnotE = numS - cE;
		double countnotF = numS - cF;
		double countEnotF = cE - cEF;
		double countnotEF = cF - cEF;
		double countnotEnotF = numS - (countnotEF + countEnotF + cEF);
		double pEF = 1.0;
		double pEnotF= 1.0;
		double pnotEF = 1.0;
		double pnotEnotF = 1.0;
		//Note that if the counts are zero the log will be multiplied by zero so it doesn't matter
		//Also note the pEF etc. are the RATIOS p(E|F)/p(E) etc. not the acutal probabilities

		if (countEF > 0) pEF = cEF * (numS/(cE*cF));
		if (countEnotF > 0) pEnotF = (countEnotF * numS)/(cE * countnotF);
		if (countnotEF > 0) pnotEF = (countnotEF * numS)/(countnotE * cF);
		if (countnotEnotF > 0) pnotEnotF = (countnotEnotF * numS) /(countnotE*countnotF);
		double llrScore = cEF * Math.log(pEF) + countEnotF * Math.log(pEnotF) + countnotEF * Math.log(pnotEF) + countnotEnotF * Math.log(pnotEnotF);
		if(llrScore < 0.0) llrScore = 0.0;
		return llrScore;
	}
	
	static double lgamma(double x) {
		if (x < 0.0) return lgamma(1.0);
		double tmp = (x - 0.5) * Math.log(x + 4.5) - (x + 4.5);
		double ser = 1.0 + 76.18009173    / (x + 0)   - 86.50532033    / (x + 1)
	    	+ 24.01409822    / (x + 2)   -  1.231739516   / (x + 3)
	    	+  0.00120858003 / (x + 4)   -  0.00000536382 / (x + 5);
		final double res = tmp + Math.log(ser * Math.sqrt(2 * Math.PI));
		System.err.println("lg("+x+")="+res);
		return res;
	}

	public static double fishersExact(int countEF, int countE, int countF, int N) {
	    double a = countEF;
	    double b = (countF - countEF);
	    double c = (countE - countEF);
	    double d = (N - countE - countF + countEF);
	    double n = a + b + c + d;
	    
	    double xp = lgamma(1.0+a+c) + lgamma(1.0+b+d) + lgamma(1.0+a+b) + lgamma(1.0+c+d) - lgamma(1.0+n) - lgamma(1.0+a) - lgamma(1.0+b) - lgamma(1.0+c) - lgamma(1.0+d);
	    System.err.println("xp="+xp);
	    double cp = Math.exp(lgamma(1.0+a+c) + lgamma(1.0+b+d) + lgamma(1.0+a+b) + lgamma(1.0+c+d) - lgamma(1.0+n) - lgamma(1.0+a) - lgamma(1.0+b) - lgamma(1.0+c) - lgamma(1.0+d));
	    System.err.println("cp="+cp);
	    double total_p = 0.0;
	    int tc = Math.min((int)b,(int)c);
	    for (int i=0; i<=tc; i++) {
	      total_p += cp;
	      double coef = b*c/(a+1.0)/(d+1.0);
	      cp *= coef;
	      ++a; --c; ++d; --b;
	    }
	  return total_p;

	}

	
}