Changeset 3469


Ignore:
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:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/rnd.cc

    r3465 r3469  
    7272  hg2(11, 100, 4);
    7373
     74  NegativeHyperGeometric nhg;
     75  nhg(10, 100, 4);
     76  HyperGeometric nhg2(10, 100, 4);
     77  nhg2();
     78  // test that symmetric case doesn't get stuck
     79  nhg2(11, 100, 6);
     80
    7481  return suite.return_value();
    7582}
  • trunk/yat/random/random.cc

    r3465 r3469  
    349349
    350350
     351  NegativeHyperGeometric::NegativeHyperGeometric(void)
     352  {}
     353
     354
     355  NegativeHyperGeometric::NegativeHyperGeometric(unsigned int n1,
     356                                                 unsigned int n2, unsigned int t)
     357    : n1_(n1), n2_(n2), t_(t)
     358  {}
     359
     360
     361  unsigned long int NegativeHyperGeometric::operator()(void) const
     362  {
     363    return (*this)(n1_, n2_, t_);
     364  }
     365
     366
     367  unsigned long int NegativeHyperGeometric::operator()(unsigned int n1,
     368                                                       unsigned int n2,
     369                                                       unsigned int t) const
     370  {
     371    assert(t <= n2);
     372
     373    // NHG can be described as an array with n1 true and n2 false, and
     374    // NHG(n1, n2, t) is the number of true left of the t:th false. By
     375    // symmetry number of true right of the t:th false is NHG(n1, n2,
     376    // n2-t+1) since t:th false counting from left is the (n2-t+1):th
     377    // false counting from right.
     378
     379    // When t is larger than midpoint (2*t > n2+1) we use this
     380    // symmetry to speed things up.
     381    if (t > (n2+1)/2)
     382      return n1 - (*this)(n1, n2, n2-t+1);
     383
     384    ContinuousUniform uniform;
     385    unsigned long int k = 0;
     386    while (t) {
     387      assert(n1 + n2);
     388      double x = uniform();
     389      if (x * (n1+n2) < n1) {
     390        --n1;
     391        ++k;
     392        if (!n1)
     393          return k;
     394      }
     395      else {
     396        --t;
     397        --n2;
     398      }
     399    }
     400    return k;
     401  }
     402
     403
    351404  Poisson::Poisson(const double m)
    352405    : m_(m)
  • trunk/yat/random/random.h

    r3465 r3469  
    524524
    525525  /**
     526     We have \a n1 samples of type 1 and \a n2 samples of type 2.
     527     Samples are drawn with replacement until \a t samles of type 2
     528     are drawn. Then \a k, number of drawn samples of type 1, follows
     529     negative hypergeometric distribution.
     530
     531     \since New in yat 0.14
     532
     533     \see HyperGeometric
     534   */
     535  class NegativeHyperGeometric : public Discrete
     536  {
     537  public:
     538    /**
     539       \brief Default constructor
     540     */
     541    NegativeHyperGeometric(void);
     542
     543    /**
     544       \brief Constructor
     545       \param n1 number of samples of type 1
     546       \param n2 number of samples of type 2
     547       \param t number of samples of type 2 to draw
     548     */
     549    NegativeHyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);
     550
     551    /**
     552       \return random number from negative hypergeometric distribution
     553       using parameters set in constructor.
     554     */
     555    unsigned long int operator()(void) const;
     556
     557    /**
     558       \return random number from negative hypergeometric distribution
     559       using parameters passed.
     560     */
     561    unsigned long int operator()(unsigned int n1, unsigned int n2,
     562                                 unsigned int t) const;
     563  private:
     564    unsigned int n1_;
     565    unsigned int n2_;
     566    unsigned int t_;
     567  };
     568
     569
     570  /**
    526571     @brief Poisson Distribution
    527572
Note: See TracChangeset for help on using the changeset viewer.