CLHEP::RandBinomial Class Reference

#include <RandBinomial.h>

Inheritance diagram for CLHEP::RandBinomial:

Collaboration diagram for CLHEP::RandBinomial:

Public Member Functions
	RandBinomial (HepRandomEngine &anEngine, long n=1, double p=0.5)
	RandBinomial (HepRandomEngine *anEngine, long n=1, double p=0.5)
virtual	~RandBinomial ()
double	fire ()
double	fire (long n, double p)
void	fireArray (const int size, double *vect)
void	fireArray (const int size, double *vect, long n, double p)
double	operator() ()
double	operator() (long n, double p)
std::ostream &	put (std::ostream &os) const
std::istream &	get (std::istream &is)
std::string	name () const
HepRandomEngine &	engine ()
Public Member Functions inherited from CLHEP::HepRandom
	HepRandom ()
	HepRandom (long seed)
	HepRandom (HepRandomEngine &algorithm)
	HepRandom (HepRandomEngine *algorithm)
virtual	~HepRandom ()
double	flat ()
void	flatArray (const int size, double *vect)
double	flat (HepRandomEngine *theNewEngine)
void	flatArray (HepRandomEngine *theNewEngine, const int size, double *vect)

Static Public Member Functions
static double	shoot ()
static double	shoot (long n, double p)
static void	shootArray (const int size, double *vect, long n=1, double p=0.5)
static double	shoot (HepRandomEngine *anEngine)
static double	shoot (HepRandomEngine *anEngine, long n, double p)
static void	shootArray (HepRandomEngine *anEngine, const int size, double *vect, long n=1, double p=0.5)
static std::string	distributionName ()
Static Public Member Functions inherited from CLHEP::HepRandom
static void	setTheSeed (long seed, int lux=3)
static long	getTheSeed ()
static void	setTheSeeds (const long *seeds, int aux=-1)
static const long *	getTheSeeds ()
static void	getTheTableSeeds (long *seeds, int index)
static HepRandom *	getTheGenerator ()
static void	setTheEngine (HepRandomEngine *theNewEngine)
static HepRandomEngine *	getTheEngine ()
static void	saveEngineStatus (const char filename[]="Config.conf")
static void	restoreEngineStatus (const char filename[]="Config.conf")
static std::ostream &	saveFullState (std::ostream &os)
static std::istream &	restoreFullState (std::istream &is)
static std::ostream &	saveDistState (std::ostream &os)
static std::istream &	restoreDistState (std::istream &is)
static std::ostream &	saveStaticRandomStates (std::ostream &os)
static std::istream &	restoreStaticRandomStates (std::istream &is)
static void	showEngineStatus ()
static int	createInstance ()
static std::string	distributionName ()

Static Private Member Functions
static double	genBinomial (HepRandomEngine *anEngine, long n, double p)

Private Attributes
std::shared_ptr< HepRandomEngine >	localEngine
long	defaultN
double	defaultP

Additional Inherited Members
Static Protected Attributes inherited from CLHEP::HepRandom
static const long	seedTable [215][2]

Detailed Description

Author

Definition at line 37 of file RandBinomial.h.

Constructor & Destructor Documentation

RandBinomial() [1/2]

CLHEP::RandBinomial::RandBinomial ( HepRandomEngine &  anEngine, long  n = 1, double  p = 0.5  )

inline

RandBinomial() [2/2]

CLHEP::RandBinomial::RandBinomial ( HepRandomEngine *  anEngine, long  n = 1, double  p = 0.5  )

inline

~RandBinomial()

CLHEP::RandBinomial::~RandBinomial ( )

virtual

Definition at line 30 of file RandBinomial.cc.

                            {
}

Member Function Documentation

distributionName()

static std::string CLHEP::RandBinomial::distributionName ( )

inlinestatic

Definition at line 99 of file RandBinomial.h.

{return "RandBinomial";}  

Here is the call graph for this function:

engine()

HepRandomEngine & CLHEP::RandBinomial::engine ( )

virtual

Reimplemented from CLHEP::HepRandom.

Definition at line 28 of file RandBinomial.cc.

{return *localEngine;}

fire() [1/2]

double CLHEP::RandBinomial::fire ( )

inline

Here is the caller graph for this function:

fire() [2/2]

double CLHEP::RandBinomial::fire	(	long	n,
		double	p
	)

Definition at line 43 of file RandBinomial.cc.

                                            {
  return genBinomial( localEngine.get(), n, p );
}

Here is the call graph for this function:

fireArray() [1/2]

void CLHEP::RandBinomial::fireArray	(	const int	size,
		double *	vect
	)

Definition at line 62 of file RandBinomial.cc.

{
  for( double* v = vect; v != vect+size; ++v )
    *v = fire(defaultN,defaultP);
}

Here is the call graph for this function:

fireArray() [2/2]

void CLHEP::RandBinomial::fireArray	(	const int	size,
		double *	vect,
		long	n,
		double	p
	)

Definition at line 68 of file RandBinomial.cc.

{
  for( double* v = vect; v != vect+size; ++v )
    *v = fire(n,p);
}

Here is the call graph for this function:

genBinomial()

double CLHEP::RandBinomial::genBinomial ( HepRandomEngine *  anEngine, long  n, double  p  )

staticprivate

Definition at line 126 of file RandBinomial.cc.

{
/******************************************************************
 *                                                                *
 *     Binomial-Distribution - Acceptance Rejection/Inversion     *
 *                                                                *
 ******************************************************************
 *                                                                *
 * Acceptance Rejection method combined with Inversion for        *
 * generating Binomial random numbers with parameters             *
 * n (number of trials) and p (probability of success).           *
 * For  min(n*p,n*(1-p)) < 10  the Inversion method is applied:   *
 * The random numbers are generated via sequential search,        *
 * starting at the lowest index k=0. The cumulative probabilities *
 * are avoided by using the technique of chop-down.               *
 * For  min(n*p,n*(1-p)) >= 10  Acceptance Rejection is used:     *
 * The algorithm is based on a hat-function which is uniform in   *
 * the centre region and exponential in the tails.                *
 * A triangular immediate acceptance region in the centre speeds  *
 * up the generation of binomial variates.                        *
 * If candidate k is near the mode, f(k) is computed recursively  *
 * starting at the mode m.                                        *
 * The acceptance test by Stirling's formula is modified          *
 * according to W. Hoermann (1992): The generation of binomial    *
 * random variates, to appear in J. Statist. Comput. Simul.       *
 * If  p < .5  the algorithm is applied to parameters n, p.       *
 * Otherwise p is replaced by 1-p, and k is replaced by n - k.    *
 *                                                                *
 ******************************************************************
 *                                                                *
 * FUNCTION:    - btpec samples a random number from the binomial *
 *                distribution with parameters n and p  and is    *
 *                valid for  n*min(p,1-p)  >  0.                  *
 * REFERENCE:   - V. Kachitvichyanukul, B.W. Schmeiser (1988):    *
 *                Binomial random variate generation,             *
 *                Communications of the ACM 31, 216-222.          *
 * SUBPROGRAMS: - StirlingCorrection()                            *
 *                            ... Correction term of the Stirling *
 *                                approximation for log(k!)       *
 *                                (series in 1/k or table values  *
 *                                for small k) with long int k    *
 *              - anEngine    ... Pointer to a (0,1)-Uniform      * 
 *                                engine                          *
 *                                                                *
 * Implemented by H. Zechner and P. Busswald, September 1992      *
 ******************************************************************/

#define C1_3     0.33333333333333333
#define C5_8     0.62500000000000000
#define C1_6     0.16666666666666667
#define DMAX_KM  20L

  static CLHEP_THREAD_LOCAL long int      n_last = -1L,  n_prev = -1L;
  static CLHEP_THREAD_LOCAL double        par,np,p0,q,p_last = -1.0, p_prev = -1.0;
  static CLHEP_THREAD_LOCAL long          b,m,nm;
  static CLHEP_THREAD_LOCAL double        pq, rc, ss, xm, xl, xr, ll, lr, c,
                 p1, p2, p3, p4, ch;

  long                 bh,i, K, Km, nK;
  double               f, rm, U, V, X, T, E;

  if (n != n_last || p != p_last)                 // set-up 
    {
     n_last = n;
     p_last = p;
     par=std::min(p,1.0-p);
     q=1.0-par;
     np = n*par;

// Check for invalid input values

     if( np <= 0.0 ) return (-1.0);

     rm = np + par;
     m  = (long int) rm;                      // mode, integer 
     if (np<10)
    {
     p0=std::exp(n*std::log(q));              // Chop-down
     bh=(long int)(np+10.0*std::sqrt(np*q));
     b=std::min(n,bh);
    }
     else
         {
    rc = (n + 1.0) * (pq = par / q);          // recurr. relat.
    ss = np * q;                              // variance  
    i  = (long int) (2.195*std::sqrt(ss) - 4.6*q); // i = p1 - 0.5
    xm = m + 0.5;
    xl = (double) (m - i);                    // limit left 
    xr = (double) (m + i + 1L);               // limit right
    f  = (rm - xl) / (rm - xl*par);  ll = f * (1.0 + 0.5*f);
    f  = (xr - rm) / (xr * q);     lr = f * (1.0 + 0.5*f);
    c  = 0.134 + 20.5/(15.3 + (double) m);    // parallelogram
                          // height
    p1 = i + 0.5;
    p2 = p1 * (1.0 + c + c);                  // probabilities
    p3 = p2 + c/ll;                           // of regions 1-4
    p4 = p3 + c/lr;
         }
  }
  if( np <= 0.0 ) return (-1.0);
  if (np<10)                                      //Inversion Chop-down
     {
      double pk;

      K=0;
      pk=p0;
      U=anEngine->flat();
      while (U>pk)
        {
         ++K;
         if (K>b)
             {
        U=anEngine->flat();
        K=0;
        pk=p0;
             }
         else
             {
        U-=pk;
        pk=(double)(((n-K+1)*par*pk)/(K*q));
             }
        }
      return ((p>0.5) ? (double)(n-K):(double)K);
     }

  for (;;)
    {
     V = anEngine->flat();
     if ((U = anEngine->flat() * p4) <= p1)  // triangular region
        {
         K=(long int) (xm - U + p1*V);
    return ((p>0.5) ? (double)(n-K):(double)K);  // immediate accept
        }
     if (U <= p2)                                // parallelogram
        {
         X = xl + (U - p1)/c;
         if ((V = V*c + 1.0 - std::fabs(xm - X)/p1) >= 1.0)  continue;
         K = (long int) X;
        }
     else if (U <= p3)                           // left tail
        {
         if ((X = xl + std::log(V)/ll) < 0.0)  continue;
         K = (long int) X;
         V *= (U - p2) * ll;
        }
     else                                         // right tail
        {
         if ((K = (long int) (xr - std::log(V)/lr)) > n)  continue;
         V *= (U - p3) * lr;
        }

 // acceptance test :  two cases, depending on |K - m|
     if ((Km = std::labs(K - m)) <= DMAX_KM || Km + Km + 2L >= ss)
      {

 // computation of p(K) via recurrence relationship from the mode
        f = 1.0;                              // f(m)
        if (m < K)
     {
      for (i = m; i < K; )
        {
        if ((f *= (rc / ++i - pq)) < V)  break;  // multiply  f
        }
     }
        else
     {
      for (i = K; i < m; )
         {
          if ((V *= (rc / ++i - pq)) > f)  break; // multiply  V
         }
     }
        if (V <= f)  break;                       // acceptance test
     }
  else
     {

 // lower and upper squeeze tests, based on lower bounds for log p(K)
        V = std::log(V);
        T = - Km * Km / (ss + ss);
        E =  (Km / ss) * ((Km * (Km * C1_3 + C5_8) + C1_6) / ss + 0.5);
        if (V <= T - E)  break;
        if (V <= T + E)
         {
    if (n != n_prev || par != p_prev)
     {
      n_prev = n;
      p_prev = par;

      nm = n - m + 1L;
      ch = xm * std::log((m + 1.0)/(pq * nm)) +
           StirlingCorrection(m + 1L) + StirlingCorrection(nm);
     }
    nK = n - K + 1L;

 // computation of log f(K) via Stirling's formula
 // final acceptance-rejection test
    if (V <= ch + (n + 1.0)*std::log((double) nm / (double) nK) +
                 (K + 0.5)*std::log(nK * pq / (K + 1.0)) -
                 StirlingCorrection(K + 1L) - StirlingCorrection(nK))  break;
        }
     }
  }
  return ((p>0.5) ? (double)(n-K):(double)K);
}

Here is the call graph for this function:

Here is the caller graph for this function:

get()

std::istream & CLHEP::RandBinomial::get ( std::istream &  is )

virtual

Reimplemented from CLHEP::HepRandom.

Definition at line 342 of file RandBinomial.cc.

                                                 {
  std::string inName;
  is >> inName;
  if (inName != name()) {
    is.clear(std::ios::badbit | is.rdstate());
    std::cerr << "Mismatch when expecting to read state of a "
              << name() << " distribution\n"
          << "Name found was " << inName
          << "\nistream is left in the badbit state\n";
    return is;
  }
  if (possibleKeywordInput(is, "Uvec", defaultN)) {
    std::vector<unsigned long> t(2);
    is >> defaultN >> defaultP;
    is >> t[0] >> t[1]; defaultP = DoubConv::longs2double(t); 
    return is;
  }
  // is >> defaultN encompassed by possibleKeywordInput
  is >> defaultP;
  return is;
}

Here is the call graph for this function:

name()

std::string CLHEP::RandBinomial::name ( ) const

virtual

Reimplemented from CLHEP::HepRandom.

Definition at line 27 of file RandBinomial.cc.

{return "RandBinomial";}

Here is the caller graph for this function:

operator()() [1/2]

double CLHEP::RandBinomial::operator() ( )

inlinevirtual

Reimplemented from CLHEP::HepRandom.

operator()() [2/2]

double CLHEP::RandBinomial::operator() ( long  n, double  p  )

inline

put()

std::ostream & CLHEP::RandBinomial::put ( std::ostream &  os ) const

virtual

Reimplemented from CLHEP::HepRandom.

Definition at line 331 of file RandBinomial.cc.

                                                       {
  int pr=os.precision(20);
  std::vector<unsigned long> t(2);
  os << " " << name() << "\n";
  os << "Uvec" << "\n";
  t = DoubConv::dto2longs(defaultP);
  os << defaultN << " " << defaultP << " " << t[0] << " " << t[1] << "\n";
  os.precision(pr);
  return os;
}

Here is the call graph for this function:

shoot() [1/4]

static double CLHEP::RandBinomial::shoot ( )

inlinestatic

Here is the caller graph for this function:

shoot() [2/4]

double CLHEP::RandBinomial::shoot ( long  n, double  p  )

static

Definition at line 38 of file RandBinomial.cc.

                                             {
  HepRandomEngine *anEngine = HepRandom::getTheEngine();
  return genBinomial( anEngine, n, p );
}

Here is the call graph for this function:

shoot() [3/4]

static double CLHEP::RandBinomial::shoot ( HepRandomEngine *  anEngine )

inlinestatic

shoot() [4/4]

double CLHEP::RandBinomial::shoot ( HepRandomEngine *  anEngine, long  n, double  p  )

static

Definition at line 33 of file RandBinomial.cc.

                                                                     {
  return genBinomial( anEngine, n, p );
}

Here is the call graph for this function:

shootArray() [1/2]

void CLHEP::RandBinomial::shootArray ( const int  size, double *  vect, long  n = 1, double  p = 0.5  )

static

Definition at line 47 of file RandBinomial.cc.

{
  for( double* v = vect; v != vect+size; ++v )
    *v = shoot(n,p);
}

Here is the call graph for this function:

shootArray() [2/2]

void CLHEP::RandBinomial::shootArray ( HepRandomEngine *  anEngine, const int  size, double *  vect, long  n = 1, double  p = 0.5  )

static

Definition at line 54 of file RandBinomial.cc.

{
  for( double* v = vect; v != vect+size; ++v )
    *v = shoot(anEngine,n,p);
}

Here is the call graph for this function:

Member Data Documentation

defaultN

long CLHEP::RandBinomial::defaultN

private

Definition at line 107 of file RandBinomial.h.

defaultP

double CLHEP::RandBinomial::defaultP

private

Definition at line 108 of file RandBinomial.h.

localEngine

std::shared_ptr<HepRandomEngine> CLHEP::RandBinomial::localEngine

private

Definition at line 106 of file RandBinomial.h.

---

The documentation for this class was generated from the following files:

• RandBinomial.h
• RandBinomial.cc