Changeset 3469

Ignore:
Timestamp:
Feb 29, 2016, 2:55:17 AM (6 years ago)
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

Unmodified
Removed
• trunk/test/rnd.cc

 r3465 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

 r3465 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

 r3465 /** 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
Note: See TracChangeset for help on using the changeset viewer.