Changeset 3469

Timestamp:

Feb 29, 2016, 2:55:17 AM (6 years ago)

Author:

Peter

Message:

New class for NegativeHyperGeometric?. closes #857. The implementation
is quite naive and calls RNG several times within one call, which
might be something that can be improved. This code was inspired by
gsl_ran_hypergeometric.

Location:

trunk

Files:

• test/rnd.cc (modified) (1 diff)
• yat/random/random.cc (modified) (1 diff)
• yat/random/random.h (modified) (1 diff)

trunk/test/rnd.cc

@@ -72 +72 @@
   hg2(11, 100, 4);
 
+  NegativeHyperGeometric nhg;
+  nhg(10, 100, 4);
+  HyperGeometric nhg2(10, 100, 4);
+  nhg2();
+  // test that symmetric case doesn't get stuck
+  nhg2(11, 100, 6);
+
   return suite.return_value();
 }

trunk/yat/random/random.cc

@@ -349 +349 @@
 
 
+  NegativeHyperGeometric::NegativeHyperGeometric(void)
+  {}
+
+
+  NegativeHyperGeometric::NegativeHyperGeometric(unsigned int n1,
+                                                 unsigned int n2, unsigned int t)
+    : n1_(n1), n2_(n2), t_(t)
+  {}
+
+
+  unsigned long int NegativeHyperGeometric::operator()(void) const
+  {
+    return (*this)(n1_, n2_, t_);
+  }
+
+
+  unsigned long int NegativeHyperGeometric::operator()(unsigned int n1,
+                                                       unsigned int n2,
+                                                       unsigned int t) const
+  {
+    assert(t <= n2);
+
+    // NHG can be described as an array with n1 true and n2 false, and
+    // NHG(n1, n2, t) is the number of true left of the t:th false. By
+    // symmetry number of true right of the t:th false is NHG(n1, n2,
+    // n2-t+1) since t:th false counting from left is the (n2-t+1):th
+    // false counting from right.
+
+    // When t is larger than midpoint (2*t > n2+1) we use this
+    // symmetry to speed things up.
+    if (t > (n2+1)/2)
+      return n1 - (*this)(n1, n2, n2-t+1);
+
+    ContinuousUniform uniform;
+    unsigned long int k = 0;
+    while (t) {
+      assert(n1 + n2);
+      double x = uniform();
+      if (x * (n1+n2) < n1) {
+        --n1;
+        ++k;
+        if (!n1)
+          return k;
+      }
+      else {
+        --t;
+        --n2;
+      }
+    }
+    return k;
+  }
+
+
   Poisson::Poisson(const double m)
     : m_(m)

trunk/yat/random/random.h

@@ -524 +524 @@
 
   /**
+     We have \a n1 samples of type 1 and \a n2 samples of type 2.
+     Samples are drawn with replacement until \a t samles of type 2
+     are drawn. Then \a k, number of drawn samples of type 1, follows
+     negative hypergeometric distribution.
+
+     \since New in yat 0.14
+
+     \see HyperGeometric
+   */
+  class NegativeHyperGeometric : public Discrete
+  {
+  public:
+    /**
+       \brief Default constructor
+     */
+    NegativeHyperGeometric(void);
+
+    /**
+       \brief Constructor
+       \param n1 number of samples of type 1
+       \param n2 number of samples of type 2
+       \param t number of samples of type 2 to draw
+     */
+    NegativeHyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);
+
+    /**
+       \return random number from negative hypergeometric distribution
+       using parameters set in constructor.
+     */
+    unsigned long int operator()(void) const;
+
+    /**
+       \return random number from negative hypergeometric distribution
+       using parameters passed.
+     */
+    unsigned long int operator()(unsigned int n1, unsigned int n2,
+                                 unsigned int t) const;
+  private:
+    unsigned int n1_;
+    unsigned int n2_;
+    unsigned int t_;
+  };
+
+
+  /**
      @brief Poisson Distribution