wnchyppr.cpp

很好用的随机数产生器

    if (++iter > 20) FatalError("Search for inflection point failed in function CWalleniusNCHypergeometric::search_inflect");
    }
  while (fabs(t - t1) > 1E-5);
  return t;}


double CWalleniusNCHypergeometric::probability(long int x_) {
  // calculate probability function. choosing best method
  x = x_;
  if (x < xmin || x > xmax) return 0.;
  if (xmin == xmax) return 1.;
  if (omega == 1.) { // hypergeometric
    return exp(lnbico() + LnFac(n) + LnFac(N-n) - LnFac(N));}

  long int x2 = n - x;
  long int x0 = x < x2 ? x : x2;
  int em = (x == m || x2 == N-m);

  if (x0 == 0 && n > 500) {
    return binoexpand();}

  if (double(n)*x0 < 1000 || double(n)*x0 < 10000 && (N > 1000.*n || em)) {
    return recursive();}

  if (x0 <= 1 && N-n < 3) { 
    return binoexpand();}

  findpars();
  if (w < 0.04 && E < 10 && (!em || w > 0.004)) {
    return laplace();}

  return integrate();}


int CWalleniusNCHypergeometric::MakeTable(double * table, long int MaxLength, long int * start, long int * end) {
  // Makes a table of Wallenius noncentral hypergeometric probability distribution.
  // Returns true if table is long enough. Otherwise it fills the table
  // with as many correct values as possible and returns false.
  // The tails are cut off where the values are < 0.01*accuracy, so that 
  // *start may be > xmin and *end may be < xmax.
  // The first and last x value represented in the table are returned in 
  // *start and *end. The resulting probability values are returned in the 
  // first (*end - *start + 1) positions of table. Any unused part of table 
  // is overwritten with garbage. It is the responsibility of the calling 
  // program to allocate sufficient space for the table. 
  
  const int BUFSIZE = 2048;            // length of internal table
  double accuracya;                    // absolute accuracy
  double p[BUFSIZE+2];                 // probabilities
  double * p1, * p2;                   // offset into p
  double mxo;                          // (m-x)*omega
  double Nmnx;                         // N-m-nu+x
  double y, y1;                        // probability. Save old p[x] before it is overwritten
  double d1, d2, dcom;                 // divisors in probability formula
  long int xi, nu;                     // xi, nu = recursion values of x, n
  long int x1, x2;                     // lowest and highest x or xi
  long int i1, i2;                     // index into table
  double area;                         // estimate of area needed for recursion method

  area = double(n)*(xmax - xmin + 1);
  if (area < 4000 || area < 10000 && N > 1000.*n) {
    // use recursion method
    accuracya = 0.005 * accuracy;      // absolute accuracy
    p1 = p2 = p + 1;                   // make space for p1[-1]
    p1[-1] = 0.;  p1[0]  = 1.;         // initialize for recursion
    x1 = x2 = 0;
    for (nu=1; nu<=n; nu++) {
      if (n - nu < xmin - x1 || p1[x1] < accuracya) {
        x1++;                          // increase lower limit when breakpoint passed or probability negligible
        p2--;}                         // compensate buffer offset in order to reduce storage space
      if (x2 < xmax && p1[x2] >= accuracya) {
        x2++;  y1 = 0.;}               // increase upper limit until x has been reached
      else {
        y1 = p1[x2];}
      if (p2 + x2 - p > BUFSIZE || x1 > x2) {
        goto ONE_BY_ONE;}              // Error: table length exceeded. Use other method

      mxo = (m-x2)*omega;
      Nmnx = N-m-nu+x2+1;
      for (xi = x2; xi >= x1; xi--) {  // backwards loop
        d2 = mxo + Nmnx;
        mxo += omega; Nmnx--;
        d1 = mxo + Nmnx;
        dcom = 1. / (d1 * d2);         // save a division by making common divisor
        y  = p1[xi-1]*mxo*d2*dcom + y1*(Nmnx+1)*d1*dcom;
        y1 = p1[xi-1];                 // (warning: pointer alias, can't swap instruction order)
        p2[xi] = y;}
      p1 = p2;}

    // return results
    i1 = i2 = x2 - x1 + 1;             // desired table length
    if (i2 > MaxLength) i2 = MaxLength;// limit table length
    *start = x1;  *end = x1 + i2 - 1;
    if (i2>0) memcpy(table, p+1, i2*sizeof(table[0]));// copy into external table
    return i1==i2;}                    // true if table size not reduced

  else {
    // area too big for recursion method
    // Calculate values one by one
    ONE_BY_ONE:
    accuracya = 0.01 * accuracy;       // absolute accuracy

    // start to fill table from the end and down. start with x = floor(mean)
    x2 = (long int)mean();
    x1 = x2 + 1;  i1 = MaxLength;
    while (x1 > xmin) {                // loop for left tail
      x1--;  i1--;
      y = probability(x1);
      table[i1] = y;
      if (y < accuracya) break;
      if (i1 == 0) break;}
    *start = x1;
    i2 = x2-x1+1; 
    if (i1>0 && i2>0) { // move numbers down to beginning of table
      memcpy(table, table+i1, i2*sizeof(table[0]));}
    // fill rest of table from mean and up
    i2--;
    while (x2 < xmax) {                // loop for right tail
      if (i2 == MaxLength-1) {
        *end = x2; return 0;}          // table full
      x2++;  i2++;
      y = probability(x2);
      table[i2] = y;
      if (y < accuracya) break;}
    *end = x2;
    return 1;}}


/***********************************************************************
    calculation methods in class CMultiWalleniusNCHypergeometric
***********************************************************************/
  
CMultiWalleniusNCHypergeometric::CMultiWalleniusNCHypergeometric(long int n_, long int * m_, double * odds_, int colors_, float accuracy_) {
  // constructor
  long int N1;
  int i;
  n = n_;  m = m_;  omega = odds_;  colors = colors_;  accuracy = accuracy_;
  r = 1.;
  for (N=N1=0, i=0; i < colors; i++) {
    if (m[i] < 0 || omega[i] < 0) FatalError("Parameter negative in constructor for CMultiWalleniusNCHypergeometric");
    N += m[i];
    if (omega[i]) N1 += m[i];}
  if (N < n) FatalError("Not enough items in constructor for CMultiWalleniusNCHypergeometric");
  if (N1< n) FatalError("Not enough items with nonzero weight in constructor for CMultiWalleniusNCHypergeometric");
  }


void CMultiWalleniusNCHypergeometric::mean(double * mu) {
  // calculate approximate mean of multivariate Wallenius noncentral hypergeometric 
  // distribution. Result is returned in mu[0..colors-1]
  double omeg[MAXCOLORS];              // scaled weights
  double omr;                          // reciprocal mean weight
  double t, t1;                        // independent variable in iteration
  double To, To1;                      // exp(t*omega[i]), 1-exp(t*omega[i])
  double H;                            // function to find root of
  double HD;                           // derivative of H
  int i;                               // color index
  int iter;                            // number of iterations

  // calculate mean weight
  for (omr=0., i=0; i < colors; i++) omr += omega[i]*m[i];
  omr = N / omr;
  // scale weights to make mean = 1
  for (i = 0; i < colors; i++) omeg[i] = omega[i] * omr;
  // Newton Raphson iteration
  iter = 0;  t = -1.;                  // first guess
  do {
    t1 = t;
    H = HD = 0.;
    // calculate H and HD
    for (i = 0; i < colors; i++) {
      if (omeg[i] != 0.) {
        To1 = pow2_1(t * (1./M_LN2) * omeg[i], &To);
        H += m[i] * To1;
        HD -= m[i] * omeg[i] * To;}}
    t -= (H-n) / HD;
    if (t >= 0) t = 0.5 * t1;
    if (++iter > 20) {
      FatalError("Search for mean failed in function CMultiWalleniusNCHypergeometric::mean");}}
  while (fabs(H - n) > 1E-3);
  // finished iteration. Get all mu[i]
  for (i = 0; i < colors; i++) {
    if (omeg[i] != 0.) {
      To1 = pow2_1(t * (1./M_LN2) * omeg[i]);
      mu[i] = m[i] * To1;}
    else {
      mu[i] = 0.;}}}


// implementations of different calculation methods
double CMultiWalleniusNCHypergeometric::binoexpand(void) {
  // binomial expansion of integrand
  // only implemented for x[i] = 0 for all but one i
  int i, j, k;
  double W = 0.;                       // total weight
  for (i=j=k=0; i<colors; i++) {
    W += omega[i] * m[i];
    if (x[i]) {
      j=i; k++;}}                      // find the nonzero x[i]
  if (k > 1) FatalError("More than one x[i] nonzero in CMultiWalleniusNCHypergeometric::binoexpand");
  return exp(FallingFactorial(m[j],n) - FallingFactorial(W/omega[j],n));}


double CMultiWalleniusNCHypergeometric::lnbico(void) {
  // natural log of binomial coefficients
  bico = 0.;
  int i;
  for (i=0; i<colors; i++) {
    if (x[i] < m[i] && omega[i]) {
      bico += LnFac(m[i]) - LnFac(x[i]) - LnFac(m[i]-x[i]);}}
  return bico;}


void CMultiWalleniusNCHypergeometric::findpars(void) {
  // calculate r, w, E
  // calculate d, E, r, w

  // find r to center peak of integrand at 0.5
  double dd;                           // scaled d
  double dr;                           // 1/d
  
  double z, zd, rr, lastr, rrc, rt, r2, r21, a, b, ro, k1;
  double omax;                         // highest omega
  double omaxr;                        // 1/omax
  double omeg[MAXCOLORS];              // scaled weights
  int i, j = 0;

  // find highest omega
  for (omax=0., i=0; i < colors; i++) {
    if (omega[i] > omax) omax = omega[i];}
  omaxr = 1. / omax;
  dd = E = 0.;
  for (i = 0; i < colors; i++) {
    // scale weights to make max = 1
    omeg[i] = omega[i] * omaxr;
    // calculate d and E
    dd += omeg[i] * (m[i]-x[i]);
    E  += omeg[i] * m[i];}
  dr = 1. / dd;
  E *= dr;
  rr = r * omax;
  if (rr <= dr) rr = 1.2 * dr;  // initial guess
  // Newton-Raphson iteration to find r
  do {
    lastr = rr;
	rrc = 1. / rr;
	z = dd - rrc;                      // z(r)
	zd = rrc * rrc;                    // z'(r)
    for (i=0; i<colors; i++) {
      rt = rr * omeg[i];
      if (rt < 100. && rt > 0.) {      // avoid overflow and division by 0
        r21 = pow2_1(rt, &r2);         // r2=2^r, r21=1.-2^r
	    a = omeg[i] / r21;             // omegai/(1.-2^r)
	    b = x[i] * a;                  // x*omegai/(1.-2^r)
	    z  += b;
	    zd += b * a * r2 * M_LN2;}}
    if (zd == 0) FatalError("can't find r in function CMultiWalleniusNCHypergeometric::findpars");
	rr -= z / zd;                      // next r
	if (rr <= dr) rr = lastr * 0.125 + dr * 0.875;
    if (++j == 70) FatalError("convergence problem searching for r in function CMultiWalleniusNCHypergeometric::findpars");
    }
  while (fabs(rr-lastr) > rr * 1.E-5);
  rd = rr * dd;
  r = rr * omaxr;

  // find peak width
  phi2d = 0.;
  for (i=0; i<colors; i++) {
    ro = rr * omeg[i];
    if (ro < 300 && ro > 0.) {         // avoid overflow and division by 0
      k1 = pow2_1(ro);
      k1 = -1. / k1;
      k1 = omeg[i] * omeg[i] * (k1 + k1*k1);}
    else k1 = 0.;
    phi2d += x[i] * k1;}
  phi2d *= -4. * rr * rr;
  if (phi2d > 0.) FatalError("peak width undefined in function CMultiWalleniusNCHypergeometric::findpars");
  w = 1./sqrt(-phi2d);}


double CMultiWalleniusNCHypergeometric::laplace(void) {
  // Laplace's method with narrow integration interval, 
  // using error function residue 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 MAXDEG = 40;        // arraysize
  int degree;                   // max expansion degree
  double accur;                 // stop expansion when terms below this threshold
  double f0;                    // factor outside integral
  double rho[MAXCOLORS];        // r*omegai
  double qi;                    // 2^(-rho)
  double qi1;                   // 1-qi
  double qq[MAXCOLORS];         // qi / qi1
  double eta[MAXCOLORS+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 * omega[i];
    if (rho[i] == 0.) continue;
    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] += x[i] * log1mx(qi, qi1);
    // eta coefficients
    eta[i][0] = 0.;
    eta[i][1] = eta[i][2] = rho[i]*rho[i];}

  // d, r, 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 CMultiWalleniusNCHypergeometric::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;