gt: 	dn = n - k;
context: /**
   * papaya: A collection of utilities for Statistics and Matrix-related manipulations
   * http://adilapapaya.com/papayastatistics.html, 1
   * Created by Adila Faruk, http://adilapapaya.com, May 1st 2012
   *
   *
   * Copyright (C) 2012 Adila Faruk http://adilapapaya.com 
   * Copyright 1999 CERN - European Organization for Nuclear Research.
   * Copyright 1995 Stephen L.Moshier.
   * 
   * 
   * This library is free software; you can redistribute it and/or
   * modify it under the terms of the GNU Lesser General Public
   * License as published by the Free Software Foundation; either
   * version 2.1 of the License, or (at your option) any later version.
   * 
   * This library 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
   * Lesser General Public License for more details.
   * 
   * You should have received a copy of the GNU Lesser General
   * Public License along with this library; if not, write to the
   * Free Software Foundation, Inc., 59 Temple Place, Suite 330,
   * Boston, MA  02111-1307  USA
   */
  package papaya;
  /**
   * Cumulative distribution functions and corresponding inverses of certain probability distributions.
   * <p>
   * <b>Description:</b>
   * <br><code>...cdf</code> signifies the cumulative distribution function and is similar to 
   * the <code>q...</code> functions in <code>R</code> and the <code>...cdf</code> functions in 
   * <code>MATLAB</code>. The cdf represent the integral, from -infinity to <code>0</code> of the
   * corresponding continuous probability density function (pdf), or, the sum from <code>0</code> to </code>x</code>
   * of an associated finite pdf.</br>
   * 
   * <br><code>...inv</code> signifies the inverse of the cumulative distribution function and is similar to 
   * the <code>p...</code> functions in <code>R</code> and the <code>...inv</code> functions in 
   * <code>MATLAB</code>.</br>
   * 
   * <p><b>Implementation:</b>
   * <ul>
   * <li>Direct application examples: 
   * <br><code>normcdf(1.96) = 0.975</code> --> 97.5% of the curve is within +/-1.96 of 0.
   * <br><code>norminv(0.975) = 1.96</code> --> finds the point x (=1.96) such that <code>normcdf(x) = .975</code>.
   * <br><code>norminv(  normcdf(x) ) = x </code> 
   * <li>Significance tests:
   * <br>Two-tailed significance test with <code>alpha = 0.05</code> and <code>H_0: &mu;=0</code>
   * --> <code>p = 1-alpha/2 = 0.975</code>
   * --> Reject null hypothesis <code>&mu;=0</code> if
   * the test statistic is greater or less than than +/-1.96 <code>(= +/- norminv(p) )</code>.
   * 
   * <br>One-tailed significance test with <code>alpha = 0.05</code> and <code>H_0: &mu; &lt; 0</code>
   * --> <code>p = 1-alpha = 0.95</code>
   * --> Reject null hypothesis <code>&mu; &lt; 0</code> if
   * the test statistic is greater than 1.645 <code>(= norminv(p) )</code>.
   * <br> One-tailed significance test with <code>alpha = 0.05</code> and <code>H_0: &mu; &gt; 0</code>
   * --> <code>p = 1-alpha = 0.95</code>
   * --> Reject null hypothesis <code>&mu; &gt; 0</code> if
   * the test statistic is less than -1.645 <code>(= -norminv(p) )</code>.
   * </ul>
   * </p>
   * 
   * 
   * The code is partially adapted from the CERN Jet Java libraries, which in turn was adapted from the 
   * Cephes Mathematical Library. As far a as I can tell, Stephen L. Moshier wrote the original C++ code for 
   * CEPHES (1989, moshier@na-net.ornl.gov),  Gedeck (at Novartis) and Wolfgang Hoschek (at CERN, hats off to you mate!) 
   * then adapted it for the Java platform, making some significant changes along the way. 
   * <br>
   * 2012 rolls along, and Adila (hello!) finds the library and makes even more changes, increasing the
   * probability (ha.ha.) of things going massively wrong (pValue<.0001).
   *
   */
  
  //public class Probability extends cern.jet.math.Constants {
  public class Probability implements PapayaConstants {
  	/*************************************************
  	 *    COEFFICIENTS FOR METHOD  normalInverse()   *
  	 *************************************************/
  	/* approximation for 0 <= |y - 0.5| <= 3/8 */
  	protected static final double P0[] = {
  		-5.99633501014107895267E1,
  		9.80010754185999661536E1,
  		-5.66762857469070293439E1,
  		1.39312609387279679503E1,
  		-1.23916583867381258016E0,
  	};
  	protected static final double Q0[] = {
  		1.00000000000000000000E0,
  		1.95448858338141759834E0,
  		4.67627912898881538453E0,
  		8.63602421390890590575E1,
  		-2.25462687854119370527E2,
  		2.00260212380060660359E2,
  		-8.20372256168333339912E1,
  		1.59056225126211695515E1,
  		-1.18331621121330003142E0,
  	};
  
  
  	/* Approximation for interval z = sqrt(-2 log y ) for z between 2 and 8
  	 * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
  	 */
  	protected static final double P1[] = {
  		4.05544892305962419923E0,
  		3.15251094599893866154E1,
  		5.71628192246421288162E1,
  		4.40805073893200834700E1,
  		1.46849561928858024014E1,
  		2.18663306850790267539E0,
  		-1.40256079171354495875E-1,
  		-3.50424626827848203418E-2,
  		-8.57456785154685413611E-4,
  	};
  	protected static final double Q1[] = {
  		1.00000000000000000000E0,
  		1.57799883256466749731E1,
  		4.53907635128879210584E1,
  		4.13172038254672030440E1,
  		1.50425385692907503408E1,
  		2.50464946208309415979E0,
  		-1.42182922854787788574E-1,
  		-3.80806407691578277194E-2,
  		-9.33259480895457427372E-4,
  	};
  
  	/* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
  	 * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
  	 */
  	protected static final double  P2[] = {
  		3.23774891776946035970E0,
  		6.91522889068984211695E0,
  		3.93881025292474443415E0,
  		1.33303460815807542389E0,
  		2.01485389549179081538E-1,
  		1.23716634817820021358E-2,
  		3.01581553508235416007E-4,
  		2.65806974686737550832E-6,
  		6.23974539184983293730E-9,
  	};
  	protected static final double  Q2[] = {
  		1.00000000000000000000E0,
  		6.02427039364742014255E0,
  		3.67983563856160859403E0,
  		1.37702099489081330271E0,
  		2.16236993594496635890E-1,
  		1.34204006088543189037E-2,
  		3.28014464682127739104E-4,
  		2.89247864745380683936E-6,
  		6.79019408009981274425E-9,
  	};
  
  	/**
  	 * Makes this class non instantiable, but still let's others inherit from it.
  	 */
  	protected Probability() {}
  	/**
  	 * Returns the area from zero to <code>x</code> under the beta density
  	 * function.
  	 * <pre>
  	 * 
  	 * P(x) = &Gamma;(a+b) / [ &Gamma;(a) * &Gamma;(b)  ] * Integral_(0 to x) [  t^(a-1) * (1-t)^(b-1) ] dt
  	 *
  	 * </pre>
  	 * This function is identical to the incomplete beta
  	 * integral function 
  	 * 
  	 * <code>Gamma.incompleteBeta(a, b, x)</code>.
  	 *
  	 * The complemented function is {@link betacdfComplemented}, specified by
  	 *
  	 * <code>1 - P(1-x)  =  Gamma.incompleteBeta( b, a, x )</code>;
  	 *
  	 */
  	static public double betacdf(double x, double a, double b) {
  		return Gamma.incompleteBeta( a, b, x );
  	}
  	/**
  	 * Returns the area under the right hand tail (from <code>x</code> to
  	 * infinity) of the beta density function.
  	 * 
  	 * This function is identical to the incomplete beta
  	 * integral function <code>Gamma.incompleteBeta(b, a, x)</code>.
  	 */
  	static public double betacdfComplemented(double x, double a, double b ) {
  		return Gamma.incompleteBeta( b, a, x );
  	}
  	
  	/** 
  	 * Returns the inverse of the beta cumulative distribution function.
  	 * That is, given <code>p</code>, the function finds <code>x</code> such that
  	 *	<pre>
  	 *  Gamma.incompleteBeta( a, b, x ) = beta(x, a, b) = p;
  	 *</pre>
  	 *
  	 * This function is identical to the incomplete beta inverse
  	 * function <code>Gamma.incompleteBetaInverse(aa, bb, yy0)</code>.
  	 *
  	 * @param a the alpha parameter of the beta distribution.
  	 * @param b the beta parameter of the beta distribution.
  	 * @param p the value for which to solve for the corresponding <code>x</code> in the 
  	 * incomplete Beta integral.
  	 * @return the value <code>x</code> such that <code>beta( a, b, x ) = p</code>.
  	 */
  	static public double betainv( double p, double a, double b) {
  		return Gamma.incompleteBetaInverse( a, b, p);
  	}
  	/**
  	 * Returns the sum of the terms <code>0</code> through <code>k</code> of the Binomial
  	 * probability density.
  	 * <pre>
  	 * 
  	 * Sum_(j=0 to k)  [ (n!)/(j!*(n-j)! ]  *  [ p^j * (1-p)^(n-j) ]
  	 *
  	 * </pre>
  	 * The terms are not summed directly; instead the incomplete
  	 * beta integral is employed, according to the formula
  	 * <p>
  	 * <code>y = binomcdf( p, k, n) = Gamma.incompleteBeta( n-k, k+1, 1-p )</code>.
  	 * <p>
  	 * All arguments must be positive, 
  	 * @param k end term.
  	 * @param n the number of trials.
  	 * @param p the probability of success (must be in <code>(0.0,1.0)</code>).
  	 */
  	static public double binomcdf(double p, int k, int n) {
  		if( (p < 0.0) || (p > 1.0) ) throw new IllegalArgumentException();
  		if( (k < 0) || (n < k) ) throw new IllegalArgumentException();
  
  		if( k == n ) return( 1.0 );
  		if( k == 0 ) return Math.pow( 1.0-p, n-k );
  
  		return Gamma.incompleteBeta( n-k, k+1, 1.0 - p );
  	}
  	/**
  	 * Returns the sum of the terms <code>k+1</code> through <code>n</code> of the Binomial
  	 * probability density.
  	 * <pre>
  	 * 
  	 * Sum_(j=k+1 to n)  [ (n!)/(j!*(n-j)! ]  *  [ p^j * (1-p)^(n-j) ]
  	 *
  	 * </pre>
  	 * The terms are not summed directly; instead the incomplete
  	 * beta integral is employed, according to the formula
  	 * <p>
  	 * <code>y = binomcdfComplemented( p, k, n ) = Gamma.incompleteBeta( k+1, n-k, p )</code>.
  	 * <p>
  	 * All arguments must be positive, 
  	 * @param k (start-1) term.
  	 * @param n the number of trials.
  	 * @param p the probability of success (must be in <code>(0.0,1.0)</code>).
  	 */
  	static public double binomcdfComplemented( double p , int k, int n) {
  		if( (p < 0.0) || (p > 1.0) ) throw new IllegalArgumentException();
  		if( (k < 0) || (n < k) ) throw new IllegalArgumentException();
  
  		if( k == n ) return( 0.0 );
  		if( k == 0 ) return 1.0 - Math.pow( 1.0-p, n-k );
  
  		return Gamma.incompleteBeta( k+1, n-k, p );
  	}
  	
  	/** 
  	 * Finds the event probability p such that the sum of the
  	 * terms 0 through k of the Binomial probability density
  	 * is equal to the given cumulative probability y. That is,
  	 *<pre>
  	 * binomcdf( p, k, n ) = y;
  	 *</pre>
  	 * This is accomplished using the inverse beta integral
  	 * function and the relation
  	 * <p>
  	 * 1 - p = Gamma.incompleteBetaInverse( n-k, k+1, y ).
  	 * </p>
  	 * All arguments must be positive.
  	 * @param k end term.
  	 * @param n the number of trials.
  	 * @param y the value to solve for such that <code>binomcdf( p, k, n ) = y</code>.
  	 */	
  	static public double binominv( double y, int k, int n){
  	double dk, dn, p;
  
  	if( (y < 0.0) || (y > 1.0) ) throw new IllegalArgumentException("Test statistic not within possible range.");
  	if( (k < 0) || (n < k) ) throw new IllegalArgumentException();
  
  