wnchyppr.cpp

很好用的随机数产生器

  gammal2 = sqrtpi * 0.75;
  accur = accuracy * sum;

  // summation loop
  for (k=3; k <= degree; k++) {
    phideri[k] = 0.;
    
    // loop for all colors
    for (i=0; i<colors; i++) {
      if (rho[i] == 0.) continue;
      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] += x[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 CMultiWalleniusNCHypergeometric::integrate(void) {
  // 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) {
    // 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 CMultiWalleniusNCHypergeometric::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, j;

  // 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 (j = 0; j < IPOINTS; j++) {
    tau = ab + delta * xval[j];
    ltau = log(tau);
	taur = r * ltau;
    y = 0.;
    for (i = 0; i < colors; i++) {
      // possible loss of precision due to subtraction here:
      if (omega[i]) {
        y += pow1pow(taur*omega[i],x[i]);}} // ln((1-e^taur*omegai)^xi)
    y += rdm1*ltau + bico;
	if (y > -50.) sum += weights[j] * exp(y);}
  return delta * sum;}


double CMultiWalleniusNCHypergeometric::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
  double t, t1;                        // independent variable
  double rho[MAXCOLORS];               // r*omega[i]
  double q;                            // t^rho[i] / (1-t^rho[i])
  double q1;                           // 1-t^rho[i]
  double zeta[MAXCOLORS][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
  t = 0.5 * (t_from + t_to);
  for (i=0; i<colors; i++) {           // calculate zeta coefficients
    rho[i] = r * omega[i];
    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)
      if (rho[i] == 0.) continue;
      q1 = pow2_1(rho[i]*log2t,&q);
      q /= q1;
      phi[1] -= x[i] * zeta[i][1][1] * q;
      phi[2] -= x[i] * q * (zeta[i][1][2] + q * zeta[i][2][2]);
      phi[3] -= x[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;
    if (++iter > 20) FatalError("Search for inflection point failed in function CMultiWalleniusNCHypergeometric::search_inflect");
    }
  while (fabs(t - t1) > 1E-5);
  return t;}


double CMultiWalleniusNCHypergeometric::probability(long int * x_) {
  // calculate probability function. choosing best method
  int i, j, em;
  long int xsum;
  x = x_;

  for (xsum=i=0; i<colors; i++)  xsum += x[i];
  if (xsum != n) 
    FatalError("sum of x values not equal to n in function CMultiWalleniusNCHypergeometric::probability");

  if (colors < 3) { 
    if (colors <= 0) return 1.;
    if (colors == 1) return x[0] == m[0];
    // colors = 2
    if (omega[1] == 0.) return x[0] == m[0];
    return CWalleniusNCHypergeometric(n,m[0],N,omega[0]/omega[1],accuracy).probability(x[0]);}

  for (i=j=em=0; i<colors; i++) {
    if (x[i] > m[i]) return 0.;
    if (omega[i]==0. && x[i]) return 0.;
    if (x[i] > 0) j++;
    if (x[i] == m[i]) em++;}

  if (n == 0 || n == N) return 1.;

  if (j == 1) 
    return binoexpand();

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

  return integrate();}


/***********************************************************************
         Methods for CMultiWalleniusNCHypergeometricMoments
***********************************************************************/

double CMultiWalleniusNCHypergeometricMoments::moments(double * mu, double * variance, long int * combinations) {
  // calculates mean and variance of multivariate Wallenius noncentral 
  // hypergeometric distribution by calculating all combinations of x-values.
  // Return value = sum of all probabilities. The deviation of this value 
  // from 1 is a measure of the accuracy.
  // Returns the mean to mean[0...colors-1]
  // Returns the variance to variance[0...colors-1]
  double sumf;               // sum of all f(x) values
  long int msum;             // temporary sum
  int i;                     // loop counter

  // get approximate mean
  mean(sx);
  // round mean to integers
  for (i=0; i < colors; i++) {
    xm[i] = (long int)(sx[i]+0.4999999);}

  // set up for recursive loops
  for (i=colors-1, msum=0; i >= 0; i--) {
    remaining[i] = msum;  msum += m[i];}
  for (i=0; i<colors; i++)  sx[i] = sxx[i] = 0.;
  sn = 0;

  // recursive loops to calculate sums  
  sumf = loop(n, 0);

  // calculate mean and variance
  for (i = 0; i < colors; i++) {
    mu[i] = sx[i]/sumf;
    variance[i] = sxx[i]/sumf - sx[i]*sx[i]/(sumf*sumf);}

  // return combinations and sum
  if (combinations) *combinations = sn;
  return sumf;}
    

double CMultiWalleniusNCHypergeometricMoments::loop(long int n, int c) {
  // recursive function to loop through all combinations of x-values.
  // used by moments()
  long int x, x0;                      // x of color c
  long int xmin, xmax;                 // min and max of x[c]
  double s1, s2, sum = 0.;             // sum of f(x) values
  int i;                               // loop counter

  if (c < colors-1) {
    // not the last color
    // calculate min and max of x[c] for given x[0]..x[c-1]
    xmin = n - remaining[c];  if (xmin < 0) xmin = 0;
    xmax = m[c];  if (xmax > n) xmax = n;
    x0 = xm[c];  if (x0 < xmin) x0 = xmin;  if (x0 > xmax) x0 = xmax;
    // loop for all x[c] from mean and up
    for (x = x0, s2 = 0.; x <= xmax; x++) {
      xi[c] = x;
      sum += s1 = loop(n-x, c+1); // recursive loop for remaining colors
      if (s1 < accuracy && s1 < s2) break; // stop when values become negligible
      s2 = s1;}
    // loop for all x[c] from mean and down
    for (x = x0-1; x >= xmin; x--) {
      xi[c] = x;
      sum += s1 = loop(n-x, c+1); // recursive loop for remaining colors
      if (s1 < accuracy && s1 < s2) break; // stop when values become negligible
      s2 = s1;}}
  else {
    // last color
    xi[c] = n;
    s1 = probability(xi);
    for (i=0; i < colors; i++) {
      sx[i]  += s1 * xi[i];
      sxx[i] += s1 * xi[i] * xi[i];}
    sn++;
    sum = s1;}
  return sum;}