CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandBinomial.cc
Go to the documentation of this file.
00001 // $Id: RandBinomial.cc,v 1.5 2010/06/16 17:24:53 garren Exp $
00002 // -*- C++ -*-
00003 //
00004 // -----------------------------------------------------------------------
00005 //                             HEP Random
00006 //                        --- RandBinomial ---
00007 //                      class implementation file
00008 // -----------------------------------------------------------------------
00009 
00010 // =======================================================================
00011 // John Marraffino - Created: 12th May 1998
00012 // M Fischler     - put and get to/from streams 12/10/04
00013 // M Fischler         - put/get to/from streams uses pairs of ulongs when
00014 //                      + storing doubles avoid problems with precision 
00015 //                      4/14/05
00016 //
00017 // =======================================================================
00018 
00019 #include "CLHEP/Random/RandBinomial.h"
00020 #include "CLHEP/Random/defs.h"
00021 #include "CLHEP/Random/DoubConv.hh"
00022 #include <algorithm>    // for min() and max()
00023 #include <cmath>        // for exp()
00024 
00025 using namespace std;
00026 
00027 namespace CLHEP {
00028 
00029 std::string RandBinomial::name() const {return "RandBinomial";}
00030 HepRandomEngine & RandBinomial::engine() {return *localEngine;}
00031 
00032 RandBinomial::~RandBinomial() {
00033 }
00034 
00035 double RandBinomial::shoot( HepRandomEngine *anEngine, long n,
00036                                                           double p ) {
00037   return genBinomial( anEngine, n, p );
00038 }
00039 
00040 double RandBinomial::shoot( long n, double p ) {
00041   HepRandomEngine *anEngine = HepRandom::getTheEngine();
00042   return genBinomial( anEngine, n, p );
00043 }
00044 
00045 double RandBinomial::fire( long n, double p ) {
00046   return genBinomial( localEngine.get(), n, p );
00047 }
00048 
00049 void RandBinomial::shootArray( const int size, double* vect,
00050                             long n, double p )
00051 {
00052   for( double* v = vect; v != vect+size; ++v )
00053     *v = shoot(n,p);
00054 }
00055 
00056 void RandBinomial::shootArray( HepRandomEngine* anEngine,
00057                             const int size, double* vect,
00058                             long n, double p )
00059 {
00060   for( double* v = vect; v != vect+size; ++v )
00061     *v = shoot(anEngine,n,p);
00062 }
00063 
00064 void RandBinomial::fireArray( const int size, double* vect)
00065 {
00066   for( double* v = vect; v != vect+size; ++v )
00067     *v = fire(defaultN,defaultP);
00068 }
00069 
00070 void RandBinomial::fireArray( const int size, double* vect,
00071                            long n, double p )
00072 {
00073   for( double* v = vect; v != vect+size; ++v )
00074     *v = fire(n,p);
00075 }
00076 
00077 /*************************************************************************
00078  *                                                                       *
00079  *  StirlingCorrection()                                                 *
00080  *                                                                       *
00081  *  Correction term of the Stirling approximation for log(k!)            *
00082  *  (series in 1/k, or table values for small k)                         *
00083  *  with long int parameter k                                            *
00084  *                                                                       *
00085  *************************************************************************
00086  *                                                                       *
00087  * log k! = (k + 1/2)log(k + 1) - (k + 1) + (1/2)log(2Pi) +              *
00088  *          StirlingCorrection(k + 1)                                    *
00089  *                                                                       *
00090  * log k! = (k + 1/2)log(k)     -  k      + (1/2)log(2Pi) +              *
00091  *          StirlingCorrection(k)                                        *
00092  *                                                                       *
00093  *************************************************************************/
00094 
00095 static double StirlingCorrection(long int k)
00096 {
00097   #define   C1               8.33333333333333333e-02     //  +1/12 
00098   #define   C3              -2.77777777777777778e-03     //  -1/360
00099   #define   C5               7.93650793650793651e-04     //  +1/1260
00100   #define   C7              -5.95238095238095238e-04     //  -1/1680
00101 
00102   static double  c[31] = {   0.0,
00103                              8.106146679532726e-02, 4.134069595540929e-02,
00104                              2.767792568499834e-02, 2.079067210376509e-02,
00105                              1.664469118982119e-02, 1.387612882307075e-02,
00106                              1.189670994589177e-02, 1.041126526197209e-02,
00107                              9.255462182712733e-03, 8.330563433362871e-03,
00108                              7.573675487951841e-03, 6.942840107209530e-03,
00109                              6.408994188004207e-03, 5.951370112758848e-03,
00110                              5.554733551962801e-03, 5.207655919609640e-03,
00111                              4.901395948434738e-03, 4.629153749334029e-03,
00112                              4.385560249232324e-03, 4.166319691996922e-03,
00113                              3.967954218640860e-03, 3.787618068444430e-03,
00114                              3.622960224683090e-03, 3.472021382978770e-03,
00115                              3.333155636728090e-03, 3.204970228055040e-03,
00116                              3.086278682608780e-03, 2.976063983550410e-03,
00117                              2.873449362352470e-03, 2.777674929752690e-03,
00118   };
00119   double    r, rr;
00120 
00121   if (k > 30L) {
00122     r = 1.0 / (double) k;  rr = r * r;
00123     return(r*(C1 + rr*(C3 + rr*(C5 + rr*C7))));
00124         }
00125         else  return(c[k]);
00126 }
00127 
00128 double RandBinomial::genBinomial( HepRandomEngine *anEngine, long n, double p )
00129 {
00130 /******************************************************************
00131  *                                                                *
00132  *     Binomial-Distribution - Acceptance Rejection/Inversion     *
00133  *                                                                *
00134  ******************************************************************
00135  *                                                                *
00136  * Acceptance Rejection method combined with Inversion for        *
00137  * generating Binomial random numbers with parameters             *
00138  * n (number of trials) and p (probability of success).           *
00139  * For  min(n*p,n*(1-p)) < 10  the Inversion method is applied:   *
00140  * The random numbers are generated via sequential search,        *
00141  * starting at the lowest index k=0. The cumulative probabilities *
00142  * are avoided by using the technique of chop-down.               *
00143  * For  min(n*p,n*(1-p)) >= 10  Acceptance Rejection is used:     *
00144  * The algorithm is based on a hat-function which is uniform in   *
00145  * the centre region and exponential in the tails.                *
00146  * A triangular immediate acceptance region in the centre speeds  *
00147  * up the generation of binomial variates.                        *
00148  * If candidate k is near the mode, f(k) is computed recursively  *
00149  * starting at the mode m.                                        *
00150  * The acceptance test by Stirling's formula is modified          *
00151  * according to W. Hoermann (1992): The generation of binomial    *
00152  * random variates, to appear in J. Statist. Comput. Simul.       *
00153  * If  p < .5  the algorithm is applied to parameters n, p.       *
00154  * Otherwise p is replaced by 1-p, and k is replaced by n - k.    *
00155  *                                                                *
00156  ******************************************************************
00157  *                                                                *
00158  * FUNCTION:    - btpec samples a random number from the binomial *
00159  *                distribution with parameters n and p  and is    *
00160  *                valid for  n*min(p,1-p)  >  0.                  *
00161  * REFERENCE:   - V. Kachitvichyanukul, B.W. Schmeiser (1988):    *
00162  *                Binomial random variate generation,             *
00163  *                Communications of the ACM 31, 216-222.          *
00164  * SUBPROGRAMS: - StirlingCorrection()                            *
00165  *                            ... Correction term of the Stirling *
00166  *                                approximation for log(k!)       *
00167  *                                (series in 1/k or table values  *
00168  *                                for small k) with long int k    *
00169  *              - anEngine    ... Pointer to a (0,1)-Uniform      * 
00170  *                                engine                          *
00171  *                                                                *
00172  * Implemented by H. Zechner and P. Busswald, September 1992      *
00173  ******************************************************************/
00174 
00175 #define C1_3     0.33333333333333333
00176 #define C5_8     0.62500000000000000
00177 #define C1_6     0.16666666666666667
00178 #define DMAX_KM  20L
00179 
00180   static long int      n_last = -1L,  n_prev = -1L;
00181   static double        par,np,p0,q,p_last = -1.0, p_prev = -1.0;
00182   static long          b,m,nm;
00183   static double        pq, rc, ss, xm, xl, xr, ll, lr, c,
00184                                  p1, p2, p3, p4, ch;
00185 
00186   long                 bh,i, K, Km, nK;
00187   double               f, rm, U, V, X, T, E;
00188 
00189   if (n != n_last || p != p_last)                 // set-up 
00190         {
00191          n_last = n;
00192          p_last = p;
00193          par=min(p,1.0-p);
00194          q=1.0-par;
00195          np = n*par;
00196 
00197 // Check for invalid input values
00198 
00199          if( np <= 0.0 ) return (-1.0);
00200 
00201          rm = np + par;
00202          m  = (long int) rm;                      // mode, integer 
00203          if (np<10)
00204         {
00205          p0=exp(n*log(q));                        // Chop-down
00206          bh=(long int)(np+10.0*sqrt(np*q));
00207          b=min(n,bh);
00208         }
00209          else
00210                  {
00211         rc = (n + 1.0) * (pq = par / q);          // recurr. relat.
00212         ss = np * q;                              // variance  
00213         i  = (long int) (2.195*sqrt(ss) - 4.6*q); // i = p1 - 0.5
00214         xm = m + 0.5;
00215         xl = (double) (m - i);                    // limit left 
00216         xr = (double) (m + i + 1L);               // limit right
00217         f  = (rm - xl) / (rm - xl*par);  ll = f * (1.0 + 0.5*f);
00218         f  = (xr - rm) / (xr * q);     lr = f * (1.0 + 0.5*f);
00219         c  = 0.134 + 20.5/(15.3 + (double) m);    // parallelogram
00220                                                   // height
00221         p1 = i + 0.5;
00222         p2 = p1 * (1.0 + c + c);                  // probabilities
00223         p3 = p2 + c/ll;                           // of regions 1-4
00224         p4 = p3 + c/lr;
00225                  }
00226   }
00227   if (np<10)                                      //Inversion Chop-down
00228          {
00229           double pk;
00230 
00231           K=0;
00232           pk=p0;
00233           U=anEngine->flat();
00234           while (U>pk)
00235                 {
00236                  ++K;
00237                  if (K>b)
00238                          {
00239                 U=anEngine->flat();
00240                 K=0;
00241                 pk=p0;
00242                          }
00243                  else
00244                          {
00245                 U-=pk;
00246                 pk=(double)(((n-K+1)*par*pk)/(K*q));
00247                          }
00248                 }
00249           return ((p>0.5) ? (double)(n-K):(double)K);
00250          }
00251 
00252   for (;;)
00253         {
00254          V = anEngine->flat();
00255          if ((U = anEngine->flat() * p4) <= p1)  // triangular region
00256                 {
00257                  K=(long int) (xm - U + p1*V);
00258         return ((p>0.5) ? (double)(n-K):(double)K);  // immediate accept
00259                 }
00260          if (U <= p2)                                // parallelogram
00261                 {
00262                  X = xl + (U - p1)/c;
00263                  if ((V = V*c + 1.0 - fabs(xm - X)/p1) >= 1.0)  continue;
00264                  K = (long int) X;
00265                 }
00266          else if (U <= p3)                           // left tail
00267                 {
00268                  if ((X = xl + log(V)/ll) < 0.0)  continue;
00269                  K = (long int) X;
00270                  V *= (U - p2) * ll;
00271                 }
00272          else                                         // right tail
00273                 {
00274                  if ((K = (long int) (xr - log(V)/lr)) > n)  continue;
00275                  V *= (U - p3) * lr;
00276                 }
00277 
00278  // acceptance test :  two cases, depending on |K - m|
00279          if ((Km = labs(K - m)) <= DMAX_KM || Km + Km + 2L >= ss)
00280           {
00281 
00282  // computation of p(K) via recurrence relationship from the mode
00283                 f = 1.0;                              // f(m)
00284                 if (m < K)
00285          {
00286           for (i = m; i < K; )
00287                 {
00288                 if ((f *= (rc / ++i - pq)) < V)  break;  // multiply  f
00289                 }
00290          }
00291                 else
00292          {
00293           for (i = K; i < m; )
00294                  {
00295                   if ((V *= (rc / ++i - pq)) > f)  break; // multiply  V
00296                  }
00297          }
00298                 if (V <= f)  break;                       // acceptance test
00299          }
00300   else
00301          {
00302 
00303  // lower and upper squeeze tests, based on lower bounds for log p(K)
00304                 V = log(V);
00305                 T = - Km * Km / (ss + ss);
00306                 E =  (Km / ss) * ((Km * (Km * C1_3 + C5_8) + C1_6) / ss + 0.5);
00307                 if (V <= T - E)  break;
00308                 if (V <= T + E)
00309                  {
00310         if (n != n_prev || par != p_prev)
00311          {
00312           n_prev = n;
00313           p_prev = par;
00314 
00315           nm = n - m + 1L;
00316           ch = xm * log((m + 1.0)/(pq * nm)) +
00317                StirlingCorrection(m + 1L) + StirlingCorrection(nm);
00318          }
00319         nK = n - K + 1L;
00320 
00321  // computation of log f(K) via Stirling's formula
00322  // final acceptance-rejection test
00323         if (V <= ch + (n + 1.0)*log((double) nm / (double) nK) +
00324                  (K + 0.5)*log(nK * pq / (K + 1.0)) -
00325                  StirlingCorrection(K + 1L) - StirlingCorrection(nK))  break;
00326                 }
00327          }
00328   }
00329   return ((p>0.5) ? (double)(n-K):(double)K);
00330 }
00331 
00332 std::ostream & RandBinomial::put ( std::ostream & os ) const {
00333   int pr=os.precision(20);
00334   std::vector<unsigned long> t(2);
00335   os << " " << name() << "\n";
00336   os << "Uvec" << "\n";
00337   t = DoubConv::dto2longs(defaultP);
00338   os << defaultN << " " << defaultP << " " << t[0] << " " << t[1] << "\n";
00339   os.precision(pr);
00340   return os;
00341 #ifdef REMOVED
00342   int pr=os.precision(20);
00343   os << " " << name() << "\n";
00344   os << defaultN << " " << defaultP << "\n";
00345   os.precision(pr);
00346   return os;
00347 #endif
00348 }
00349 
00350 std::istream & RandBinomial::get ( std::istream & is ) {
00351   std::string inName;
00352   is >> inName;
00353   if (inName != name()) {
00354     is.clear(std::ios::badbit | is.rdstate());
00355     std::cerr << "Mismatch when expecting to read state of a "
00356               << name() << " distribution\n"
00357               << "Name found was " << inName
00358               << "\nistream is left in the badbit state\n";
00359     return is;
00360   }
00361   if (possibleKeywordInput(is, "Uvec", defaultN)) {
00362     std::vector<unsigned long> t(2);
00363     is >> defaultN >> defaultP;
00364     is >> t[0] >> t[1]; defaultP = DoubConv::longs2double(t); 
00365     return is;
00366   }
00367   // is >> defaultN encompassed by possibleKeywordInput
00368   is >> defaultP;
00369   return is;
00370 }
00371 
00372 
00373 }  // namespace CLHEP