The natural logarithm of the beta function.

@param p require p>0
@param q require q>0
@return 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
@author Jaco van Kooten
352337
private static function _logBeta($p, $q) {
		if ($p != self::$_logBetaCache_p || $q != self::$_logBetaCache_q) {
			self::$_logBetaCache_p = $p;
			self::$_logBetaCache_q = $q;
			if (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > LOG_GAMMA_X_MAX_VALUE)) {
				self::$_logBetaCache_result = 0.0;
			} else {
				self::$_logBetaCache_result = self::_logGamma($p) + self::_logGamma($q) - self::_logGamma($p + $q);
			}
		}
		return self::$_logBetaCache_result;
	}
883389
protected static double regularizedIncBetaQuadrature(double alpha, double beta, double x) {
    final double alphapbeta = alpha + beta;
    final double a1 = alpha - 1.0;
    final double b1 = beta - 1.0;
    final double mu = alpha / alphapbeta;
    final double lnmu = FastMath.log(mu);
    final double lnmuc = FastMath.log1p(-mu);
    double t = FastMath.sqrt(alpha * beta / (alphapbeta * alphapbeta * (alphapbeta + 1.0)));
    final double xu;
    if(x > alpha / alphapbeta) {
      if(x >= 1.0) {
        return 1.0;
      }
      xu = Math.min(1.0, Math.max(mu + 10.0 * t, x + 5.0 * t));
    }
    else {
      if(x <= 0.0) {
        return 0.0;
      }
      xu = Math.max(0.0, Math.min(mu - 10.0 * t, x - 5.0 * t));
    }
    double sum = 0.0;
    for(int i = 0; i < GAUSSLEGENDRE_Y.length; i++) {
      t = x + (xu - x) * GAUSSLEGENDRE_Y[i];
      sum += GAUSSLEGENDRE_W[i] * FastMath.exp(a1 * (FastMath.log(t) - lnmu) + b1 * (FastMath.log1p(-t) - lnmuc));
    }
    double ans = sum * (xu - x) * FastMath.exp(a1 * lnmu - GammaDistribution.logGamma(alpha) + b1 * lnmuc - GammaDistribution.logGamma(beta) + GammaDistribution.logGamma(alphapbeta));
    return ans > 0 ? 1.0 - ans : -ans;
  }
650585
public static double KullbackLeiblerDivergence(double[] p, double[] q) {
        boolean intersection = false;
        double k = 0;

        for (int i = 0; i < p.length; i++) {
            if (p[i] != 0 && q[i] != 0) {
                intersection = true;
                k += p[i] * Math.log(p[i] / q[i]);
            }
        }

        if (intersection)
            return k;
        else
            return Double.POSITIVE_INFINITY;
    }
941654
def BetaPrime(alpha, beta, tag=None):
    
    assert (
        alpha > 0 and beta > 0
    ), 'BetaPrime "alpha" and "beta" parameters must be greater than zero'
    x = Beta(alpha, beta, tag)
    return x / (1 - x)
437930
private static boolean needPlural(final int n)
    {
        int tmp = 0;
        int an = Math.abs(n);

        while (an > 0)
        {
            tmp = an % 10;

            if (tmp > 1)
            {
                return true;
            }

            an /= 10;
        }

        return false;
    }
870485
public static double Correlation(double[] p, double[] q) {

        double x = 0;
        double y = 0;

        for (int i = 0; i < p.length; i++) {
            x += -p[i];
            y += -q[i];
        }

        x /= p.length;
        y /= q.length;

        double num = 0;
        double den1 = 0;
        double den2 = 0;
        for (int i = 0; i < p.length; i++) {
            num += (p[i] + x) * (q[i] + y);

            den1 += Math.abs(Math.pow(p[i] + x, 2));
            den2 += Math.abs(Math.pow(q[i] + x, 2));
        }

        return 1 - (num / (Math.sqrt(den1) * Math.sqrt(den2)));

    }
276826
function invErfc(p) {
  if (p < 0.0 || p > 2.0) {
    throw RangeError('Argument must be betweeen 0 and 2');
  }

  else if (p === 0.0) {
    return Infinity;
  }
  
  else if (p === 2.0) {
    return -Infinity;
  }
  
  else {
    var pp = p < 1.0 ? p : 2.0 - p;
    var t = Math.sqrt(-2.0 * Math.log(pp / 2.0));
    var x = -0.70711 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t);

    var err1 = erfc(x) - pp;
    x += err1 / (M_2_SQRTPI * Math.exp(-x * x) - x * err1);
    var err2 = erfc(x) - pp;
    x += err2 / (M_2_SQRTPI * Math.exp(-x * x) - x * err2);

    return p < 1.0 ? x : -x;
  }
}
1286929
double betaFunction(List<Double> x, int type) {
    double beta;
    beta = 0;
    int dim = x.size();

    if (dim == 0) {
      beta = 0;
    }

    if (type == 1) {
      beta = 0;
      for (int i = 0; i < dim; i++) {
        beta += x.get(i) * x.get(i);
      }
      beta = 2.0 * beta / dim;
    }

    if (type == 2) {
      beta = 0;
      for (int i = 0; i < dim; i++) {
        beta += Math.sqrt(i + 1) * x.get(i) * x.get(i);
      }
      beta = 2.0 * beta / dim;
    }

    if (type == 3) {
      double sum = 0, xx;
      for (int i = 0; i < dim; i++) {
        xx = 2 * x.get(i);
        sum += (xx * xx - Math.cos(4 * Math.PI * xx) + 1);
      }
      beta = 2.0 * sum / dim;
    }

    if (type == 4) {
      double sum = 0, prod = 1, xx;
      for (int i = 0; i < dim; i++) {
        xx = 2 * x.get(i);
        sum += xx * xx;
        prod *= Math.cos(10 * Math.PI * xx / Math.sqrt(i + 1));
      }
      beta = 2.0 * (sum - 2 * prod + 2) / dim;
    }

    return beta;
  }
842966
@GwtIncompatible("TODO")
  public static BigInteger divide(BigInteger p, BigInteger q, RoundingMode mode) {
    BigDecimal pDec = new BigDecimal(p);
    BigDecimal qDec = new BigDecimal(q);
    return pDec.divide(qDec, 0, mode).toBigIntegerExact();
  }
352148
function qrPolynomial(num, shift) {

		if (typeof num.length == 'undefined') {
			throw new Error(num.length + '/' + shift);
		}

		var _num = function() {
			var offset = 0;
			while (offset < num.length && num[offset] == 0) {
				offset += 1;
			}
			var _num = new Array(num.length - offset + shift);
			for (var i = 0; i < num.length - offset; i += 1) {
				_num[i] = num[i + offset];
			}
			return _num;
		}();

		var _this = {};

		_this.getAt = function(index) {
			return _num[index];
		};

		_this.getLength = function() {
			return _num.length;
		};

		_this.multiply = function(e) {

			var num = new Array(_this.getLength() + e.getLength() - 1);

			for (var i = 0; i < _this.getLength(); i += 1) {
				for (var j = 0; j < e.getLength(); j += 1) {
					num[i + j] ^= QRMath.gexp(QRMath.glog(_this.getAt(i) ) + QRMath.glog(e.getAt(j) ) );
				}
			}

			return qrPolynomial(num, 0);
		};

		_this.mod = function(e) {

			if (_this.getLength() - e.getLength() < 0) {
				return _this;
			}

			var ratio = QRMath.glog(_this.getAt(0) ) - QRMath.glog(e.getAt(0) );

			var num = new Array(_this.getLength() );
			for (var i = 0; i < _this.getLength(); i += 1) {
				num[i] = _this.getAt(i);
			}

			for (var i = 0; i < e.getLength(); i += 1) {
				num[i] ^= QRMath.gexp(QRMath.glog(e.getAt(i) ) + ratio);
			}

			// recursive call
			return qrPolynomial(num, 0).mod(e);
		};

		return _this;
	}
1216712
private static int estimateQuotient(int u2, int u1, int u0, int v1, int v0)
    {
        // estimate qhat based on the first 2 digits of divisor divided by the first digit of a dividend
        long u21 = combineInts(u2, u1);
        long qhat;
        if (u2 == v1) {
            qhat = INT_BASE - 1;
        }
        else if (u21 >= 0) {
            qhat = u21 / (v1 & LONG_MASK);
        }
        else {
            qhat = divideUnsignedLong(u21, v1);
        }

        if (qhat == 0) {
            return 0;
        }

        // Check if qhat is greater than expected considering only first 3 digits of a dividend
        // This step help to eliminate all the cases when the estimation is greater than q by 2
        // and eliminates most of the cases when qhat is greater than q by 1
        //
        // u2 * b * b + u1 * b + u0 >= (v1 * b + v0) * qhat
        // u2 * b * b + u1 * b + u0 >= v1 * b * qhat + v0 * qhat
        // u2 * b * b + u1 * b - v1 * b * qhat >=  v0 * qhat - u0
        // (u21 - v1 * qhat) * b >=  v0 * qhat - u0
        // (u21 - v1 * qhat) * b + u0 >=  v0 * qhat
        // When ((u21 - v1 * qhat) * b + u0) is less than (v0 * qhat) decrease qhat by one

        int iterations = 0;
        long rhat = u21 - (v1 & LONG_MASK) * qhat;
        while (Long.compareUnsigned(rhat, INT_BASE) < 0 && Long.compareUnsigned((v0 & LONG_MASK) * qhat, combineInts(lowInt(rhat), u0)) > 0) {
            iterations++;
            qhat--;
            rhat += (v1 & LONG_MASK);
        }

        if (iterations > 2) {
            throw new IllegalStateException("qhat is greater than q by more than 2: " + iterations);
        }

        return (int) qhat;
    }
591041
public static double klDivergence(double[] p, double[] q) {
        if (p.length != q.length) {
            throw new IllegalStateException("The length of p and q must be the same.");
        }
        double delta = 1e-8;
        if (!Multinomials.isMultinomial(p, delta)) {
            throw new IllegalStateException("p is not a multinomial");
        }
        if (!Multinomials.isMultinomial(q, delta)) {
            throw new IllegalStateException("q is not a multinomial");
        }
        
        double kl = 0.0;
        for (int i=0; i<p.length; i++) {
            if (p[i] == 0.0 || q[i] == 0.0) {
                continue;
            }
            kl += p[i] * FastMath.log(p[i] / q[i]);
        }
        return kl;
    }
281119
public static double logSubtractExact(double x, double y) {
        if (x < y) {
            throw new IllegalStateException("x must be >= y. x=" + x + " y=" + y);
        }
        
        // p = 0 or q = 0, where x = log(p), y = log(q)
        if (Double.NEGATIVE_INFINITY == y) {
            return x;
        } else if (Double.NEGATIVE_INFINITY == x) {
            return y;
        } else if (x == y) {
            return Double.NEGATIVE_INFINITY;
        }
        
        // p != 0 && q != 0
        //return x + Math.log1p(-FastMath.exp(y - x));
        // 
        // The method below is more numerically stable for small differences in x and y.
        // See paper: http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf
        return x + log1mexp(x - y);
    }
1142491
@Override
  public double value(double alpha) {
    double logLikelihood = 0.0;
    for (ReadCounts rc : readCounts.values()) {
      double refProb = rc.getRefFreq();

      double pPosition = 0.0;
      for (Genotype trueGenotype1 : Genotype.values()) {
        double pGenotype1 = pGenotype(trueGenotype1, refProb);
        for (Genotype trueGenotype2 : Genotype.values()) {
          double pGenotype2 = pGenotype(trueGenotype2, refProb);

          double pObsGivenGenotype = 1.0;

          for (ReadQualityCount rqc : rc.getReadQualityCounts()) {
            Base base = rqc.getBase();
            double pErr = phredToProb(rqc.getQuality());
            double pObs
                = ((1.0 - alpha)
                    * probObsGivenTruth(base, trueGenotype1, ReadStatus.CORRECT)
                    + (alpha)
                    * probObsGivenTruth(base, trueGenotype2, ReadStatus.CORRECT)
                  ) * (1.0 - pErr)
                + ((1.0 - alpha)
                    * probObsGivenTruth(base, trueGenotype1, ReadStatus.ERROR)
                    + (alpha)
                    * probObsGivenTruth(base, trueGenotype2, ReadStatus.ERROR)
                  ) * pErr;
            pObsGivenGenotype *= Math.pow(pObs, rqc.getCount());
          }
          pPosition += pObsGivenGenotype * pGenotype1 * pGenotype2;
        }
      }
      logLikelihood += Math.log(pPosition);
    }
    return logLikelihood;
  }
661120
public static /*@pure@*/ double probToLogOdds(double prob) {

    if (gr(prob, 1) || (sm(prob, 0))) {
      throw new IllegalArgumentException("probToLogOdds: probability must " +
				     "be in [0,1] "+prob);
    }
    double p = SMALL + (1.0 - 2 * SMALL) * prob;
    return Math.log(p / (1 - p));
  }
1164797
public double quantile(double p) {
        if (nd > 0) {
            update();
        }
        int jl = 0, jh = nq - 1, j;
        while (jh - jl > 1) {
            j = (jh + jl) >> 1;
            if (p > pval[j]) {
                jl = j;
            } else {
                jh = j;
            }
        }
        j = jl;
        double q = qile[j] + (qile[j + 1] - qile[j]) * (p - pval[j]) / (pval[j + 1] - pval[j]);
        return Math.max(qile[0], Math.min(qile[nq - 1], q));
    }
519119
public static double sumInvLog1p(int s, int e) {
    double sum = 0.;
    // Iterate e + 1 .. s + 1, descending for better precision
    for(int i = e + 1; i > s; i--) {
      sum += 1. / FastMath.log(i);
    }
    return sum;
  }