wnchyppr.cpp

很好用的随机数产生器

/************************** WNCHYPPR.CPP ********************** 2002-10-20 AF *
*
* Calculation of univariate and multivariate Wallenius noncentral 
* hypergeometric probability distribution.
*
* This file contains source code for the class CWalleniusNCHypergeometric 
* and CMultiWalleniusNCHypergeometricMoments defined in stocc.h.
*
* Documentation:
* ==============
* The file stocc.h contains class definitions.
* The file stocc.htm contains further instructions.
* The file nchyp.pdf, available from www.agner.org/random/theory 
* describes the theory of the calculation methods.
*
*******************************************************************************/

#include <math.h>      // math functions
#include <string.h>    // memcpy function
#include "stocc.h"     // class definition
#include "erfres.cpp"  // table of error function residues for Laplace's method

/***********************************************************************
         Shared functions
***********************************************************************/

double pow2_1(double q, double * y0 = 0) {
// calculate 2^q and (1-2^q) without loss of precision.
// return value is (1-2^q). 2^q is returned in *y0
  double y, y1, y2, qn, i, ifac;
  q *= M_LN2;
  if (fabs(q) > 0.1) {
    y = exp(q);
    y1 = 1. - y;}
  else { // expand 1-e^q = -summa(q^n/n!) to avoid loss of precision
    y1 = 0;  qn = i = ifac = 1;
    do {
      y2 = y1;
      qn *= q;  ifac *= i++;
      y1 -= qn / ifac;}
    while (y1 != y2);
    y = 1.-y1;}
  if (y0) *y0 = y;
  return y1;}
  

double pow1pow(double q, double x) {
// calculate ln((1-e^q)^x) without loss of precision
// Uses various Taylor expansions to avoid loss of precision
  double y, y1, y2, z1, z2, qn, i, ifac;

  if (fabs(q) > 0.1) {
    y = exp(q);  y1 = 1. - y;}
  else { // expand 1-e^q = -summa(q^n/n!) to avoid loss of precision
    y1 = 0;  qn = i = ifac = 1;
    do {
      y2 = y1;
      qn *= q;  ifac *= i++;
      y1 -= qn / ifac;}
    while (y1 != y2);
    y = 1. - y1;}
  if (y > 0.1) { // (1-y)^x calculated without problem
    return x * log(y1);}
  else { // expand ln(1-y) = -summa(y^n/n)
    y1 = i = 1.;  z1 = 0;
    do {
      z2 = z1;
      y1 *= y;
      z1 -= y1 / i++;}
    while (z1 != z2);
    return x * z1;}}
    

double log1mx(double x, double x1) {
// calculate log(1-x) without loss of precision when x is small
// parameter x1 must be = 1-x
  if (x > 0.03) {
    return log(x1);}
  else { // expand ln(1-x) = -summa(x^n/n)
    double y, z1, z2, i;
    y = i = 1.;  z1 = 0;
    do {
      z2 = z1;
      y *= x;
      z1 -= y / i++;}
    while (z1 != z2);
    return z1;}}
    

double LnFacr(double x) {
  // log factorial of non-integer x
  long int ix = (long int)(x);
  if (x == ix) return LnFac(ix); // x is integer
  double r, r2, D = 1., f;
  static const double             
    C0 =  0.918938533204672722,   // ln(sqrt(2*pi))
    C1 =  1./12.,
    C3 = -1./360.,
    C5 =  1./1260.,
    C7 = -1./1680.;
  if (x < 6.) {
    if (x == 0 || x == 1) return 0;
    while (x < 6) D *= ++x;}
  r  = 1. / x;  r2 = r*r;
  f = (x + 0.5)*log(x) - x + C0 + r*(C1 + r2*(C3 + r2*(C5 + r2*C7)));
  if (D != 1.) f -= log(D);
  return f;}
  

double FallingFactorial(double a, double b) {
  // calculates ln(a*(a-1)*(a-2)* ... * (a-b+1))

  if (b < 30 && int(b) == b && a < 1E10) {
    // direct calculation
    double f = 1.;
    for (int i=0; i<b; i++) f *= a--;
    return log(f);}
    
  if (a > 100.*b && b > 1.) {
    // combine Stirling formulas for a and (a-b) to avoid loss of precision
    double ar = 1./a;
    double cr = 1./(a-b);
    // calculate -log(1-b/a) by Taylor expansion
    double s = 0., lasts, n = 1., ba = b*ar, f = ba;
    do {
      lasts = s;
      s += f/n;
      f *= ba;
      n++;}
    while (s != lasts);
    return (a+0.5)*s + b*log(a-b) - b + (1./12.)*(ar-cr)    
    //- (1./360.)*(ar*ar*ar-cr*cr*cr)
    ;}

  // use LnFacr function
  return LnFacr(a)-LnFacr(a-b);}


/***********************************************************************
             Methods for class CWalleniusNCHypergeometric
***********************************************************************/

CWalleniusNCHypergeometric::CWalleniusNCHypergeometric(long int n_, long int m_, long int N_, double odds, double accuracy_) {
  // constructor
  n = n_; m = m_; N = N_; omega = odds; accuracy = accuracy_; // set parameters
  xmin = m + n - N;  if (xmin < 0) xmin = 0;  // calculate xmin and xmax
  xmax = n;  if (xmax > m) xmax = m;
  ParametersChanged = 1;  r = 1.;}            // initialize


double CWalleniusNCHypergeometric::mean(void) {
  // find approximate mean
  int iter;                            // number of iterations
  double a, b;                         // temporaries in calculation of first guess
  double mean, mean1;                  // iteration value of mean
  double m1r, m2r;                     // 1/m, 1/m2
  double e1, e2;                       // temporaries
  double g;                            // function to find root of
  double gd;                           // derivative of g
  double omegar;                       // 1/omega

  if (omega == 1.) { // simple hypergeometric
    return double(m)*n/N;}

  // calculate Cornfield mean of Fisher noncentral hypergeometric distribution as first guess
  a = (m+n)*omega + (N-m-n); 
  b = a*a - 4.*omega*(omega-1.)*m*n;
  b = b > 0. ? sqrt(b) : 0.;
  mean = (a-b)/(2.*(omega-1.));
  if (mean < xmin) mean = xmin;
  if (mean > xmax) mean = xmax;

  iter = 0;
  m1r = 1./m;  m2r = 1./(N-m);  omegar = 1./omega; 

  if (omega > 1.) {
    do { // Newton Raphson iteration
      mean1 = mean;
      e1 = 1.-(n-mean)*m2r;
      if (e1 < 1E-14) {
        e2 = 0.;} // avoid underflow
      else {
        e2 = pow(e1,omega-1.);}
      g = e2*e1 + (mean-m)*m1r;
      gd = e2*omega*m2r + m1r;
      mean -= g / gd;
      if (mean < xmin) mean = xmin;
      if (mean > xmax) mean = xmax;
      if (++iter > 40) {
        FatalError("Search for mean failed in function CWalleniusNCHypergeometric::mean");}}
    while (fabs(mean1 - mean) > 2E-6);}
  else { // omega < 1
    do { // Newton Raphson iteration
      mean1 = mean;
      e1 = 1.-mean*m1r;
      if (e1 < 1E-14) {
        e2 = 0.;} // avoid underflow
      else {
        e2 = pow(e1,omegar-1.);}
      g = 1.-(n-mean)*m2r-e2*e1;
      gd = e2*omegar*m1r + m2r;
      mean -= g / gd;
      if (mean < xmin) mean = xmin;
      if (mean > xmax) mean = xmax;
      if (++iter > 40) {
        FatalError("Search for mean failed in function CWalleniusNCHypergeometric::mean");}}
    while (fabs(mean1 - mean) > 2E-6);}
  return mean;}


double CWalleniusNCHypergeometric::moments(double * mean_, double * var_) {
  // calculate exact mean and variance
  // return value = sum of f(x), expected = 1.
  double y, sy=0, sxy=0, sxxy=0, me1;
  long int x, xm, x1;
  const float accuracy = 1E-10f;  // accuracy of calculation
  xm = (long int)mean();  // approximation to mean
  for (x=xm; x<=xmax; x++) {
    y = probability(x);
    x1 = x - xm;  // subtract approximate mean to avoid loss of precision in sums
    sy += y; sxy += x1 * y; sxxy += x1 * x1 * y;
    if (y < accuracy) break;}
  for (x=xm-1; x>=xmin; x--) {
    y = probability(x);
    x1 = x - xm;  // subtract approximate mean to avoid loss of precision in sums
    sy += y; sxy += x1 * y; sxxy += x1 * x1 * y;
    if (y < accuracy) break;}

  me1 = sxy / sy;
  *mean_ = me1 + xm;
  y = sxxy / sy - me1 * me1;
  if (y < 0) y=0;
  *var_ = y;
  return sy;}


double CWalleniusNCHypergeometric::lnbico() {
  // natural log of binomial coefficients.
  // returns lambda = log(m!*x!/(m-x)!*m2!*x2!/(m2-x2)!)
  long int x2 = n-x, m2 = N-m;
  if (ParametersChanged) {
    mFac = LnFac(m) + LnFac(m2);
    xLast = -99; ParametersChanged = 0;}
  if (m < FAK_LEN && m2 < FAK_LEN)  goto DEFLT;
  switch (x - xLast) {
  case 0: // x unchanged
    break;
  case 1: // x incremented. calculate from previous value
    xFac += log (double(x) * (m2-x2) / (double(x2+1)*(m-x+1)));
    break;
  case -1: // x decremented. calculate from previous value
    xFac += log (double(x2) * (m-x) / (double(x+1)*(m2-x2+1)));
    break;
  default: DEFLT: // calculate all
    xFac = LnFac(x) + LnFac(x2) + LnFac(m-x) + LnFac(m2-x2);
    }
  xLast = x;
  return bico = mFac - xFac;}
  

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

  // find r to center peak of integrand at 0.5
  double dd, d1, z, zd, rr, lastr, rrc, rt, r2, r21, a, b;
  double oo[2];
  double xx[2] = {x, n-x};
  int i, j = 0;
  if (omega > 1.) { // make both omegas <= 1 to avoid overflow
    oo[0] = 1.;  oo[1] = 1./omega;}
  else {
    oo[0] = omega;  oo[1] = 1.;}

  dd = oo[0]*(m-x) + oo[1]*(N-m-xx[1]);
  d1 = 1./dd;
  E = (oo[0]*m + oo[1]*(N-m)) * d1;
  rr = r;
  if (rr <= d1) rr = 1.2*d1;           // initial guess
  // Newton-Raphson iteration to find r
  do {
    lastr = rr;
	rrc = 1. / rr;
	z = dd - rrc;
	zd = rrc * rrc;
    for (i=0; i<2; i++) {
      rt = rr * oo[i];
      if (rt < 100.) {                 // avoid overflow if rt big
        r21 = pow2_1(rt, &r2);         // r2=2^r, r21=1.-2^r
	    a = oo[i] / r21;           // omegai/(1.-2^r)
	    b = xx[i] * a;             // x*omegai/(1.-2^r)
	    z  += b;
	    zd += b * a * M_LN2 * r2;}}
    if (zd == 0) FatalError("can't find r in function CWalleniusNCHypergeometric::findpars");
	rr -= z / zd;
	if (rr <= d1) rr = lastr * 0.125 + d1*0.875;
    if (++j == 70) FatalError("convergence problem searching for r in function CWalleniusNCHypergeometric::findpars");
    }
  while (fabs(rr-lastr) > rr * 1.E-6);
  if (omega > 1) {
    dd *= omega;  rr *= oo[1];}
  r = rr;  rd = rr * dd;

  // find peak width
  double ro, k1, k2;
  ro = r * omega;
  if (ro < 300) {                      // avoid overflow
    k1 = pow2_1(ro);
    k1 = -1. / k1;
    k1 = omega*omega*(k1+k1*k1);}
  else k1 = 0.;
  if (r < 300) {                       // avoid overflow
    k2 = pow2_1(r);
    k2 = -1. / k2;
    k2 = (k2+k2*k2);}
  else k2 = 0.;
  phi2d = -4.*r*r*(x*k1 + (n-x)*k2);
  if (phi2d > 0.) FatalError("peak width undefined in function CWalleniusNCHypergeometric::findpars");
  w = 1./sqrt(-phi2d);}


/***********************************************************************
        calculation methods in class CWalleniusNCHypergeometric
***********************************************************************/

double CWalleniusNCHypergeometric::recursive() {
  // recursive calculation
  // Wallenius noncentral hypergeometric distribution by recursion formula
  // Approximate by ignoring probabilities < accuracy and minimize storage requirement
  const int BUFSIZE = 512;            // buffer size
  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;                       // 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;                    // xi_min, xi_max

  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 < x - 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 < x && p1[x2] >= accuracya) {
      x2++;  y1 = 0.;}    // increase upper limit until x has been reached
    else {
      y1 = p1[x2];}
    if (x1 > x2) return 0.;
    if (p2+x2-p > BUFSIZE) FatalError("buffer overrun in function CWalleniusNCHypergeometric::recursive");

    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;}

  if (x < x1 || x > x2) return 0.;
  return p1[x];}


double CWalleniusNCHypergeometric::binoexpand() {
  // calculate by binomial expansion of integrand
  // only for x < 2 or n-x < 2 (not implemented for higher x because of loss of precision)
  long int x1, m1, m2;