Start line:  
End line:  

Snippet Preview

Snippet HTML Code

Stack Overflow Questions
  /*
  Copyright � 1999 CERN - European Organization for Nuclear Research.
  Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
  is hereby granted without fee, provided that the above copyright notice appear in all copies and
  that both that copyright notice and this permission notice appear in supporting documentation.
  CERN makes no representations about the suitability of this software for any purpose.
  It is provided "as is" without expressed or implied warranty.
  */
  package org.apache.mahout.math.jet.random;
 
 
Partially deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
 
 public final class Poisson extends AbstractDiscreteDistribution {
 
   private final double mean;
 
   // precomputed and cached values (for performance only)
   // cache for < SWITCH_MEAN
   private double myOld = -1.0;
   private double p;
   private double q;
   private double p0;
   private final double[] pp = new double[36];
   private int llll;
 
   // cache for >= SWITCH_MEAN
   private double myLast = -1.0;
   private double ll;
   private int k2;
   private int k4;
   private int k1;
   private int k5;
   private double dl;
   private double dr;
   private double r1;
   private double r2;
   private double r4;
   private double r5;
   private double lr;
   private double lMy;
   private double cPm;
   private double f1;
   private double f2;
   private double f4;
   private double f5;
   private double p1;
   private double p2;
   private double p3;
   private double p4;
   private double p5;
   private double p6;
 
   // cache for both;
 
 
   private static final double MEAN_MAX = .;
   // for all means larger than that, we don't try to compute a poisson deviation, but return the mean.
   private static final double SWITCH_MEAN = 10.0; // switch from method A to method B
 

  
Constructs a poisson distribution. Example: mean=1.0.
 
   public Poisson(double meanRandom randomGenerator) {
     setRandomGenerator(randomGenerator);
     this. = mean;
   }
 
   private static double f(int kdouble lNudouble cPm) {
     return Math.exp(k * lNu - Arithmetic.logFactorial(k) - cPm);
   }
 
   @Override
   public int nextInt() {
     return nextInt();
   }

  
Returns a random number from the distribution; bypasses the internal state.
 
   public int nextInt(double theMean) {
    
* Poisson Distribution - Patchwork Rejection/Inversion * * * For parameter my < 10 Tabulated Inversion is applied. * For my >= 10 Patchwork Rejection is employed: * The area below the histogram function f(x) is rearranged in * its body by certain point reflections. Within a large center * interval variates are sampled efficiently by rejection from * uniform hats. Rectangular immediate acceptance regions speed * up the generation. The remaining tails are covered by * exponential functions. * * /
 
     Random gen = getRandomGenerator();
 
     //double t, g, my_k;
    //double gx, gy, px, py, e, x, xx, delta, v;
    //int sign;
    //static double p,q,p0,pp[36];
    //static long ll,m;
    int m;
    if (theMean < ) { // CASE B: Inversion- start new table and calculate p0
      if (theMean != ) {
         = theMean;
         = 0;
         = Math.exp(-theMean);
         = ;
         = ;
        //for (k=pp.length; --k >=0;) pp[k] = 0;
      }
      m = theMean > 1.0 ? (inttheMean : 1;
      while (true) {
        double u = gen.nextDouble();
        int k = 0;
        if (u <= ) {
          return k;
        }
        if ( != 0) {              // Step T. Table comparison
          int i = u > 0.458 ? Math.min(m) : 1;
          for (k = ik <= k++) {
            if (u <= [k]) {
              return k;
            }
          }
          if ( == 35) {
            continue;
          }
        }
        for (k =  + 1; k <= 35; k++) { // Step C. Creation of new prob.
           *= theMean / k;
           += ;
          [k] = ;
          if (u <= ) {
             = k;
            return k;
          }
        }
         = 35;
      }
      // end my < SWITCH_MEAN
    } else if (theMean < ) { // CASE A: acceptance complement
      //static double        my_last = -1.0;
      //static long int      m,  k2, k4, k1, k5;
      //static double        dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
      //             f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
      m = (inttheMean;
      if (theMean != ) { //  set-up
         = theMean;
        // approximate deviation of reflection points k2, k4 from my - 1/2
        double Ds = Math.sqrt(theMean + 0.25);
        // mode m, reflection points k2 and k4, and points k1 and k5, which
        // delimit the centre region of h(x)
         = (int) Math.ceil(theMean - 0.5 - Ds);
         = (int) (theMean - 0.5 + Ds);
         =  +  - m + 1;
         =  +  - m;
        // range width of the critical left and right centre region
         =  - ;
         =  - ;
        // recurrence constants r(k) = p(k)/p(k-1) at k = k1, k2, k4+1, k5+1
         = theMean / ;
         = theMean / ;
         = theMean / ( + 1);
         = theMean / ( + 1);
        // reciprocal values of the scale parameters of expon. tail envelopes
         = Math.log();                     // expon. tail left
         = -Math.log();                     // expon. tail right
        // Poisson constants, necessary for computing function values f(k)
         = Math.log(theMean);
         = m *  - Arithmetic.logFactorial(m);
        // function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5
         = f();
         = f();
         = f();
         = f();
        // area of the two centre and the two exponential tail regions
        // area of the two immediate acceptance regions between k2, k4
         =  * ( + 1.0);                    // immed. left
         =  *  + ;               // centre left
         =  * ( + 1.0) + ;               // immed. right
         =  *  + ;               // centre right
         =  /  + ;               // expon. tail left
         =  /  + ;               // expon. tail right
      } // end set-up
      while (true) {
        // generate uniform number U -- U(0, p6)
        // case distinction corresponding to U
        double W;
        double V;
        double U;
        int Y;
        int X;
        int Dk;
        if ((U = gen.nextDouble() * ) < ) {         // centre left
          // immediate acceptance region R2 = [k2, m) *[0, f2),  X = k2, ... m -1
          if ((V = U - ) < 0.0) {
            return  + (int) (U / );
          }
          // immediate acceptance region R1 = [k1, k2)*[0, f1),  X = k1, ... k2-1
          if ((W = V / ) < ) {
            return  + (int) (V / );
          }
          // computation of candidate X < k2, and its counterpart Y > k2
          // either squeeze-acceptance of X or acceptance-rejection of Y
          Dk = gen.nextInt((int) + 1;
          if (W <=  - Dk * ( -  / )) {            // quick accept of
            return  - Dk;                          // X = k2 - Dk
          }
          if ((V =  +  - W) < 1.0) {                // quick reject of Y
            Y =  + Dk;
            if (V <=  + Dk * (1.0 - ) / ( + 1.0)) { // quick accept of
              return Y;                             // Y = k2 + Dk
            }
            if (V <= f(Y)) {
              return Y;
            }    // final accept of Y
          }
          X =  - Dk;
        } else if (U < ) {                                 // centre right
          // immediate acceptance region R3 = [m, k4+1)*[0, f4), X = m, ... k4
          if ((V = U - ) < 0.0) {
            return  - (int) ((U - ) / );
          }
          // immediate acceptance region R4 = [k4+1, k5+1)*[0, f5)
          if ((W = V / ) < ) {
            return  - (int) (V / );
          }
          // computation of candidate X > k4, and its counterpart Y < k4
          // either squeeze-acceptance of X or acceptance-rejection of Y
          Dk = gen.nextInt((int) + 1;
          if (W <=  - Dk * ( -  * )) {             // quick accept of
            return  + Dk;                           // X = k4 + Dk
          }
          if ((V =  +  - W) < 1.0) {                 // quick reject of Y
            Y =  - Dk;
            if (V <=  + Dk * (1.0 - ) / ) {       // quick accept of
              return Y;                             // Y = k4 - Dk
            }
            if (V <= f(Y)) {
              return Y;
            }    // final accept of Y
          }
          X =  + Dk;
        } else {
          W = gen.nextDouble();
          if (U < ) {                                  // expon. tail left
            Dk = (int) (1.0 - Math.log(W) / );
            if ((X =  - Dk) < 0) {
              continue;
            }          // 0 <= X <= k1 - 1
            W *= (U - ) * ;                        // W -- U(0, h(x))
            if (W <=  - Dk * ( -  / )) {
              return X;
            } // quick accept of X
          } else {                                         // expon. tail right
            Dk = (int) (1.0 - Math.log(W) / );
            X =  + Dk;                              // X >= k5 + 1
            W *= (U - ) * ;                        // W -- U(0, h(x))
            if (W <=  - Dk * ( -  * )) {
              return X;
            } // quick accept of X
          }
        }
        // acceptance-rejection test of candidate X from the original area
        // test, whether  W <= f(k),    with  W = U*h(x)  and  U -- U(0, 1)
        // log f(X) = (X - m)*log(my) - log X! + log m!
        if (Math.log(W) <= X *  - Arithmetic.logFactorial(X) - ) {
          return X;
        }
      }
    } else { // mean is too large
      return (inttheMean;
    }
  }
New to GrepCode? Check out our FAQ X