wnchyppr.cpp

很好用的随机数产生器

  double o;
  if (x > n/2) { // invert
    x1 = n-x; m1 = N-m; m2 = m; o = 1./omega;}
  else {
    x1 = x; m1 = m; m2 = N-m; o = omega;}

  if (x1 == 0) {
    return exp(FallingFactorial(m2,n) - FallingFactorial(m2+o*m1,n));}
    
  if (x1 == 1) {
    double d, e, q, q0, q1;
    q = FallingFactorial(m2,n-1);
    e = o*m1+m2;
    q1 = q - FallingFactorial(e,n);
    e -= o;
    q0 = q - FallingFactorial(e,n);
    d = e - (n-1);
    return m1*d*(exp(q0) - exp(q1));}

  FatalError("x > 1 not supported by function CWalleniusNCHypergeometric::binoexpand");
  return 0;}


double CWalleniusNCHypergeometric::laplace() {
  // Laplace's method with narrow integration interval, 
  // using error function residues table, defined in erfres.cpp
  // NOTE: This function can only be used when the integrand peak is narrow.
  // findpars() must be called before this function.

  const int COLORS = 2;         // number of colors
  const int MAXDEG = 40;        // arraysize
  int degree;                   // max expansion degree
  double accur;                 // stop expansion when terms below this threshold
  double omegai[COLORS] = {omega, 1.}; // weights for each color
  double xi[COLORS] = {x, n-x}; // number of each color sampled
  double f0;                    // factor outside integral
  double rho[COLORS];           // r*omegai
  double qi;                    // 2^(-rho)
  double qi1;                   // 1-qi
  double qq[COLORS];            // qi / qi1
  double eta[COLORS+1][MAXDEG+1]; // eta coefficients
  double phideri[MAXDEG+1];     // derivatives of phi
  double PSIderi[MAXDEG+1];     // derivatives of PSI
  double epsilon;               // factor for integration interval
  double epsilonmax;            // max possible epsilon
  double * erfres;              // pointer to table of error function expansion residues
  int erfreslen;                // length of table

  // variables in asymptotic summation
  static const double sqrtpi = 1.772453850905515993; // sqrt(pi)
  static const double sqrt8  = 2.828427124746190098; // sqrt(8)
  double qqpow;                 // qq^j
  double pow2k;                 // 2^k
  double mfac;                  // (-1)^(k-1)*(k-1)!
  double kfac;                  // k!
  double bino;                  // binomial coefficient  
  double wr;                    // 1/w, w = integrand peak width
  double vr;                    // 1/v, v = integration interval
  double v2m2;                  // (2*v)^(-2)
  double v2mk1;                 // (2*v)^(-k-1)
  double gammal2;               // Gamma((k+1)/2);
  double s;                     // summation term
  double sum;                   // Taylor sum

  int i;                        // loop counter for color
  int j;                        // loop counter for derivative
  int k;                        // loop counter for expansion degree
  int ll;                       // k/2
  int converg=0;                // number of consequtive terms below accuracy


  // initialize
  for (k=0; k<=2; k++)  phideri[k] = PSIderi[k] = 0;

  // find rho[i], qq[i], first eta coefficients, and zero'th derivative of phi
  for (i=0; i<COLORS; i++) {
    rho[i] = r * omegai[i];
    if (rho[i] > 40.) {
      qi=0.;  qi1 = 1.;}  // avoid underflow
    else {
      qi1 = pow2_1(-rho[i], &qi);}   // qi=2^(-rho), qi1=1.-2^(-rho)
    qq[i] = qi / qi1;    // 2^(-r*omegai)/(1.-2^(-r*omegai))
    // peak = zero'th derivative
    phideri[0] += xi[i] * log1mx(qi, qi1);
    // eta coefficients
    eta[i][0] = 0.;
    eta[i][1] = eta[i][2] = rho[i]*rho[i];}

  // r, rd, and w must be calculated by findpars()
  // zero'th derivative
  phideri[0] -= (rd - 1.) * M_LN2;
  // scaled factor outside integral
  f0 = rd * exp(phideri[0] + lnbico());
  // calculate narrowed integration interval
  vr = sqrt8 * w;
  wr = 1. / w;
  phideri[2] = phi2d;
  epsilonmax = 0.5 * wr;
  if (accuracy < 1.E-10 && epsilonmax >= 8.) {
    epsilon = 8.; erfres = erfres80; erfreslen = sizeof(erfres80)/sizeof(erfres80[0]);}
  else if (accuracy < 1.E-8 && epsilonmax >= 7.) {
    epsilon = 7.; erfres = erfres70; erfreslen = sizeof(erfres70)/sizeof(erfres70[0]);}
  else if (accuracy < 1.E-6 && epsilonmax >= 6.) {
    epsilon = 6.; erfres = erfres60; erfreslen = sizeof(erfres60)/sizeof(erfres60[0]);}
  else {
    epsilon = 5.; erfres = erfres50; erfreslen = sizeof(erfres50)/sizeof(erfres50[0]);
    }
  if (epsilon > epsilonmax) FatalError("Laplace failed. peak width too high in function CWalleniusNCHypergeometric::laplace");
  degree = MAXDEG;
  if (degree >= erfreslen*2) degree = erfreslen*2-2;

  // variables used in summation
  v2m2 = 0.25*vr*vr;        // (2*v)^(-2)

  // set up for starting loop at k=3
  PSIderi[0] = 1.;
  pow2k = 8.;
  mfac = 2.;
  kfac = 2.;
  sum = 0.5 * vr * sqrtpi * erfres[0];
  v2mk1 = 0.5*vr*v2m2*v2m2;
  gammal2 = sqrtpi * 0.75;
  accur = accuracy * sum;

  // summation loop
  for (k=3; k <= degree; k++) {
    phideri[k] = 0.;
    
    // loop for all (2) colors
    for (i=0; i<COLORS; i++) {
      eta[i][k] = 0.;
      // backward loop for all powers
      for (j=k; j>0; j--) {
        // find coefficients recursively from previous coefficients
        eta[i][j]  =  eta[i][j]*(j*rho[i]-(k-2)) +  eta[i][j-1]*rho[i]*(j-1);}
      qqpow = 1.;
      // forward loop for all powers
      for (j=1; j<=k; j++) {
        qqpow *= qq[i];   // qq^j
        // contribution to derivative
        phideri[k] += xi[i] * eta[i][j] * qqpow;}}

    // finish calculation of derivatives
    phideri[k] = -pow2k * phideri[k] + 2*(1-k)*phideri[k-1];

    pow2k *= 2.;    // 2^k
    mfac *= -k;     // (-1)^(k-1)*(k-1)!

    // loop to calculate derivatives of PSI from derivatives of psi.
    // terms # 0, 1, 2, k-2, and k-1 are zero and not included in loop.
    // The j'th derivatives of psi are identical to the derivatives of phi for j>2, and
    // zero for j=1,2. Hence we are using phideri[j] for j>2 here.
    PSIderi[k] = phideri[k]; // this is term # k
    bino = 0.5*(k-1)*(k-2); // binomial coefficient for term # 3
    for (j=3; j < k-2; j++) { // loop for remaining nonzero terms (if k>5)
      PSIderi[k] += PSIderi[k-j] * phideri[j] * bino;
      bino *= double(k-j)/double(j);}

    kfac *= k;

    if ((k & 1) == 0) { // only for even k
      ll = k/2;
      s = (PSIderi[k]/kfac)*v2mk1* gammal2 * erfres[ll];
      sum += s;

      // check for convergence of Taylor expansion
      if (fabs(s) < accur) converg++; else converg = 0;
      if (converg > 1) break;

      // update recursive expressions
      v2mk1 *= v2m2;
      gammal2 *= k/2 + 0.5;}}

  // multiply by terms outside integral  
  return f0 * sum;}


double CWalleniusNCHypergeometric::integrate() {
  // Wallenius non-central hypergeometric distribution function
  // calculation by numerical integration with variable-length steps
  // NOTE: findpars() must be called before this function.
  double s;                            // result of integration step
  double sum;                          // integral
  double ta, tb;                       // subinterval for integration step

  lnbico();                            // compute log of binomial coefficients

  // choose method:
  if (w < 0.02 || w < 0.1 && (x==m || n-x==N-m) && accuracy > 1E-6) {
    // normal method. Step length determined by peak width w
    double delta, s1;
    s1 = accuracy < 1E-9 ? 0.5 : 1.;
    delta = s1 * w;                    // integration steplength
    ta = 0.5 + 0.5 * delta;
    sum = integrate_step(1.-ta, ta);   // first integration step around center peak
    do {
      tb = ta + delta;
      if (tb > 1.) tb = 1.;
      s  = integrate_step(ta, tb);     // integration step to the right of peak
      s += integrate_step(1.-tb,1.-ta);// integration step to the left of peak
      sum += s;
      if (s < accuracy * sum) break;   // stop before interval finished if accuracy reached
      ta = tb;
      if (tb > 0.5 + w) delta *= 2.;}  // increase step length far from peak
    while (tb < 1.);}

  else {
    // difficult situation. Step length determined by inflection points
    double t1, t2, tinf, delta, delta1;
    sum = 0.;
    // do left and right half of integration interval separately:
    for (t1=0., t2=0.5; t1 < 1.; t1+=0.5, t2+=0.5) { 
      // integrate from 0 to 0.5 or from 0.5 to 1
      tinf = search_inflect(t1, t2);   // find inflection point
      delta = tinf - t1; if (delta > t2 - tinf) delta = t2 - tinf; // distance to nearest endpoint
      delta *= 1./7.;                  // 1/7 will give 3 steps to nearest endpoint
      if (delta < 1E-4) delta = 1E-4;
      delta1 = delta;
      // integrate from tinf forwards to t2
      ta = tinf;
      do {
        tb = ta + delta1;
        if (tb > t2 - 0.25*delta1) tb = t2; // last step of this subinterval
        s = integrate_step(ta, tb);         // integration step
        sum += s;
        delta1 *= 2;                        // double steplength
        if (s < sum * 1E-4) delta1 *= 8.;   // large step when s small
        ta = tb;}
      while (tb < t2);
      if (tinf) {
        // integrate from tinf backwards to t1
        tb = tinf;
        do {
          ta = tb - delta;
          if (ta < t1 + 0.25*delta) ta = t1; // last step of this subinterval
          s = integrate_step(ta, tb);        // integration step
          sum += s;
          delta *= 2;                        // double steplength
          if (s < sum * 1E-4) delta *= 8.;   // large step when s small
          tb = ta;}
        while (ta > t1);}}}

  return sum * rd;}


double CWalleniusNCHypergeometric::integrate_step(double ta, double tb) {
  // integration subprocedure used by integrate()
  // makes one integration step from ta to tb using Gauss-Legendre method.
  // result is scaled by multiplication with exp(bico)
  double ab, delta, tau, ltau, y, sum, taur, rdm1;
  int i;

  // define constants for Gauss-Legendre integration with IPOINTS points
  #define IPOINTS  8  // number of points in each integration step

  #if IPOINTS == 3
  static const double xval[3]    = {-.774596669241,0,0.774596668241};
  static const double weights[3] = {.5555555555555555,.88888888888888888,.55555555555555};
  #elif IPOINTS == 4
  static const double xval[4]    = {-0.861136311594,-0.339981043585,0.339981043585,0.861136311594},
  static const double weights[4] = {0.347854845137,0.652145154863,0.652145154863,0.347854845137};
  #elif IPOINTS == 5
  static const double xval[5]    = {-0.906179845939,-0.538469310106,0,0.538469310106,0.906179845939};
  static const double weights[5] = {0.236926885056,0.478628670499,0.568888888889,0.478628670499,0.236926885056};
  #elif IPOINTS == 6
  static const double xval[6]    = {-0.932469514203,-0.661209386466,-0.238619186083,0.238619186083,0.661209386466,0.932469514203};
  static const double weights[6] = {0.171324492379,0.360761573048,0.467913934573,0.467913934573,0.360761573048,0.171324492379};
  #elif IPOINTS == 8
  static const double xval[8]    = {-0.960289856498,-0.796666477414,-0.525532409916,-0.183434642496,0.183434642496,0.525532409916,0.796666477414,0.960289856498};
  static const double weights[8] = {0.10122853629,0.222381034453,0.313706645878,0.362683783378,0.362683783378,0.313706645878,0.222381034453,0.10122853629};
  #elif IPOINTS == 12
  static const double xval[12]   = {-0.981560634247,-0.90411725637,-0.769902674194,-0.587317954287,-0.367831498998,-0.125233408511,0.125233408511,0.367831498998,0.587317954287,0.769902674194,0.90411725637,0.981560634247};
  static const double weights[12]= {0.0471753363866,0.106939325995,0.160078328543,0.203167426723,0.233492536538,0.249147045813,0.249147045813,0.233492536538,0.203167426723,0.160078328543,0.106939325995,0.0471753363866};
  #elif IPOINTS == 16
  static const double xval[16]   = {-0.989400934992,-0.944575023073,-0.865631202388,-0.755404408355,-0.617876244403,-0.458016777657,-0.281603550779,-0.0950125098376,0.0950125098376,0.281603550779,0.458016777657,0.617876244403,0.755404408355,0.865631202388,0.944575023073,0.989400934992};
  static const double weights[16]= {0.027152459411,0.0622535239372,0.0951585116838,0.124628971256,0.149595988817,0.169156519395,0.182603415045,0.189450610455,0.189450610455,0.182603415045,0.169156519395,0.149595988817,0.124628971256,0.0951585116838,0.0622535239372,0.027152459411};
  #else
  #error // IPOINTS must be a value for which the tables are defined
  #endif

  delta = 0.5 * (tb - ta);
  ab = 0.5 * (ta + tb);
  rdm1 = rd - 1.;
  sum = 0;

  for (i = 0; i < IPOINTS; i++) {
    tau = ab + delta * xval[i];
    ltau = log(tau);
    taur = r * ltau;
    // possible loss of precision due to subtraction here:
    y = pow1pow(taur*omega,x) + pow1pow(taur,n-x) + rdm1*ltau + bico;
      if (y > -50.) sum += weights[i] * exp(y);}
  return delta * sum;}


double CWalleniusNCHypergeometric::search_inflect(double t_from, double t_to) {
  // search for an inflection point of the integrand PHI(t) in the interval
  // t_from < t < t_to
  const int COLORS = 2;                // number of colors
  double t, t1;                        // independent variable
  double rho[COLORS];                  // r*omega[i]
  double q;                            // t^rho[i] / (1-t^rho[i])
  double q1;                           // 1-t^rho[i]
  double xx[COLORS];                   // x[i]
  double zeta[COLORS][4][4];           // zeta[i,j,k] coefficients
  double phi[4];                       // derivatives of phi(t) = log PHI(t)
  double Z2;                           // PHI''(t)/PHI(t)
  double Zd;                           // derivative in Newton Raphson iteration
  double rdm1;                         // r * d - 1
  double tr;                           // 1/t
  double log2t;                        // log2(t)
  double method;                       // 0 for z2'(t) method, 1 for z3(t) method
  int i;                               // color
  int iter;                            // count iterations

  rdm1 = rd - 1.;
  if (t_from == 0 && rdm1 <= 1.) return 0.; //no inflection point
  rho[0] = r*omega;  rho[1] = r;
  xx[0] = x;  xx[1] = n - x;
  t = 0.5 * (t_from + t_to);
  for (i=0; i<COLORS; i++) {           // calculate zeta coefficients
    zeta[i][1][1] = rho[i];
    zeta[i][1][2] = rho[i] * (rho[i] - 1.);
    zeta[i][2][2] = rho[i] * rho[i];
    zeta[i][1][3] = zeta[i][1][2] * (rho[i] - 2.);
    zeta[i][2][3] = zeta[i][1][2] * rho[i] * 3.;
    zeta[i][3][3] = zeta[i][2][2] * rho[i] * 2.;}
  iter = 0;

  do {
    t1 = t;
    tr = 1. / t;
    log2t = log(t)*(1./M_LN2);
    phi[1] = phi[2] = phi[3] = 0.;
    for (i=0; i<COLORS; i++) {         // calculate first 3 derivatives of phi(t)
      q1 = pow2_1(rho[i]*log2t,&q);
      q /= q1;
      phi[1] -= xx[i] * zeta[i][1][1] * q;
      phi[2] -= xx[i] * q * (zeta[i][1][2] + q * zeta[i][2][2]);
      phi[3] -= xx[i] * q * (zeta[i][1][3] + q * (zeta[i][2][3] + q * zeta[i][3][3]));}
    phi[1] += rdm1;
    phi[2] -= rdm1;
    phi[3] += 2. * rdm1;
    phi[1] *= tr;
    phi[2] *= tr * tr;
    phi[3] *= tr * tr * tr;
    method = (iter & 2) >> 1;          // alternate between the two methods
    Z2 = phi[1]*phi[1] + phi[2];
    Zd = method*phi[1]*phi[1]*phi[1] + (2.+method)*phi[1]*phi[2] + phi[3];

    if (t < 0.5) {
      if (Z2 > 0) {
        t_from = t;}
      else {
        t_to = t;}
      if (Zd >= 0) { 
        // use binary search if Newton-Raphson iteration makes problems
        t = (t_from ? 0.5 : 0.2) * (t_from + t_to);}
      else {
        // Newton-Raphson iteration
        t -= Z2 / Zd;}}
    else {
      if (Z2 < 0) {
        t_from = t;}
      else {
        t_to = t;}
      if (Zd <= 0) {
        // use binary search if Newton-Raphson iteration makes problems
        t = 0.5 * (t_from + t_to);}
      else {
        // Newton-Raphson iteration
        t -= Z2 / Zd;}}
    if (t >= t_to) t = (t1 + t_to) * 0.5;
    if (t <= t_from) t = (t1 + t_from) * 0.5;