import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Uses a sum-of-uniform-0-1 random variable approximation to the U statistic in order to return an approximate
 * p-value. See Buckle, Kraft, van Eeden [1969] (approx) and Billingsly [1995] or Stephens, MA [1966, biometrika] (sum of uniform CDF)
 *
 * @param n - The number of entries in the stochastically smaller (dominant) set
 * @param m - The number of entries in the stochastically larger (dominated) set
 * @param u - mann-whitney u value
 * @return p-value according to sum of uniform approx
 * todo -- this is currently not called due to not having a good characterization of where it is significantly more accurate than the
 * todo -- normal approxmation (e.g. enough to merit the runtime hit)
 */
public static double calculatePUniformApproximation(int n, int m, long u) {
    long R = u + (n * (n + 1)) / 2;
    double a = Math.sqrt(m * (n + m + 1));
    double b = (n / 2.0) * (1 - Math.sqrt((n + m + 1) / m));
    double z = b + ((double) R) / a;
    if (z < 0) {
        return 1.0;
    } else if (z > n) {
        return 0.0;
    } else {
        if (z > ((double) n) / 2) {
            return 1.0 - 1 / (Arithmetic.factorial(n)) * uniformSumHelper(z, (int) Math.floor(z), n, 0);
        } else {
            return 1 / (Arithmetic.factorial(n)) * uniformSumHelper(z, (int) Math.floor(z), n, 0);
        }
    }
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
private static double computePValue(int[][] table) {

        int[] rowSums = {sumRow(table, 0), sumRow(table, 1)};
        int[] colSums = {sumColumn(table, 0), sumColumn(table, 1)};
        int N = rowSums[0] + rowSums[1];

        // calculate in log space so we don't die with high numbers
        double pCutoff = Arithmetic.logFactorial(rowSums[0])
                + Arithmetic.logFactorial(rowSums[1])
                + Arithmetic.logFactorial(colSums[0])
                + Arithmetic.logFactorial(colSums[1])
                - Arithmetic.logFactorial(table[0][0])
                - Arithmetic.logFactorial(table[0][1])
                - Arithmetic.logFactorial(table[1][0])
                - Arithmetic.logFactorial(table[1][1])
                - Arithmetic.logFactorial(N);
        return Math.exp(pCutoff);
    }
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
private static double computePValue(int[][] table) {

        int[] rowSums = { sumRow(table, 0), sumRow(table, 1) };
        int[] colSums = { sumColumn(table, 0), sumColumn(table, 1) };
        int N = rowSums[0] + rowSums[1];

        // calculate in log space so we don't die with high numbers
        double pCutoff = Arithmetic.logFactorial(rowSums[0])
                         + Arithmetic.logFactorial(rowSums[1])
                         + Arithmetic.logFactorial(colSums[0])
                         + Arithmetic.logFactorial(colSums[1])
                         - Arithmetic.logFactorial(table[0][0])
                         - Arithmetic.logFactorial(table[0][1])
                         - Arithmetic.logFactorial(table[1][0])
                         - Arithmetic.logFactorial(table[1][1])
                         - Arithmetic.logFactorial(N);
        return Math.exp(pCutoff);
    }
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
public static double computePValue(int[][] table) {

        int[] rowSums = { sumRow(table, 0), sumRow(table, 1) };
        int[] colSums = { sumColumn(table, 0), sumColumn(table, 1) };
        int N = rowSums[0] + rowSums[1];

        // calculate in log space for better precision
        double pCutoff = Arithmetic.logFactorial(rowSums[0])
                + Arithmetic.logFactorial(rowSums[1])
                + Arithmetic.logFactorial(colSums[0])
                + Arithmetic.logFactorial(colSums[1])
                - Arithmetic.logFactorial(table[0][0])
                - Arithmetic.logFactorial(table[0][1])
                - Arithmetic.logFactorial(table[1][0])
                - Arithmetic.logFactorial(table[1][1])
                - Arithmetic.logFactorial(N);
        return Math.exp(pCutoff);
    }
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Uses a sum-of-uniform-0-1 random variable approximation to the U statistic in order to return an approximate
 * p-value. See Buckle, Kraft, van Eeden [1969] (approx) and Billingsly [1995] or Stephens, MA [1966, biometrika] (sum of uniform CDF)
 * @param n - The number of entries in the stochastically smaller (dominant) set
 * @param m - The number of entries in the stochastically larger (dominated) set
 * @param u - mann-whitney u value
 * @return p-value according to sum of uniform approx
 * todo -- this is currently not called due to not having a good characterization of where it is significantly more accurate than the
 * todo -- normal approxmation (e.g. enough to merit the runtime hit)
 */
public static double calculatePUniformApproximation(int n, int m, long u) {
    long R = u + (n*(n+1))/2;
    double a = Math.sqrt(m*(n+m+1));
    double b = (n/2.0)*(1-Math.sqrt((n+m+1)/m));
    double z = b + ((double)R)/a;
    if ( z < 0 ) { return 1.0; }
    else if ( z > n ) { return 0.0; }
    else {
        if ( z > ((double) n) /2 ) {
            return 1.0-1/(Arithmetic.factorial(n))*uniformSumHelper(z, (int) Math.floor(z), n, 0);
        } else {
            return 1/(Arithmetic.factorial(n))*uniformSumHelper(z, (int) Math.floor(z), n, 0);
        }
    }
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
public static double getBinomialCoeff(int n, int k) throws ExtensionException {
  // Returns "n choose k" as a double. Note the "integerization" of
  // the double return value.
  try {
  return Math.rint(Arithmetic.binomial((long) n, (long) k));
  } catch (IllegalArgumentException | ArithmeticException ex) {
    throw new ExtensionException("colt .Arithmetic.binomial reports: " + ex);
  }
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Constructs an approximate quantile finder with b buffers, each having k elements.
 * @param b the number of buffers
 * @param k the number of elements per buffer
 * @param N the total number of elements over which quantiles are to be computed.
 * @param samplingRate 1.0 --> all elements are consumed. 10.0 --> Consumes one random element from successive blocks of 10 elements each. Etc.
 * @param generator a uniform random number generator.
 */
public KnownDoubleQuantileEstimator(int b, int k, long N, double samplingRate, RandomEngine generator) {
	this.samplingRate = samplingRate;
	this.N = N;

	if (this.samplingRate <= 1.0) {
		this.samplingAssistant = null;
	}
	else {
		this.samplingAssistant = new RandomSamplingAssistant(Arithmetic.floor(N/samplingRate), N, generator);
	}
	
	setUp(b,k);
	this.clear();
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Removes all elements from the receiver.  The receiver will
 * be empty after this call returns, and its memory requirements will be close to zero.
 */
public void clear() {
	super.clear();
	this.beta=1.0;
	this.weHadMoreThanOneEmptyBuffer = false;
	//this.setSamplingRate(samplingRate,N);

	RandomSamplingAssistant assist = this.samplingAssistant;
	if (assist != null) {
		this.samplingAssistant = new RandomSamplingAssistant(Arithmetic.floor(N/samplingRate), N, assist.getRandomGenerator());
	}
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Sets the parameters number of trials and the probability of success.
 * @param n the number of trials
 * @param p the probability of success.
 * @throws IllegalArgumentException if <tt>n*Math.min(p,1-p) &lt;= 0.0</tt>
 */
public void setNandP(int n, double p) {
	if (n*Math.min(p,1-p) <= 0.0) throw new IllegalArgumentException();
	this.n = n;
	this.p = p;
	
	this.log_p = Math.log(p);
	this.log_q = Math.log(1.0-p);
	this.log_n = Arithmetic.logFactorial(n);
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Returns the probability distribution function.
 */
public double pdf(int k) {
	return Math.exp(k*Math.log(this.mean) - Arithmetic.logFactorial(k) - this.mean);
	
	// Overflow sensitive:
	// return (Math.pow(mean,k) / cephes.Arithmetic.factorial(k)) * Math.exp(-this.mean);
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Computes the number of buffers and number of values per buffer such that
 * quantiles can be determined with an approximation error no more than epsilon with a certain probability.
 * Assumes that quantiles are to be computed over N values.
 * The required sampling rate is computed and stored in the first element of the provided <tt>returnSamplingRate</tt> array, which, therefore must be at least of length 1.
 * @param N the anticipated number of values over which quantiles shall be computed (e.g 10^6).
 * @param epsilon the approximation error which is guaranteed not to be exceeded (e.g. <tt>0.001</tt>) (<tt>0 &lt;= epsilon &lt;= 1</tt>). To get exact result, set <tt>epsilon=0.0</tt>;
 * @param delta the probability that the approximation error is more than than epsilon (e.g. <tt>0.0001</tt>) (<tt>0 &lt;= delta &lt;= 1</tt>). To avoid probabilistic answers, set <tt>delta=0.0</tt>.
 * @param quantiles the number of quantiles to be computed (e.g. <tt>100</tt>) (<tt>quantiles &gt;= 1</tt>). If unknown in advance, set this number large, e.g. <tt>quantiles &gt;= 10000</tt>.
 * @param samplingRate a <tt>double[1]</tt> where the sampling rate is to be filled in.
 * @return <tt>long[2]</tt> - <tt>long[0]</tt>=the number of buffers, <tt>long[1]</tt>=the number of elements per buffer, <tt>returnSamplingRate[0]</tt>=the required sampling rate.
 */
protected static long[] known_N_compute_B_and_K_slow(long N, double epsilon, double delta, int quantiles, double[] returnSamplingRate) {
	final int maxBuffers = 50;
	final int maxHeight = 50;
	final double N_double = N;

	// One possibility is to use one buffer of size N
	//
	long ret_b = 1;
	long ret_k = N;
	double sampling_rate = 1.0;
	long memory = N;


	// Otherwise, there are at least two buffers (b >= 2)
	// and the height of the tree is at least three (h >= 3)
	//
	// We restrict the search for b and h to MAX_BINOM, a large enough value for
	// practical values of    epsilon >= 0.001   and    delta >= 0.00001
	//
	final double logarithm = Math.log(2.0*quantiles/delta);
	final double c = 2.0 * epsilon * N_double;
	for (long b=2 ; b<maxBuffers ; b++)
		for (long h=3 ; h<maxHeight ; h++) {
			double binomial = Arithmetic.binomial(b+h-2, h-1);
			long tmp = (long) Math.ceil(N_double / binomial);
			if ((b * tmp < memory) && 
					((h-2) * binomial - Arithmetic.binomial(b+h-3, h-3) + Arithmetic.binomial(b+h-3, h-2)
					<= c) ) {
				ret_k = tmp ;
				ret_b = b ;
				memory = ret_k * b;
				sampling_rate = 1.0 ;
			}
			if (delta > 0.0) {
				double t = (h-2) * Arithmetic.binomial(b+h-2, h-1) - Arithmetic.binomial(b+h-3, h-3) + Arithmetic.binomial(b+h-3, h-2) ;
				double u = logarithm / epsilon ;
				double v = Arithmetic.binomial (b+h-2, h-1) ;
				double w = logarithm / (2.0*epsilon*epsilon) ;

				// From our SIGMOD 98 paper, we have two equantions to satisfy:
				// t  <= u * alpha/(1-alpha)^2
				// kv >= w/(1-alpha)^2
				//
				// Denoting 1/(1-alpha)    by x,
				// we see that the first inequality is equivalent to
				// t/u <= x^2 - x
				// which is satisfied by x >= 0.5 + 0.5 * sqrt (1 + 4t/u)
				// Plugging in this value into second equation yields
				// k >= wx^2/v

				double x = 0.5 + 0.5 * Math.sqrt(1.0 + 4.0*t/u) ;
				long k = (long) Math.ceil(w*x*x/v) ;
				if (b * k < memory) {
					ret_k = k ;
					ret_b = b ;
					memory = b * k ;
					sampling_rate = N_double*2.0*epsilon*epsilon / logarithm ;
				}
			}
		}
		
	long[] result = new long[2];
	result[0]=ret_b;
	result[1]=ret_k;
	returnSamplingRate[0]=sampling_rate;
	return result;
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
private static double fc_lnpk(int k, int N_Mn, int M, int n) {
	return(Arithmetic.logFactorial(k) + Arithmetic.logFactorial(M - k) + Arithmetic.logFactorial(n - k) + Arithmetic.logFactorial(N_Mn + k));
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Returns a random number from the distribution.
 */
protected int hmdu(int N, int M, int n, RandomEngine randomGenerator) {

  int            I, K;
  double              p, nu, c, d, U;

  if (N != N_last || M != M_last || n != n_last) {   // set-up           */
		N_last = N;
	 M_last = M;
	 n_last = n;

	 Mp = (double) (M + 1);
	 np = (double) (n + 1);  N_Mn = N - M - n;

	 p  = Mp / (N + 2.0);
	 nu = np * p;                             /* mode, real       */
	 if ((m = (int) nu) == nu && p == 0.5) {     /* mode, integer    */
		mp = m--;
	 }
	 else {
		mp = m + 1;                           /* mp = m + 1       */
		}

 /* mode probability, using the external function flogfak(k) = ln(k!)    */
	 fm = Math.exp(Arithmetic.logFactorial(N - M) - Arithmetic.logFactorial(N_Mn + m) - Arithmetic.logFactorial(n - m)
		 + Arithmetic.logFactorial(M)     - Arithmetic.logFactorial(M - m)    - Arithmetic.logFactorial(m)
		 - Arithmetic.logFactorial(N)     + Arithmetic.logFactorial(N - n)    + Arithmetic.logFactorial(n)   );

 /* safety bound  -  guarantees at least 17 significant decimal digits   */
 /*                  b = min(n, (long int)(nu + k*c')) */
		b = (int) (nu + 11.0 * Math.sqrt(nu * (1.0 - p) * (1.0 - n/(double)N) + 1.0));	
		if (b > n) b = n;
	}

	for (;;) {
		if ((U = randomGenerator.raw() - fm) <= 0.0)  return(m);
		c = d = fm;

 /* down- and upward search from the mode                                */
		for (I = 1; I <= m; I++) {
			K  = mp - I;                                   /* downward search  */
			c *= (double)K/(np - K) * ((double)(N_Mn + K)/(Mp - K));
			if ((U -= c) <= 0.0)  return(K - 1);

			K  = m + I;                                    /* upward search    */
			d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
			if ((U -= d) <= 0.0)  return(K);
		}

 /* upward search from K = 2m + 1 to K = b                               */
		for (K = mp + m; K <= b; K++) {
			d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
			if ((U -= d) <= 0.0)  return(K);
		}
	}
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Returns the probability distribution function.
 */
public double pdf(int k) {
	return Arithmetic.binomial(my_s, k) * Arithmetic.binomial(my_N - my_s, my_n - k) 
		/ Arithmetic.binomial(my_N, my_n);
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Returns the probability distribution function.
 */
public double pdf(int k) {
	if (k < 0) throw new IllegalArgumentException();
	int r = this.n - k;
	return Math.exp(this.log_n - Arithmetic.logFactorial(k) - Arithmetic.logFactorial(r) + this.log_p * k + this.log_q * r);
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
private static double f(int k, double l_nu, double c_pm) {
	return  Math.exp(k * l_nu - Arithmetic.logFactorial(k) - c_pm);
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Helper function for the sum of n uniform random variables
 *
 * @param z - value at which to compute the (un-normalized) cdf
 * @param m - a cutoff integer (defined by m <= z < m + 1)
 * @param n - the number of uniform random variables
 * @param k - holder variable for the recursion (alternatively, the index of the term in the sequence)
 * @return the (un-normalized) cdf for the sum of n random variables
 */
private static double uniformSumHelper(double z, int m, int n, int k) {
    if (k > m) {
        return 0;
    }
    int coef = (k % 2 == 0) ? 1 : -1;
    return coef * Arithmetic.binomial(n, k) * Math.pow(z - k, n) + uniformSumHelper(z, m, n, k + 1);
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Computes the log-space PDF of the Poisson distribution with the
 * specified mean. Code lifted from the COLT Poisson class so that it can
 * be implemented as a static method instead of having to instantiate a
 * COLT Poisson object
 * @param mean
 * @param k
 * @return
 */
public static double poissonLogPDF(double mean, int k) {
  return k * Math.log(mean) - Arithmetic.logFactorial(k) - mean;
}
 

import cern.jet.math.Arithmetic; //导入依赖的package包/类
/**
 * Helper function for the sum of n uniform random variables
 * @param z - value at which to compute the (un-normalized) cdf
 * @param m - a cutoff integer (defined by m <= z < m + 1)
 * @param n - the number of uniform random variables
 * @param k - holder variable for the recursion (alternatively, the index of the term in the sequence)
 * @return the (un-normalized) cdf for the sum of n random variables
 */
private static double uniformSumHelper(double z, int m, int n, int k) {
    if ( k > m ) { return 0; }
    int coef = (k % 2 == 0) ? 1 : -1;
    return coef*Arithmetic.binomial(n,k)*Math.pow(z-k,n) + uniformSumHelper(z,m,n,k+1);
}