Changeset 424 for trunk/lib/random/random.h
- Timestamp:
- Dec 7, 2005, 11:01:17 AM (16 years ago)
- File:
-
- 1 edited
-
trunk/lib/random/random.h (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/lib/random/random.h
r411 r424 12 12 13 13 namespace theplu { 14 namespace random { 14 namespace random { 15 15 16 16 /// 17 17 /// @brief Random Number Generator 18 18 /// 19 20 21 22 19 /// The RNG class is wrapper to the GSL random number generator 20 /// (rng). This class provides a single global instance of the rng, 21 /// and makes sure there is only one point of access to the 22 /// generator. 23 23 /// 24 24 /// There is information about how to change seeding and generators … … 44 44 /// thread safety. 45 45 /// 46 /// @see Design Patterns (the singleton and adapter pattern). GSL 47 /// documentation. 46 /// @see Design Patterns (the singleton and adapter pattern). GSL 47 /// documentation. 48 /// 49 /// @todo Is this class properly implemented? The underlying random 50 /// number genereator should be a singleton, while allowing for 51 /// several different distribution to be used. Jari's feeling is 52 /// that the current implementation restricts the use to only one 53 /// distribution per binary! 48 54 /// 49 55 class RNG … … 53 59 virtual ~RNG(void); 54 60 55 static RNG* instance(u_long seed=0); 61 /// 62 /// @brief Get an instance of the random number 63 /// generator/distribution. 64 /// 65 /// @return A pointer to the random number generator. 66 /// 67 static RNG* instance(u_long seed=0) 68 { if (!instance_) instance_=new RNG(seed); return instance_; } 56 69 57 70 /// … … 81 94 /// @brief Set the seed \a s for the rng. 82 95 /// 83 /// @see seed_ dev_urandom96 /// @see seed_from_devurandom 84 97 /// 85 98 inline void seed(u_long s) const { gsl_rng_set(rng_,s); } … … 93 106 94 107 private: 95 96 108 RNG(u_long seed); 97 109 … … 100 112 }; 101 113 102 103 104 /// 105 /// @brief Continuous random number distributions. 106 /// 107 /// Abstract base class for continuous random number distributions. 108 /// 109 class Continuous 110 { 111 public: 112 113 /// 114 /// @brief Constructor 115 /// 116 inline Continuous(void) { rng_=RNG::instance(); } 117 118 /// 119 /// @return A random number 120 /// 121 virtual double operator()(void) const = 0; 122 114 // --------------------- Discrete distribtuions --------------------- 115 116 /// 117 /// @brief Discrete random number distributions. 118 /// 119 /// Abstract base class for discrete random number 120 /// distributions. Given K discrete events with different 121 /// probabilities \f$ P[k] \f$, produce a random value k consistent 122 /// with its probability. 123 /// 124 class Discrete 125 { 126 public: 127 /// 128 /// @brief Constructor 129 /// 130 inline Discrete(void) { rng_=RNG::instance(); } 131 132 /// 133 /// @return A random number. 134 /// 135 virtual u_long operator()(void) const = 0; 136 123 137 protected: 124 138 RNG* rng_; 125 139 }; 126 140 141 /// 142 /// @brief General 143 /// 144 class DiscreteGeneral : public Discrete 145 { 146 public: 147 /// 148 /// @brief Constructor 149 /// 150 /// @param hist is a Histogram defining the probability distribution 151 /// 152 DiscreteGeneral(const statistics::Histogram& hist); 153 154 /// 155 /// @brief Destructor 156 /// 157 virtual ~DiscreteGeneral(void); 158 159 /// 160 /// The generated number is an integer and proportinal to the 161 /// frequency in the corresponding histogram bin. In other words, 162 /// the probability that 0 is returned is proportinal to the size 163 /// of the first bin. 164 /// 165 /// @return A random number. 166 /// 167 inline u_long 168 operator()(void) const { return gsl_ran_discrete(rng_->rng(), gen_); } 169 170 private: 171 gsl_ran_discrete_t* gen_; 172 }; 173 174 /// 175 /// @brief Discrete uniform distribution 176 /// 177 /// Discrete uniform distribution also known as the "equally likely 178 /// outcomes" distribution. Each outcome, in this case an integer 179 /// from [0,n-1] , have equal probability to occur. 180 /// 181 /// Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le 182 /// k < n \f$ \n 183 /// Expectation value: \f$ \frac{n-1}{2} \f$ \n 184 /// Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$ 185 /// 186 class DiscreteUniform : public Discrete 187 { 188 public: 189 /// 190 /// @brief Default constructor. 191 /// 192 DiscreteUniform(void) : range_(rng_->max()+1) {} 193 194 /// 195 /// @brief Constructor. 196 /// 197 /// @param n sets the range. Object will generate integers from 198 /// [0,n-1]. 199 /// 200 DiscreteUniform(const u_long n) : range_(n) 201 { if ( range_-1>rng_->max() ) range_=rng_->max()+1; } 202 203 /// 204 /// This function returns a random integer from 0 to n-1 205 /// inclusive. All integers in the range [0,n-1] are equally 206 /// likely. n is set in constructor. 207 /// 208 inline u_long 209 operator()(void) const { return gsl_rng_uniform_int(rng_->rng(), range_); } 210 211 /// 212 /// This function returns a random integer from 0 to n-1 213 /// inclusive. All integers in the range [0,n-1] are equally 214 /// likely. 215 /// 216 inline u_long operator()(const u_long n) const 217 { return gsl_rng_uniform_int(rng_->rng(), n); } 218 219 private: 220 u_long range_; 221 }; 222 223 /// 224 /// @brief Poisson Distribution 225 /// 226 /// Having a Poisson process (i.e. no memory), number of occurences 227 /// within a given time window is Poisson distributed. This 228 /// distribution is the limit of a Binomial distribution when number 229 /// of attempts is large, and the probability for one attempt to be 230 /// succesful is small (in such a way that the expected number of 231 /// succesful attempts is \f$ m \f$. 232 /// 233 /// Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le 234 /// k \f$ \n 235 /// Expectation value: \f$ m \f$ \n 236 /// Variance: \f$ m \f$ 237 /// 238 class Poisson : public Discrete 239 { 240 public: 241 /// 242 /// @brief Constructor 243 /// 244 /// @param m is expectation value 245 inline Poisson(const double m=1) : m_(m) {} 246 247 /// 248 /// @return A Poisson distributed number. 249 /// 250 inline u_long 251 operator()(void) const { return gsl_ran_poisson(rng_->rng(), m_); } 252 253 /// 254 /// @return A Poisson distributed number with expectation value \a 255 /// m 256 /// 257 /// @note this operator ignores parameters set in Constructor 258 /// 259 inline u_long 260 operator()(const double m) const { return gsl_ran_poisson(rng_->rng(), m); } 261 262 private: 263 double m_; 264 }; 265 266 // --------------------- Continuous distribtuions --------------------- 267 268 /// 269 /// @brief Continuous random number distributions. 270 /// 271 /// Abstract base class for continuous random number distributions. 272 /// 273 class Continuous 274 { 275 public: 276 277 /// 278 /// @brief Constructor 279 /// 280 inline Continuous(void) { rng_=RNG::instance(); } 281 282 /// 283 /// @return A random number 284 /// 285 virtual double operator()(void) const = 0; 286 287 protected: 288 RNG* rng_; 289 }; 127 290 128 291 /// … … 139 302 { 140 303 public: 141 142 304 inline double operator()(void) const { return gsl_rng_uniform(rng_->rng());} 143 305 }; 306 307 /// 308 /// Class to generate numbers from a histogram in a continuous manner. 309 /// 310 class ContinuousGeneral : public Continuous 311 { 312 public: 313 /// 314 /// @brief Constructor 315 /// 316 /// @param hist is a Histogram defining the probability distribution 317 /// 318 inline ContinuousGeneral(const statistics::Histogram& hist) 319 : discrete_(DiscreteGeneral(hist)), hist_(hist) {} 320 321 /// 322 /// @brief Destructor 323 /// 324 virtual ~ContinuousGeneral(void); 325 326 /// 327 /// The number is generated in a two step process. First the bin 328 /// in the histogram is randomly selected (see 329 /// DiscreteGeneral). Then a number is generated uniformly from 330 /// the interval defined by the bin. 331 /// 332 /// @return A random number. 333 /// 334 inline double operator()(void) const 335 { return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing(); } 336 337 private: 338 const DiscreteGeneral discrete_; 339 const statistics::Histogram& hist_; 340 ContinuousUniform u_; 341 }; 342 343 /// 344 /// @brief Generator of random numbers from an exponential 345 /// distribution. 346 /// 347 /// The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a) 348 /// \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and 349 /// variance \f$ m^2 \f$ 350 /// 351 class Exponential : public Continuous 352 { 353 public: 354 /// 355 /// @brief Constructor 356 /// 357 /// @param \a m is the expectation value of the distribution. 358 /// 359 inline Exponential(const double m=1) : m_(m) {} 360 361 /// 362 /// @return A random number from exponential distribution. 363 /// 364 inline double 365 operator()(void) const { return gsl_ran_exponential(rng_->rng(), m_); } 366 367 /// 368 /// @return A random number from exponential distribution, with 369 /// expectation value \a m 370 /// 371 /// @note This operator ignores parameters given in constructor. 372 /// 373 inline double operator()(const double m) const 374 { return gsl_ran_exponential(rng_->rng(), m); } 375 376 private: 377 double m_; 144 378 }; 145 379 … … 165 399 /// m is the expectation value \f$ \mu \f$ of the distribution 166 400 /// 167 inline Gaussian(const double s=1, const double m=0) 168 : m_(m), s_(s) {} 401 inline Gaussian(const double s=1, const double m=0) : m_(m), s_(s) {} 169 402 170 403 /// 171 404 /// @return A random Gaussian number 172 405 /// 173 inline double operator()(void) const174 { return gsl_ran_gaussian(rng_->rng(), s_)+m_; }406 inline double 407 operator()(void) const { return gsl_ran_gaussian(rng_->rng(), s_)+m_; } 175 408 176 409 /// … … 179 412 /// @note this operator ignores parameters given in Constructor 180 413 /// 181 inline double operator()(const double s) const182 { return gsl_ran_gaussian(rng_->rng(), s); }414 inline double 415 operator()(const double s) const { return gsl_ran_gaussian(rng_->rng(), s); } 183 416 184 417 /// … … 187 420 /// @note this operator ignores parameters given in Constructor 188 421 /// 189 inline double operator()(const double s, const double m) const 422 inline double operator()(const double s, const double m) const 190 423 { return gsl_ran_gaussian(rng_->rng(), s)+m; } 191 424 … … 195 428 }; 196 429 197 ///198 /// @brief Exponential distribution199 ///200 /// Class for generating a random number from an Exponential201 /// distribution.202 ///203 /// Distribution function \f$ f(x) = \frac{1}{m}\exp(-x/a) \f$ for204 /// \f$ x \f$ \n205 /// Expectation value: \f$ m \f$ \n206 /// Variance: \f$ m^2 \f$207 ///208 class Exponential : public Continuous209 {210 public:211 ///212 /// @brief Constructor213 /// @param m is the expectation value of the distribution.214 ///215 inline Exponential(const double m=1)216 : m_(m){}217 218 ///219 /// @return A random number from exponential distribution220 ///221 inline double operator()(void) const222 { return gsl_ran_exponential(rng_->rng(), m_); }223 224 ///225 /// @return A random number from exponential distribution, with226 /// expectation value \a m227 /// @note this operator ignores parameters given in Constructor228 ///229 inline double operator()(const double m) const230 { return gsl_ran_exponential(rng_->rng(), m); }231 232 private:233 double m_;234 };235 236 ///237 /// @brief Discrete random number distributions.238 ///239 /// Abstract base class for discrete random number240 /// distributions. Given K discrete events with different241 /// probabilities \f$ P[k] \f$, produce a random value k consistent242 /// with its probability.243 ///244 class Discrete245 {246 public:247 ///248 /// @brief Constructor249 ///250 inline Discrete(void) { rng_=RNG::instance(); }251 252 ///253 /// @return A random number.254 ///255 virtual u_long operator()(void) const = 0;256 257 protected:258 RNG* rng_;259 };260 261 ///262 /// @brief Discrete uniform distribution263 ///264 /// Discrete uniform distribution also known as the "equally likely265 /// outcomes" distribution. Each outcome, in this case an integer266 /// from [0,n-1] , have equal probability to occur.267 ///268 /// Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le269 /// k < n \f$ \n270 /// Expectation value: \f$ \frac{n-1}{2} \f$ \n271 /// Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$272 ///273 274 class DiscreteUniform : public Discrete275 {276 public:277 ///278 /// @brief Constructor.279 ///280 DiscreteUniform(void);281 282 ///283 /// @brief Constructor.284 ///285 /// @param n sets the range. Object will generate integers from286 /// [0,n-1].287 ///288 DiscreteUniform(const u_long n);289 290 ///291 /// This function returns a random integer from 0 to n-1292 /// inclusive. All integers in the range [0,n-1] are equally293 /// likely. n is set in constructor.294 ///295 inline u_long operator()(void) const296 { return gsl_rng_uniform_int(rng_->rng(), range_); }297 298 ///299 /// This function returns a random integer from 0 to n-1300 /// inclusive. All integers in the range [0,n-1] are equally301 /// likely.302 ///303 inline u_long operator()(const u_long n) const304 { return gsl_rng_uniform_int(rng_->rng(), n); }305 306 307 private:308 u_long range_;309 };310 311 312 ///313 /// @brief Poisson Distribution314 ///315 /// Having a Poisson process (no memory), number of occurences316 /// within a given time window is Poisson distributed. This317 /// distribution is the limit of a Binomial distribution when number318 /// of attempts is large, and the probability for one attempt to be319 /// succesful is small (in such a way that the expected number of320 /// succesful attempts is \f$ m \f$.321 ///322 /// Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le323 /// k \f$ \n324 /// Expectation value: \f$ m \f$ \n325 /// Variance: \f$ m \f$326 ///327 class Poisson : public Discrete328 {329 public:330 ///331 /// @brief Constructor332 ///333 /// @param m is expectation value334 inline Poisson(const double m=1)335 : m_(m){}336 337 ///338 /// @return A Poisson distributed number.339 ///340 inline u_long operator()(void) const341 { return gsl_ran_poisson(rng_->rng(), m_); }342 343 ///344 /// @return A Poisson distributed number with expectation value \a345 /// m346 /// @note this operator ignores parameters set in Constructor347 inline u_long operator()(const double m) const348 { return gsl_ran_poisson(rng_->rng(), m); }349 350 private:351 double m_;352 };353 354 ///355 /// @brief General356 ///357 class DiscreteGeneral : public Discrete358 {359 public:360 ///361 /// @brief Constructor362 ///363 /// @param hist is a Histogram defining the probability distribution364 ///365 DiscreteGeneral(const statistics::Histogram& hist) ;366 367 ///368 /// @brief Destructor369 ///370 virtual ~DiscreteGeneral(void);371 372 ///373 /// The generated number is an integer and proportinal to the374 /// frequency in the corresponding histogram bin. In other words,375 /// the probability that 0 is returned is proportinal to the size376 /// of the first bin.377 ///378 /// @return A random number.379 ///380 inline u_long operator()(void) const381 { return gsl_ran_discrete(rng_->rng(), gen_); }382 383 private:384 gsl_ran_discrete_t* gen_;385 };386 387 388 ///389 /// Class to generate numbers from a histogram in a continuous manner.390 ///391 class ContinuousGeneral : public Continuous392 {393 public:394 ///395 /// @brief Constructor396 ///397 /// @param hist is a Histogram defining the probability distribution398 ///399 inline ContinuousGeneral(const statistics::Histogram& hist)400 : discrete_(DiscreteGeneral(hist)), hist_(hist) {}401 402 ///403 /// @brief Destructor404 ///405 virtual ~ContinuousGeneral(void);406 407 ///408 /// The number is generated in a two step process. First the bin409 /// in the histogram is randomly selected (see410 /// DiscreteGeneral). Then a number is generated uniformly from411 /// the interval defined by the bin.412 ///413 /// @return A random number.414 ///415 inline double operator()(void) const416 { return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing(); }417 418 private:419 const DiscreteGeneral discrete_;420 const statistics::Histogram& hist_;421 ContinuousUniform u_;422 };423 424 430 425 431 }} // of namespace random and namespace theplu
Note: See TracChangeset
for help on using the changeset viewer.
