Java DoubleVector.isFinite方法代码示例

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double dpois_raw(double x, double lambda, boolean give_log) {
  /*       x >= 0 ; integer for dpois(), but not e.g. for pgamma()!
  lambda >= 0
   */
  if (lambda == 0) {
    return ((x == 0) ? SignRank.R_D__1(true, give_log) : SignRank.R_D__0(true, give_log));
  }
  if (!DoubleVector.isFinite(lambda)) {
    return SignRank.R_D__0(true, give_log);
  }
  if (x < 0) {
    return (SignRank.R_D__0(true, give_log));
  }
  if (x <= lambda * Double.MIN_VALUE) {
    return (SignRank.R_D_exp(-lambda, true, give_log));
  }
  if (lambda < x * Double.MIN_VALUE) {
    return (SignRank.R_D_exp(-lambda + x * Math.log(lambda) - Gamma.logGamma(x + 1), true, give_log));
  }
  return (SignRank.R_D_fexp(2 * Math.PI * x, -Binom.stirlerr(x) - Binom.bd0(x, lambda), true, give_log));
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double bd0(double x, double np) {
  double ej, s, s1, v;
  int j;

  if (!DoubleVector.isFinite(x) || !DoubleVector.isFinite(np) || np == 0.0) {
    return DoubleVector.NaN;
  }

  if (Math.abs(x - np) < 0.1 * (x + np)) {
    v = (x - np) / (x + np);
    s = (x - np) * v;/* s using v -- change by MM */
    ej = 2 * x * v;
    v = v * v;
    for (j = 1;; j++) { /* Taylor series */
      ej *= v;
      s1 = s + ej / ((j << 1) + 1);
      if (s1 == s) /* last term was effectively 0 */ {
        return (s1);
      }
      s = s1;
    }
  }
  /* else:  | x - np |  is not too small */
  return (x * Math.log(x / np) + np - x);
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double pnf(double x, double df1, double df2, double ncp, boolean lower_tail, boolean log_p) {
  double y;
  if (DoubleVector.isNaN(x) || DoubleVector.isNaN(df1) || DoubleVector.isNaN(df2) || DoubleVector.isNaN(ncp)) {
    return x + df2 + df1 + ncp;
  }
  if (df1 <= 0. || df2 <= 0. || ncp < 0) {
    return DoubleVector.NaN;
  }
  if (!DoubleVector.isFinite(ncp)) {
    return DoubleVector.NaN;
  }
  if (!DoubleVector.isFinite(df1) && !DoubleVector.isFinite(df2)) { /* both +Inf */
    return DoubleVector.NaN;
  }

  //R_P_bounds_01(x, 0., ML_POSINF);
  if (x <= 0.0) {
    return SignRank.R_DT_0(lower_tail, log_p);
  }
  if (x >= Double.POSITIVE_INFINITY) {
    return SignRank.R_DT_1(lower_tail, log_p);
  }

  if (df2 > 1e8) /* avoid problems with +Inf and loss of accuracy */ {
    return Distributions.pnchisq(x * df1, df1, ncp, lower_tail, log_p);
  }

  y = (df1 / df2) * x;
  return Beta.pnbeta2(y / (1. + y), 1. / (1. + y), df1 / 2., df2 / 2.,
          ncp, lower_tail, log_p);
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double qnf(double p, double df1, double df2, double ncp, boolean lower_tail, boolean log_p) {
  double y;

  if (DoubleVector.isNaN(p) || DoubleVector.isNaN(df1) || DoubleVector.isNaN(df2) || DoubleVector.isNaN(ncp)) {
    return p + df1 + df2 + ncp;
  }

  if (df1 <= 0. || df2 <= 0. || ncp < 0) {
    return DoubleVector.NaN;
  }
  if (!DoubleVector.isFinite(ncp)) {
    return DoubleVector.NaN;
  }
  if (!DoubleVector.isFinite(df1) && !DoubleVector.isFinite(df2)) {
    return DoubleVector.NaN;
  }

  //R_Q_P01_boundaries(p, 0, ML_POSINF);
  if ((log_p && p > 0) || (!log_p && (p < 0 || p > 1))) {
    return DoubleVector.NaN;
  }
  if (p == SignRank.R_DT_0(lower_tail, log_p)) {
    return 0;
  }
  if (p == SignRank.R_DT_1(lower_tail, log_p)) {
    return Double.POSITIVE_INFINITY;
  }
  //end of R_Q_P01_boundaries

  if (df2 > 1e8) /* avoid problems with +Inf and loss of accuracy */ {
    return Distributions.qnchisq(p, df1, ncp, lower_tail, log_p) / df1;
  }

  y = Beta.qnbeta(p, df1 / 2., df2 / 2., ncp, lower_tail, log_p);
  return y / (1 - y) * (df2 / df1);
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double pnchisq(double x, double df, double ncp, boolean lower_tail, boolean log_p) {
  double ans;

  if (DoubleVector.isNaN(x) || DoubleVector.isNaN(df) || DoubleVector.isNaN(ncp)) {
    return x + df + ncp;
  }

  if (!DoubleVector.isFinite(df) || !DoubleVector.isFinite(ncp)) {
    return DoubleVector.NaN;
  }

  if (df < 0. || ncp < 0.) {
    return DoubleVector.NaN;
  }

  ans = pnchisq_raw(x, df, ncp, 1e-12, 8 * SignRank.DBL_EPSILON, 1000000, lower_tail);
  if (ncp >= 80) {
    if (lower_tail) {
      ans = Math.min(ans, 1.0);  /* e.g., pchisq(555, 1.01, ncp = 80) */
    } else { /* !lower_tail */
      /* since we computed the other tail cancellation is likely */
      if (ans < 1e-10) {
        /*
         * Is there any function for ML_ERROR?????
         */
        //ML_ERROR(ME_PRECISION, "pnchisq");
      }
      ans = Math.max(ans, 0.0);  /* Precaution PR#7099 */
    }
  }
  if (!log_p) {
    return ans;
  }
  /* if ans is near one, we can do better using the other tail */
  if (ncp >= 80 || ans < 1 - 1e-8) {
    return Math.log(ans);
  }
  ans = pnchisq_raw(x, df, ncp, 1e-12, 8 * SignRank.DBL_EPSILON, 1000000, !lower_tail);
  return Math.log1p(-ans);
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double pgeom(double x, double p, boolean lower_tail, boolean log_p) {

    if (DoubleVector.isNaN(x) || DoubleVector.isNaN(p)) {
      return x + p;
    }

    if (p <= 0 || p > 1) {
      return DoubleVector.NaN;
    }

    if (x < 0.) {
      return SignRank.R_DT_0(lower_tail, log_p);
    }

    if (!DoubleVector.isFinite(x)) {
      return SignRank.R_DT_1(lower_tail, log_p);
    }

    x = Math.floor(x + 1e-7);

    if (p == 1.) { /* we cannot assume IEEE */
      x = lower_tail ? 1 : 0;
      return log_p ? Math.log(x) : x;
    }
    x = Math.log1p(-p) * (x + 1);
    if (log_p) {
      return SignRank.R_DT_Clog(x, lower_tail, log_p);
    } else {
      return lower_tail ? -Math.expm1(x) : Math.exp(x);
    }
  }
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double dgeom(double x, double p, boolean give_log) {
  double prob;


  if (DoubleVector.isNaN(x) || DoubleVector.isNaN(p)) {
    return x + p;
  }


  if (p <= 0 || p > 1) {
    return DoubleVector.NaN;
  }

  if (SignRank.R_D_nonint(x, true, give_log)) {
    return SignRank.R_D__0(true, give_log);
  }

  if (x < 0 || !DoubleVector.isFinite(x) || p == 0) {
    return SignRank.R_D__0(true, give_log);
  }
  x = SignRank.R_D_forceint(x);

  /* prob = (1-p)^x, stable for small p */
  prob = Binom.dbinom_raw(0., x, p, 1 - p, give_log);

  return ((give_log) ? Math.log(p) + prob : p * prob);
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double pnbinom_mu(double x, double size, double mu, boolean lower_tail, boolean log_p) {

    if (DoubleVector.isNaN(x) || DoubleVector.isNaN(size) || DoubleVector.isNaN(mu)) {
      return x + size + mu;
    }
    if (!DoubleVector.isFinite(size) || !DoubleVector.isFinite(mu)) {
      return DoubleVector.NaN;
    }

    if (size <= 0 || mu < 0) {
      return DoubleVector.NaN;
    }

    if (x < 0) {
      SignRank.R_DT_0(lower_tail, log_p);
    }
    if (!DoubleVector.isFinite(x)) {
      return SignRank.R_DT_1(lower_tail, log_p);
    }
    x = Math.floor(x + 1e-7);
    /* return
     * pbeta(pr, size, x + 1, lower_tail, log_p);  pr = size/(size + mu), 1-pr = mu/(size+mu)
     *
     *= pbeta_raw(pr, size, x + 1, lower_tail, log_p)
     *            x.  pin   qin
     *=  bratio (pin,  qin, x., 1-x., &w, &wc, &ierr, log_p),  and return w or wc ..
     *=  bratio (size, x+1, pr, 1-pr, &w, &wc, &ierr, log_p) */
    {
      int[] ierr = new int[1];
      double[] w = new double[1];
      double[] wc = new double[1];
      Utils.bratio(size, x + 1, size / (size + mu), mu / (size + mu), w, wc, ierr, log_p);
      return lower_tail ? w[0] : wc[0];
    }
  }
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double dnf(double x, double df1, double df2, double ncp, boolean give_log) {
  double y, z, f;

  if (DoubleVector.isNaN(x) || DoubleVector.isNaN(df1) || DoubleVector.isNaN(df2) || DoubleVector.isNaN(ncp)) {
    return x + df2 + df1 + ncp;
  }

  /* want to compare dnf(ncp=0) behavior with df() one, hence *NOT* :
   * if (ncp == 0)
   *   return df(x, df1, df2, give_log); */

  if (df1 <= 0. || df2 <= 0. || ncp < 0) {
    return DoubleVector.NaN;
  }
  if (x < 0.) {
    return (SignRank.R_D__0(true, give_log));
  }
  if (!DoubleVector.isFinite(ncp)) { /* ncp = +Inf -- FIXME?: in some cases, limit exists */
    return DoubleVector.NaN;
  }

  /* This is not correct for  df1 == 2, ncp > 0 - and seems unneeded:
   *  if (x == 0.) return(df1 > 2 ? R_D__0 : (df1 == 2 ? R_D__1 : ML_POSINF));
   */
  if (!DoubleVector.isFinite(df1) && !DoubleVector.isFinite(df2)) { /* both +Inf */
    /* PR: not sure about this (taken from  ncp==0)  -- FIXME ? */
    if (x == 1.) {
      return Double.POSITIVE_INFINITY;
    }
    /* else */ return SignRank.R_D__0(true, give_log);
  }
  if (!DoubleVector.isFinite(df2)) /* i.e.  = +Inf */ {
    return df1 * Distributions.dnchisq(x * df1, df1, ncp, give_log);
  }
  /*	 ==  dngamma(x, df1/2, 2./df1, ncp, give_log)  -- but that does not exist */
  if (df1 > 1e14 && ncp < 1e7) {
    /* includes df1 == +Inf: code below is inaccurate there */
    f = 1 + ncp / df1; /* assumes  ncp << df1 [ignores 2*ncp^(1/2)/df1*x term] */
    z = Distributions.dgamma(1. / x / f, df2 / 2, 2. / df2, give_log);
    return give_log ? z - 2 * Math.log(x) - Math.log(f) : z / (x * x) / f;
  }

  y = (df1 / df2) * x;
  z = Distributions.dnbeta(y / (1 + y), df1 / 2., df2 / 2., ncp, give_log);
  return give_log
          ? z + Math.log(df1) - Math.log(df2) - 2 * Math.log1p(y)
          : z * (df1 / df2) / (1 + y) / (1 + y);
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static void rmultinom(/*Session context,*/ int n, double[] prob, int K, int[] rN) /* `Return' vector  rN[1:K] {K := length(prob)}
 *  where rN[j] ~ Bin(n, prob[j]) ,  sum_j rN[j] == n,  sum_j prob[j] == 1,
 */ {
  int k;
  double pp;
  double p_tot = 0.; /* LDOUBLE */
  /* This calculation is sensitive to exact values, so we try to
  ensure that the calculations are as accurate as possible
  so different platforms are more likely to give the same
  result. */


  if (K < 1) {
    //ML_ERROR(ME_DOMAIN, "rmultinom"); 
    return;
  }
  if (n < 0) {
    return;
  }

  if (K == IntVector.NA || K < 1) {
    //ML_ERROR(ME_DOMAIN, "rmultinom"); 
    return;
  }

  if (n == IntVector.NA || n < 0) {
    return;
  }


  /* Note: prob[K] is only used here for checking  sum_k prob[k] = 1 ;
   *       Could make loop one shorter and drop that check !
   */
  for (k = 0; k < K; k++) {
    pp = prob[k];
    if (!DoubleVector.isFinite(pp) || pp < 0. || pp > 1.) {
      //ML_ERR_ret_NAN(k);
      return;
    }
    p_tot += pp;
    rN[k] = 0;
  }

  if (Math.abs(p_tot - 1.) > 1e-7) //MATHLIB_ERROR(_("rbinom: probability sum should be 1, but is %g"), (double) p_tot);
  {
    if (n == 0) {
      return;
    }
  }

  if (K == 1 && p_tot == 0.) {
    return; /* trivial border case: do as rbinom */
  }

  /* Generate the first K-1 obs. via binomials */

  for (k = 0; k < K - 1; k++) { /* (p_tot, n) are for "remaining binomial" */
    if (prob[k] != 0.0) {
      pp = prob[k] / p_tot;
      /* printf("[%d] %.17f\n", k+1, pp); */
      rN[k] = ((pp < 1.) ? (int) Binom.rbinom(/*context,*/ (double) n, pp)
              : /*>= 1; > 1 happens because of rounding */ n);
      n -= rN[k];
    } else {
      rN[k] = 0;
    }
    if (n <= 0) /* we have all*/ {
      return;
    }
    p_tot -= prob[k]; /* i.e. = sum(prob[(k+1):K]) */
  }
  rN[K - 1] = n;
  return;
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double qnt(double p, double df, double ncp, boolean lower_tail, boolean log_p) {
  final double accu = 1e-13;
  final double Eps = 1e-11; /* must be > accu */

  double ux, lx, nx, pp;


  if (DoubleVector.isNaN(p) || DoubleVector.isNaN(df) || DoubleVector.isNaN(ncp)) {
    return p + df + ncp;
  }

  if (!DoubleVector.isFinite(df)) {
    return DoubleVector.NaN;
  }

  /* Was
   * df = floor(df + 0.5);
   * if (df < 1 || ncp < 0) ML_ERR_return_NAN;
   */
  if (df <= 0.0) {
    return (DoubleVector.NaN);
  }

  if (ncp == 0.0 && df >= 1.0) {
    return Distributions.qt(p, df, lower_tail, log_p);
  }

  //R_Q_P01_boundaries(p, ML_NEGINF, ML_POSINF);
  //#define R_Q_P01_boundaries(p, _LEFT_, _RIGHT_)
  if (log_p) {
    if (p > 0) {
      return DoubleVector.NaN;
    }
    if (p == 0) /* upper bound*/ {
      return lower_tail ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
    }
    if (p == Double.NEGATIVE_INFINITY) {
      return lower_tail ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
    }
  } else { /* !log_p */
    if (p < 0 || p > 1) {
      return DoubleVector.NaN;
    }
    if (p == 0) {
      return lower_tail ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
    }
    if (p == 1) {
      return lower_tail ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
    }
  }

  p = Normal.R_DT_qIv(p, log_p ? 1.0 : 0.0, lower_tail ? 1.0 : 0.0);
  /* Invert pnt(.) :
   * 1. finding an upper and lower bound */
  if (p > 1 - SignRank.DBL_EPSILON) {
    return Double.POSITIVE_INFINITY;
  }
  pp = Math.min(1 - SignRank.DBL_EPSILON, p * (1 + Eps));
  for (ux = Math.max(1., ncp);
          ux < Double.MAX_VALUE && pnt(ux, df, ncp, true, false) < pp;
          ux *= 2);
  pp = p * (1 - Eps);
  for (lx = Math.min(-1., -ncp);
          lx > -Double.MAX_VALUE && pnt(lx, df, ncp, true, false) > pp;
          lx *= 2);

  /* 2. interval (lx,ux)  halving : */
  do {
    nx = 0.5 * (lx + ux);
    if (pnt(nx, df, ncp, true, false) > p) {
      ux = nx;
    } else {
      lx = nx;
    }
  } while ((ux - lx) / Math.abs(nx) > accu);

  return 0.5 * (lx + ux);
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double dnt(double x, double df, double ncp, boolean give_log) {
  double u;

  if (DoubleVector.isNaN(x) || DoubleVector.isNaN(df)) {
    return x + df;
  }

  /* If non-positive df then error */
  if (df <= 0.0) {
    return DoubleVector.NaN;
  }

  if (ncp == 0.0) {
    return Distributions.dt(x, df, give_log);
  }

  /* If x is infinite then return 0 */
  if (!DoubleVector.isFinite(x)) {
    return SignRank.R_D__0(true, give_log);
  }

  /* If infinite df then the density is identical to a
  normal distribution with mean = ncp.  However, the formula
  loses a lot of accuracy around df=1e9
   */
  if (!DoubleVector.isFinite(df) || df > 1e8) {
    return Distributions.dnorm(x, ncp, 1., give_log);
  }

  /* Do calculations on log scale to stabilize */

  /* Consider two cases: x ~= 0 or not */
  if (Math.abs(x) > Math.sqrt(df * SignRank.DBL_EPSILON)) {
    u = Math.log(df) - Math.log(Math.abs(x))
            + Math.log(Math.abs(Distributions.pnt(x * Math.sqrt((df + 2) / df), df + 2, ncp, true, false)
            - Distributions.pnt(x, df, ncp, true, false)));
    /* FIXME: the above still suffers from cancellation (but not horribly) */
  } else {  /* x ~= 0 : -> same value as for  x = 0 */
    u = org.apache.commons.math.special.Gamma.logGamma((df + 1) / 2) - org.apache.commons.math.special.Gamma.logGamma(df / 2)
            - .5 * (Math.log(Math.PI) + Math.log(df) + ncp * ncp);
  }

  return (give_log ? u : Math.exp(u));
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double dnbeta(double x, double a, double b, double ncp, boolean give_log) {
  final double eps = 1.e-15;

  int kMax;
  double k, ncp2, dx2, d, D;
  double sum, term, p_k, q; /* They were LDOUBLE */


  if (DoubleVector.isNaN(x) || DoubleVector.isNaN(a) || DoubleVector.isNaN(b) || DoubleVector.isNaN(ncp)) {
    return x + a + b + ncp;
  }

  if (ncp < 0 || a <= 0 || b <= 0) {
    return DoubleVector.NaN;
  }

  if (!DoubleVector.isFinite(a) || !DoubleVector.isFinite(b) || !DoubleVector.isFinite(ncp)) {
    return DoubleVector.NaN;
  }

  if (x < 0 || x > 1) {
    return (SignRank.R_D__0(true, give_log));
  }

  if (ncp == 0) {
    return Distributions.dbeta(x, a, b, give_log);
  }

  /* New algorithm, starting with *largest* term : */
  ncp2 = 0.5 * ncp;
  dx2 = ncp2 * x;
  d = (dx2 - a - 1) / 2;
  D = d * d + dx2 * (a + b) - a;
  if (D <= 0) {
    kMax = 0;
  } else {
    D = Math.ceil(d + Math.sqrt(D));
    kMax = (D > 0) ? (int) D : 0;
  }

  /* The starting "middle term" --- first look at it's log scale: */
  term = Distributions.dbeta(x, a + kMax, b, /* log = */ true);
  p_k = Poisson.dpois_raw(kMax, ncp2, true);
  if (x == 0. || !DoubleVector.isFinite(term) || !DoubleVector.isFinite(p_k)) /* if term = +Inf */ {
    return SignRank.R_D_exp(p_k + term, true, give_log);
  }

  /* Now if s_k := p_k * t_k  {here = exp(p_k + term)} would underflow,
   * we should rather scale everything and re-scale at the end:*/

  p_k += term; /* = log(p_k) + log(t_k) == log(s_k) -- used at end to rescale */
  /* mid = 1 = the rescaled value, instead of  mid = exp(p_k); */

  /* Now sum from the inside out */
  sum = term = 1. /* = mid term */;
  /* middle to the left */
  k = kMax;
  while (k > 0 && term > sum * eps) {
    k--;
    q = /* 1 / r_k = */ (k + 1) * (k + a) / (k + a + b) / dx2;
    term *= q;
    sum += term;
  }
  /* middle to the right */
  term = 1.;
  k = kMax;
  do {
    q = /* r_{old k} = */ dx2 * (k + a + b) / (k + a) / (k + 1);
    k++;
    term *= q;
    sum += term;
  } while (term > sum * eps);

  return SignRank.R_D_exp(p_k + Math.log(sum), true, give_log);
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double qnbeta(double p, double a, double b, double ncp,
        /*@Recycle(false)*/ boolean lower_tail, /*@Recycle(false)*/ boolean log_p) {
  final double accu = 1e-15;
  final double Eps = 1e-14; /* must be > accu */

  double ux, lx, nx, pp;


  if (DoubleVector.isNaN(p) || DoubleVector.isNaN(a) || DoubleVector.isNaN(b) || DoubleVector.isNaN(ncp)) {
    return p + a + b + ncp;
  }

  if (!DoubleVector.isFinite(a)) {
    return DoubleVector.NaN;
  }

  if (ncp < 0. || a <= 0. || b <= 0.) {
    return DoubleVector.NaN;
  }

  //R_Q_P01_boundaries(p, 0, 1);
  if ((log_p && p > 0) || (!log_p && (p < 0 || p > 1))) {
    return DoubleVector.NaN;
  }
  if (p == SignRank.R_DT_0(lower_tail, log_p)) {
    return 0.0;
  }
  if (p == SignRank.R_DT_1(lower_tail, log_p)) {
    return 1.0;
  }
  //end of R_Q_P01_boundaries

  p = Normal.R_DT_qIv(p, log_p ? 1.0 : 0.0, lower_tail ? 1.0 : 0.0);

  /* Invert pnbeta(.) :
   * 1. finding an upper and lower bound */
  if (p > 1 - SignRank.DBL_EPSILON) {
    return 1.0;
  }
  pp = Math.min(1 - SignRank.DBL_EPSILON, p * (1 + Eps));
  for (ux = 0.5;
          ux < 1 - SignRank.DBL_EPSILON && pnbeta(ux, a, b, ncp, true, false) < pp;
          ux = 0.5 * (1 + ux));
  pp = p * (1 - Eps);
  for (lx = 0.5;
          lx > Double.MIN_VALUE && pnbeta(lx, a, b, ncp, true, false) > pp;
          lx *= 0.5);

  /* 2. interval (lx,ux)  halving : */
  do {
    nx = 0.5 * (lx + ux);
    if (pnbeta(nx, a, b, ncp, true, false) > p) {
      ux = nx;
    } else {
      lx = nx;
    }
  } while ((ux - lx) / nx > accu);

  return 0.5 * (ux + lx);
}
 

import org.renjin.sexp.DoubleVector; //导入方法依赖的package包/类
public static double dnbinom_mu(double x, double size, double mu, boolean give_log) {
  /* originally, just set  prob :=  size / (size + mu)  and called dbinom_raw(),
   * but that suffers from cancellation when   mu << size  */
  double ans, p;


  if (DoubleVector.isNaN(x) || DoubleVector.isNaN(size) || DoubleVector.isNaN(mu)) {
    return x + size + mu;
  }

  if (mu < 0 || size < 0) {
    return DoubleVector.NaN;
  }

  //R_D_nonint_check(x);

  if (SignRank.R_D_nonint(x, true, give_log)) {
    //MATHLIB_WARNING("non-integer x = %f", x);	
    //How to warn??
    return SignRank.R_D__0(true, give_log);
  }

  if (x < 0 || !DoubleVector.isFinite(x)) {
    return SignRank.R_D__0(true, give_log);
  }
  x = SignRank.R_D_forceint(x);

  if (x == 0)/* be accerate, both for n << mu, and n >> mu :*/ {
    return SignRank.R_D_exp(size * (size < mu ? Math.log(size / (size + mu)) : Math.log1p(-mu / (size + mu))), true, give_log);
  }
  if (x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */
    /* FIXME --- 1e-8 shows problem; rather use algdiv() from ./toms708.c */
    return SignRank.R_D_exp(x * Math.log(size * mu / (size + mu)) - mu - org.apache.commons.math.special.Gamma.logGamma(x + 1)
            + Math.log1p(x * (x - 1) / (2 * size)), true, give_log);
  }
  /* else: no unnecessary cancellation inside dbinom_raw, when
   * x_ = size and n_ = x+size are so close that n_ - x_ loses accuracy
   */
  ans = dbinom_raw(size, x + size, size / (size + mu), mu / (size + mu), give_log);
  p = ((double) size) / (size + x);
  return ((give_log) ? Math.log(p) + ans : p * ans);
}