source: trunk/yat/random/random.h @ 820

Last change on this file since 820 was 820, checked in by Peter, 15 years ago

Changed name again. ROCScore is now called AUC. Also fixed several bugs and added test for these.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 13.6 KB
Line 
1#ifndef _theplu_yat_random_
2#define _theplu_yat_random_
3
4// $Id: random.h 820 2007-03-17 21:54:52Z peter $
5
6/*
7  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson
8
9  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 2 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
26
27#include "yat/statistics/Histogram.h"
28
29#include <gsl/gsl_rng.h>
30#include <gsl/gsl_randist.h>
31
32#include <string>
33
34namespace theplu {
35namespace yat {
36namespace random {
37
38  //forward declarion
39  class RNG_state;
40
41  ///
42  /// @brief Random Number Generator
43  ///
44  /// The RNG class is wrapper to the GSL random number generator
45  /// (rng). This class provides a single global instance of the rng,
46  /// and makes sure there is only one point of access to the
47  /// generator.
48  ///
49  /// There is information about how to change seeding and generators
50  /// at run time without recompilation using environment variables in
51  /// the GSL manual (Chapter on random number generators). RNG of
52  /// course support seeding at compile time if you don't want to
53  /// bother about environment variables and GSL.
54  ///
55  /// There are many different rng's available in GSL. Currently only
56  /// the default generator is implemented and no other one is
57  /// choosable through the class interface. This means that you have
58  /// to fall back to the use of environment variables as described in
59  /// the GSL documentation, or be bold and request support for other
60  /// rng's through the class interface.
61  ///
62  /// Not all GSL functionality is implemented, we'll add
63  /// functionality when needed and may do it when requested. Better
64  /// yet, supply us with code and we will probably add it to the code
65  /// (BUT remember to implement reasonable tests for your code and
66  /// follow the coding style.)
67  ///
68  /// The current implementation is NOT thread safe since the RNG is
69  /// implemented as a singleton. However, the underlying GSL rng's
70  /// support thread safety since each instance of GSL rng's keep
71  /// track of their own state accordning to GSL documentation.
72  ///
73  /// @see Design Patterns (the singleton and adapter pattern). GSL
74  /// documentation.
75  ///
76  class RNG
77  {
78  public:
79
80    virtual ~RNG(void);
81
82    ///
83    /// @brief Get an instance of the random number generator.
84    ///
85    /// Get an instance of the random number generator. If the random
86    /// number generator is not already created, the call will create
87    /// a new generator and use the default seed. The seed must be
88    /// changed with the seed or seed_from_devurandom member
89    /// functions.
90    ///
91    /// @return A pointer to the random number generator.
92    ///
93    /// @see seed and seed_from_devurandom
94    ///
95    static RNG* instance(void);
96
97    ///
98    /// @brief Returns the largest number that the random number
99    /// generator can return.
100    ///
101    u_long max(void) const;
102
103    ///
104    /// @brief Returns the smallest number that the random number
105    /// generator can return.
106    ///
107    u_long min(void) const;
108
109    ///
110    /// @brief Returns the name of the random number generator
111    ///
112    std::string name(void) const;
113
114    ///
115    /// @return const pointer to underlying GSL random generator.
116    ///
117    const gsl_rng* rng(void) const;
118
119    ///
120    /// @brief Set the seed \a s for the rng.
121    ///
122    /// Set the seed \a s for the rng. If \a s is zero, a default
123    /// value from the rng's original implementation is used (cf. GSL
124    /// documentation).
125    ///
126    /// @see seed_from_devurandom
127    ///
128    void seed(u_long s) const;
129
130    ///
131    /// @brief Seed the rng using the /dev/urandom device.
132    ///
133    /// @return The seed acquired from /dev/urandom.
134    ///
135    u_long seed_from_devurandom(void);
136
137    /**
138       \brief Set the state to \a state.
139
140       \return 0 on success, non-zero otherwise.
141
142       \see gsl_rng_memcpy
143    */
144    int set_state(const RNG_state&);
145
146  private:
147    RNG(void);
148
149    static RNG* instance_;
150    gsl_rng* rng_;
151  };
152
153
154  ///
155  /// @brief Class holding state of a random generator
156  ///
157  class RNG_state
158  {
159  public:
160    ///
161    /// @brief Constructor
162    ///
163    RNG_state(const RNG*);
164
165    ///
166    /// @brief Destructor
167    ///
168    ~RNG_state(void);
169
170    ///
171    /// @return const pointer to underlying GSL random generator.
172    ///
173    const gsl_rng* rng(void) const;
174
175  private:
176    gsl_rng* rng_;
177
178  };
179   
180
181  // --------------------- Discrete distribtuions ---------------------
182
183  ///
184  /// @brief Discrete random number distributions.
185  ///
186  /// Abstract base class for discrete random number
187  /// distributions. Given K discrete events with different
188  /// probabilities \f$ P[k] \f$, produce a random value k consistent
189  /// with its probability.
190  ///
191  class Discrete
192  {
193  public:
194    ///
195    /// @brief Constructor
196    ///
197    Discrete(void);
198
199    ///
200    /// @brief The destructor
201    ///
202    virtual ~Discrete(void);
203
204    ///
205    /// @brief Set the seed to \a s.
206    ///
207    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
208    /// default value from the rng's original implementation is used
209    /// (cf. GSL documentation).
210    ///
211    /// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
212    ///
213    void seed(u_long s) const;
214
215    ///
216    /// @brief Set the seed using the /dev/urandom device.
217    ///
218    /// @return The seed acquired from /dev/urandom.
219    ///
220    /// @see seed, RNG::seed_from_devurandom, RNG::seed
221    ///
222    u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
223
224    ///
225    /// @return A random number.
226    ///
227    virtual u_long operator()(void) const = 0;
228   
229  protected:
230    /// GSL random gererator
231    RNG* rng_;
232  };
233
234  ///
235  /// @brief General
236  ///
237  class DiscreteGeneral : public Discrete
238  {
239  public:
240    ///
241    /// @brief Constructor
242    ///
243    /// @param hist is a Histogram defining the probability distribution
244    ///
245    DiscreteGeneral(const statistics::Histogram& hist);
246   
247    ///
248    /// @brief Destructor
249    ///
250    ~DiscreteGeneral(void);
251
252    ///
253    /// The generated number is an integer and proportinal to the
254    /// frequency in the corresponding histogram bin. In other words,
255    /// the probability that 0 is returned is proportinal to the size
256    /// of the first bin.
257    ///
258    /// @return A random number.
259    ///
260    u_long operator()(void) const;
261
262  private:
263     gsl_ran_discrete_t* gen_;
264  };
265
266  ///
267  /// @brief Discrete uniform distribution
268  ///
269  /// Discrete uniform distribution also known as the "equally likely
270  /// outcomes" distribution. Each outcome, in this case an integer
271  /// from [0,n-1] , have equal probability to occur.
272  ///
273  /// Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le
274  /// k < n \f$ \n
275  /// Expectation value: \f$ \frac{n-1}{2} \f$ \n
276  /// Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$
277  ///
278  class DiscreteUniform : public Discrete
279  {
280  public:
281    /**
282       \brief Constructor.
283
284       The generator will generate integers within the range \f$
285       [0,n-1] \f$. If \a n is zero, then the whole range of the
286       underlying RNG will be used \f$ [min,max] \f$. Setting \a n to
287       zero is the preferred way to sample the whole range of the
288       underlying RNG, i.e. not setting \n to RNG.max.
289
290       \throw If \a n is larger than the maximum number the underlying
291       random number generator can return, then a GSL_error exception
292       is thrown.
293    */
294    DiscreteUniform(u_long n=0);
295
296    /**
297       \brief Get a random number
298
299       The returned integer is either in the range [RNG.min,RNG.max]
300       or [0,n-1] depending on how the random number generator was
301       created.
302
303       \see DiscreteUniform(const u_long n=0)
304    */
305    u_long operator()(void) const;
306
307    /**
308       \brief Get a random integer in the range \f$ [0,n-1] \f$.
309
310       All integers in the range [0,n-1] are equally likely. This
311       function should be avoided for sampling the whole range of the
312       underlying RNG.
313
314       \throw GSL_error if \an is larger than the range of the
315       underlying generator.
316    */
317    u_long operator()(u_long n) const;
318
319  private:
320    u_long range_;
321  };
322
323  ///
324  /// @brief Poisson Distribution
325  ///
326  /// Having a Poisson process (i.e. no memory), number of occurences
327  /// within a given time window is Poisson distributed. This
328  /// distribution is the limit of a Binomial distribution when number
329  /// of attempts is large, and the probability for one attempt to be
330  /// succesful is small (in such a way that the expected number of
331  /// succesful attempts is \f$ m \f$.
332  ///
333  /// Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
334  /// k  \f$ \n
335  /// Expectation value: \f$ m \f$ \n
336  /// Variance: \f$ m \f$
337  ///
338  class Poisson : public Discrete
339  {
340  public:
341    ///
342    /// @brief Constructor
343    ///
344    /// @param m is expectation value
345    ///
346    Poisson(const double m=1);
347
348    ///
349    /// @return A Poisson distributed number.
350    ///
351    u_long operator()(void) const;
352
353    ///
354    /// @return A Poisson distributed number with expectation value
355    /// \a m
356    ///
357    /// @note this operator ignores parameters set in Constructor
358    ///
359    u_long operator()(const double m) const;
360
361  private:
362    double m_;
363  };
364
365  // --------------------- Continuous distribtuions ---------------------
366
367  ///
368  /// @brief Continuous random number distributions.
369  ///
370  /// Abstract base class for continuous random number distributions.
371  ///
372  class Continuous
373  {
374  public:
375
376    ///
377    /// @brief Constructor
378    ///
379    Continuous(void);
380
381    ///
382    /// @brief The destructor
383    ///
384    virtual ~Continuous(void);
385
386    ///
387    /// @brief Set the seed to \a s.
388    ///
389    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
390    /// default value from the rng's original implementation is used
391    /// (cf. GSL documentation).
392    ///
393    /// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
394    ///
395    void seed(u_long s) const;
396
397    ///
398    /// @brief Set the seed using the /dev/urandom device.
399    ///
400    /// @return The seed acquired from /dev/urandom.
401    ///
402    /// @see seed, RNG::seed_from_devurandom, RNG::seed
403    ///
404    u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
405
406    ///
407    /// @return A random number
408    ///
409    virtual double operator()(void) const = 0;
410
411  protected:
412    /// pointer to GSL random generator
413    RNG* rng_;
414  };
415
416  // ContinuousUniform is declared before ContinuousGeneral to avoid
417  // forward declaration
418  ///
419  /// @brief Uniform distribution
420  ///
421  /// Class for generating a random number from a uniform distribution
422  /// in the range [0,1), i.e. zero is included but not 1.
423  ///
424  /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
425  /// Expectation value: 0.5 \n
426  /// Variance: \f$ \frac{1}{12} \f$
427  ///
428  class ContinuousUniform : public Continuous
429  {
430  public:
431    double operator()(void) const;
432  };
433
434  ///
435  /// @brief Generates numbers from a histogram in a continuous manner.
436  ///
437  class ContinuousGeneral : public Continuous
438  {
439  public:
440    ///
441    /// @brief Constructor
442    ///
443    /// @param hist is a Histogram defining the probability distribution
444    ///
445    ContinuousGeneral(const statistics::Histogram& hist);
446
447    ///
448    /// The number is generated in a two step process. First the bin
449    /// in the histogram is randomly selected (see
450    /// DiscreteGeneral). Then a number is generated uniformly from
451    /// the interval defined by the bin.
452    ///
453    /// @return A random number.
454    ///
455    double operator()(void) const;
456
457  private:
458    const DiscreteGeneral discrete_;
459    const statistics::Histogram hist_;
460    ContinuousUniform u_;
461  };
462
463  ///
464  /// @brief Generator of random numbers from an exponential
465  /// distribution.
466  ///
467  /// The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
468  /// \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
469  /// variance \f$ m^2 \f$
470  ///
471  class Exponential : public Continuous
472  {
473  public:
474    ///
475    /// @brief Constructor
476    ///
477    /// @param m is the expectation value of the distribution.
478    ///
479    Exponential(const double m=1);
480
481    ///
482    /// @return A random number from exponential distribution.
483    ///
484    double operator()(void) const;
485
486    ///
487    /// @return A random number from exponential distribution, with
488    /// expectation value \a m
489    ///
490    /// @note This operator ignores parameters given in constructor.
491    ///
492    double operator()(const double m) const;
493
494  private:
495    double m_;
496  };
497
498  ///
499  /// @brief Gaussian distribution
500  ///
501  /// Class for generating a random number from a Gaussian
502  /// distribution between zero and unity. Utilizes the Box-Muller
503  /// algorithm, which needs two calls to random generator.
504  ///
505  /// Distribution function \f$ f(x) =
506  /// \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
507  /// \f$ \n
508  /// Expectation value: \f$ \mu \f$ \n
509  /// Variance: \f$ \sigma^2 \f$
510  ///
511  class Gaussian : public Continuous
512  {
513  public:
514    ///
515    /// @brief Constructor
516    ///
517    /// @param s is the standard deviation \f$ \sigma \f$ of distribution
518    /// @param m is the expectation value \f$ \mu \f$ of the distribution
519    ///
520    Gaussian(const double s=1, const double m=0);
521
522    ///
523    /// @return A random Gaussian number
524    ///
525    double operator()(void) const;
526
527    ///
528    /// @return A random Gaussian number with standard deviation \a s
529    /// and expectation value 0.
530    ///
531    /// @note this operator ignores parameters given in Constructor
532    ///
533    double operator()(const double s) const;
534
535    ///
536    /// @return A random Gaussian number with standard deviation \a s
537    /// and expectation value \a m.
538    ///
539    /// @note this operator ignores parameters given in Constructor
540    ///
541    double operator()(const double s, const double m) const;
542
543  private:
544    double m_;
545    double s_;
546  };
547
548}}} // of namespace random, yat, and theplu
549
550#endif
Note: See TracBrowser for help on using the repository browser.