yat/random/random.h

Code
Comments
Other
Rev Date Author Line
680 11 Oct 06 jari 1 #ifndef _theplu_yat_random_
680 11 Oct 06 jari 2 #define _theplu_yat_random_
675 10 Oct 06 jari 3
368 05 Aug 05 peter 4 // $Id$
361 04 Aug 05 jari 5
561 15 Mar 06 jari 6 /*
2119 12 Dec 09 peter 7   Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 8   Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2020, 2022, 2023 Peter Johansson
561 15 Mar 06 jari 9
1437 25 Aug 08 peter 10   This file is part of the yat library, http://dev.thep.lu.se/yat
561 15 Mar 06 jari 11
675 10 Oct 06 jari 12   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 13   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 14   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 15   License, or (at your option) any later version.
561 15 Mar 06 jari 16
675 10 Oct 06 jari 17   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 18   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 19   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
561 15 Mar 06 jari 20   General Public License for more details.
561 15 Mar 06 jari 21
561 15 Mar 06 jari 22   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 23   along with yat. If not, see <http://www.gnu.org/licenses/>.
561 15 Mar 06 jari 24 */
561 15 Mar 06 jari 25
4281 29 Jan 23 peter 26 #include "yat/utility/CholeskyDecomposer.h"
3591 20 Jan 17 peter 27 #include "yat/utility/config_public.h"
4281 29 Jan 23 peter 28 #include "yat/utility/Vector.h"
3591 20 Jan 17 peter 29
820 17 Mar 07 peter 30 #include "yat/statistics/Histogram.h"
2578 04 Oct 11 peter 31 #include "yat/utility/deprecate.h"
820 17 Mar 07 peter 32
2145 15 Jan 10 peter 33 #include <boost/concept_check.hpp>
2145 15 Jan 10 peter 34
362 04 Aug 05 jari 35 #include <gsl/gsl_rng.h>
362 04 Aug 05 jari 36 #include <gsl/gsl_randist.h>
361 04 Aug 05 jari 37
1004 23 Jan 08 peter 38 #include <algorithm>
3591 20 Jan 17 peter 39 #include <atomic>
3959 30 Jul 20 peter 40 #include <memory>
3959 30 Jul 20 peter 41 #include <mutex>
4252 18 Nov 22 peter 42 #include <random>
374 05 Aug 05 jari 43 #include <string>
1616 05 Nov 08 peter 44 #include <vector>
374 05 Aug 05 jari 45
361 04 Aug 05 jari 46 namespace theplu {
680 11 Oct 06 jari 47 namespace yat {
4281 29 Jan 23 peter 48   namespace utility {
4281 29 Jan 23 peter 49     class MatrixBase;
4281 29 Jan 23 peter 50   }
424 07 Dec 05 jari 51 namespace random {
361 04 Aug 05 jari 52
2578 04 Oct 11 peter 53   //forward declaration
674 10 Oct 06 peter 54   class RNG_state;
674 10 Oct 06 peter 55
361 04 Aug 05 jari 56   ///
389 15 Aug 05 peter 57   /// @brief Random Number Generator
389 15 Aug 05 peter 58   ///
424 07 Dec 05 jari 59   /// The RNG class is wrapper to the GSL random number generator
2802 29 Jul 12 peter 60   /// (rng). In yat 0.8 (or older) this class provided a single global
2802 29 Jul 12 peter 61   /// instance of the rng, and made sure there was only one point of
2802 29 Jul 12 peter 62   /// access to the generator. Since version 0.9 this class provides
2802 29 Jul 12 peter 63   /// one rng per thread in order to avoid collisions in multi-thread
2802 29 Jul 12 peter 64   /// applications.
361 04 Aug 05 jari 65   ///
2802 29 Jul 12 peter 66   /// There are many different rng's available in GSL. RNG uses the
2802 29 Jul 12 peter 67   /// default generator, unless the global variable \c gsl_rng_default
2802 29 Jul 12 peter 68   /// has been modified (see <a
2802 29 Jul 12 peter 69   /// href=\gsl_url/Random-number-environment-variables.html>GSL
2802 29 Jul 12 peter 70   /// Manual</a>). Note, \c gsl_rng_default should be changed before
2802 29 Jul 12 peter 71   /// RNG creates its generator and safest way to achieve this is to
2802 29 Jul 12 peter 72   /// modify \c gsl_rng_default prior calling instance() the first
2802 29 Jul 12 peter 73   /// time.
2802 29 Jul 12 peter 74   ///
374 05 Aug 05 jari 75   /// There is information about how to change seeding and generators
738 07 Jan 07 jari 76   /// at run time without recompilation using environment variables in
2802 29 Jul 12 peter 77   /// the <a
2802 29 Jul 12 peter 78   /// href=\gsl_url/Random-number-environment-variables.html>GSL
2802 29 Jul 12 peter 79   /// Manual</a>. RNG supports seeding at compile time if you don't
2802 29 Jul 12 peter 80   /// want to bother about environment variables and GSL.
362 04 Aug 05 jari 81   ///
2802 29 Jul 12 peter 82   /// The class provides one generator per thread. The first generator
2802 29 Jul 12 peter 83   /// created is seeded with \c gsl_rng_default_seed and subsequent
2802 29 Jul 12 peter 84   /// generators are seeded with \c gsl_rng_default_seed + 1, \c
2802 29 Jul 12 peter 85   /// gsl_rng_default_seed + 2 etc, unless the seed has been modified
2802 29 Jul 12 peter 86   /// with seed() or seed_from_devurandom().
362 04 Aug 05 jari 87   ///
2802 29 Jul 12 peter 88   /// @see Design Patterns (the singleton and adapter pattern). <a
2802 29 Jul 12 peter 89   /// href=\gsl_url/Random-Number-Generation.html>GSL documentation</a>.
374 05 Aug 05 jari 90   ///
362 04 Aug 05 jari 91   class RNG
362 04 Aug 05 jari 92   {
361 04 Aug 05 jari 93   public:
361 04 Aug 05 jari 94
424 07 Dec 05 jari 95     ///
2802 29 Jul 12 peter 96     /// @brief Get an instance of the Random Number Generator.
561 15 Mar 06 jari 97     ///
2802 29 Jul 12 peter 98     /// Get an instance of RNG. If a random number generator is not
2802 29 Jul 12 peter 99     /// already created for current thread, the call will create a new
2802 29 Jul 12 peter 100     /// generator of type \c gsl_rng_default. If it is the first
2802 29 Jul 12 peter 101     /// generator created it will be seeded with \c
2802 29 Jul 12 peter 102     /// gsl_rng_default_seed; otherwise created generator will be
2802 29 Jul 12 peter 103     /// seeded with \c seed + \c n, where \c seed is the latest seed
2802 29 Jul 12 peter 104     /// set (with seed() or seed_from_devurandom()) The seed may be
561 15 Mar 06 jari 105     /// changed with the seed or seed_from_devurandom member
425 07 Dec 05 jari 106     /// functions.
424 07 Dec 05 jari 107     ///
424 07 Dec 05 jari 108     /// @return A pointer to the random number generator.
424 07 Dec 05 jari 109     ///
425 07 Dec 05 jari 110     /// @see seed and seed_from_devurandom
425 07 Dec 05 jari 111     ///
752 17 Feb 07 jari 112     static RNG* instance(void);
362 04 Aug 05 jari 113
374 05 Aug 05 jari 114     ///
374 05 Aug 05 jari 115     /// @brief Returns the largest number that the random number
374 05 Aug 05 jari 116     /// generator can return.
374 05 Aug 05 jari 117     ///
1271 09 Apr 08 peter 118     unsigned long max(void) const;
374 05 Aug 05 jari 119
374 05 Aug 05 jari 120     ///
374 05 Aug 05 jari 121     /// @brief Returns the smallest number that the random number
374 05 Aug 05 jari 122     /// generator can return.
374 05 Aug 05 jari 123     ///
1271 09 Apr 08 peter 124     unsigned long min(void) const;
374 05 Aug 05 jari 125
374 05 Aug 05 jari 126     ///
374 05 Aug 05 jari 127     /// @brief Returns the name of the random number generator
374 05 Aug 05 jari 128     ///
718 26 Dec 06 jari 129     std::string name(void) const;
374 05 Aug 05 jari 130
648 14 Sep 06 peter 131     ///
2802 29 Jul 12 peter 132     /// Access underlying GSL random number generator speicific to
2802 29 Jul 12 peter 133     /// current thread. Behaviour of returned generator is undefined
2802 29 Jul 12 peter 134     /// outside current thread.
2802 29 Jul 12 peter 135     ///
648 14 Sep 06 peter 136     /// @return const pointer to underlying GSL random generator.
648 14 Sep 06 peter 137     ///
718 26 Dec 06 jari 138     const gsl_rng* rng(void) const;
362 04 Aug 05 jari 139
374 05 Aug 05 jari 140     ///
561 15 Mar 06 jari 141     /// @brief Set the seed \a s for the rng.
561 15 Mar 06 jari 142     ///
374 05 Aug 05 jari 143     /// Set the seed \a s for the rng. If \a s is zero, a default
374 05 Aug 05 jari 144     /// value from the rng's original implementation is used (cf. GSL
374 05 Aug 05 jari 145     /// documentation).
374 05 Aug 05 jari 146     ///
2802 29 Jul 12 peter 147     /// This function will also effect generators created subsequently
2802 29 Jul 12 peter 148     /// in other threads. The seed \a s is saved and subsequent
2802 29 Jul 12 peter 149     /// generators will be created with seed \c s + 1, \c s + 2, etc.
2802 29 Jul 12 peter 150     ///
424 07 Dec 05 jari 151     /// @see seed_from_devurandom
374 05 Aug 05 jari 152     ///
1271 09 Apr 08 peter 153     void seed(unsigned long s) const;
374 05 Aug 05 jari 154
374 05 Aug 05 jari 155     ///
374 05 Aug 05 jari 156     /// @brief Seed the rng using the /dev/urandom device.
374 05 Aug 05 jari 157     ///
2802 29 Jul 12 peter 158     /// This function will also effect generators in other threads
2802 29 Jul 12 peter 159     /// created subsequntly (see seed()).
2802 29 Jul 12 peter 160     ///
374 05 Aug 05 jari 161     /// @return The seed acquired from /dev/urandom.
374 05 Aug 05 jari 162     ///
1271 09 Apr 08 peter 163     unsigned long seed_from_devurandom(void);
374 05 Aug 05 jari 164
738 07 Jan 07 jari 165     /**
738 07 Jan 07 jari 166        \brief Set the state to \a state.
738 07 Jan 07 jari 167
2804 31 Jul 12 peter 168        \return 0 always.
738 07 Jan 07 jari 169
2802 29 Jul 12 peter 170        \note this function only effects the RNG in current thread
2802 29 Jul 12 peter 171
2804 31 Jul 12 peter 172        \throw utility::GSL_error on error
2804 31 Jul 12 peter 173
738 07 Jan 07 jari 174        \see gsl_rng_memcpy
738 07 Jan 07 jari 175     */
2804 31 Jul 12 peter 176     // return int for backward compatibility with yat 0.8
674 10 Oct 06 peter 177     int set_state(const RNG_state&);
674 10 Oct 06 peter 178
361 04 Aug 05 jari 179   private:
425 07 Dec 05 jari 180     RNG(void);
362 04 Aug 05 jari 181
852 01 May 07 jari 182     /**
852 01 May 07 jari 183        \brief Not implemented.
852 01 May 07 jari 184
852 01 May 07 jari 185        This copy contructor is not implemented. The constructor is
852 01 May 07 jari 186        declared in order to avoid compiler generated default copy
852 01 May 07 jari 187        constructor.
852 01 May 07 jari 188      */
852 01 May 07 jari 189      RNG(const RNG&);
852 01 May 07 jari 190
1614 05 Nov 08 peter 191     /**
1614 05 Nov 08 peter 192        There can be only one RNG so assignment is always
1614 05 Nov 08 peter 193        self-assignment and we do not allow it
1614 05 Nov 08 peter 194     */
1614 05 Nov 08 peter 195     RNG& operator=(const RNG&);
1614 05 Nov 08 peter 196
2770 11 Jul 12 peter 197     virtual ~RNG(void);
2770 11 Jul 12 peter 198     void rng_alloc(void) const;
2770 11 Jul 12 peter 199
3591 20 Jan 17 peter 200     static std::atomic<RNG*> instance_;
2800 28 Jul 12 peter 201     // holds one gsl_rng per thread. Access through rng(void) so a
2800 28 Jul 12 peter 202     // gsl_rng is allocated if necessary.
3959 30 Jul 20 peter 203     static thread_local std::unique_ptr<gsl_rng, void(*)(gsl_rng*)> rng_;
2802 29 Jul 12 peter 204     mutable unsigned long seed_;
2838 16 Sep 12 peter 205     // guard needs to be mutable because major mission for it is to protect seed_ against multi-access, and seed_ is mutable...
3959 30 Jul 20 peter 206     mutable std::mutex mutex_;
3959 30 Jul 20 peter 207     static std::mutex init_mutex_;
361 04 Aug 05 jari 208   };
361 04 Aug 05 jari 209
674 10 Oct 06 peter 210
674 10 Oct 06 peter 211   ///
674 10 Oct 06 peter 212   /// @brief Class holding state of a random generator
674 10 Oct 06 peter 213   ///
674 10 Oct 06 peter 214   class RNG_state
674 10 Oct 06 peter 215   {
674 10 Oct 06 peter 216   public:
674 10 Oct 06 peter 217     ///
674 10 Oct 06 peter 218     /// @brief Constructor
674 10 Oct 06 peter 219     ///
2089 21 Oct 09 peter 220     explicit RNG_state(const RNG*);
674 10 Oct 06 peter 221
1614 05 Nov 08 peter 222     /**
1614 05 Nov 08 peter 223        Copy Constructor
1614 05 Nov 08 peter 224
1614 05 Nov 08 peter 225        \since Explicitely declared since yat 0.5
1614 05 Nov 08 peter 226      */
1614 05 Nov 08 peter 227     RNG_state(const RNG_state&);
1614 05 Nov 08 peter 228
674 10 Oct 06 peter 229     ///
674 10 Oct 06 peter 230     /// @brief Destructor
674 10 Oct 06 peter 231     ///
674 10 Oct 06 peter 232     ~RNG_state(void);
674 10 Oct 06 peter 233
674 10 Oct 06 peter 234     ///
674 10 Oct 06 peter 235     /// @return const pointer to underlying GSL random generator.
674 10 Oct 06 peter 236     ///
718 26 Dec 06 jari 237     const gsl_rng* rng(void) const;
674 10 Oct 06 peter 238
1614 05 Nov 08 peter 239     /**
1614 05 Nov 08 peter 240        Assignment operator
1614 05 Nov 08 peter 241
1614 05 Nov 08 peter 242        \since Explicitely declared since yat 0.5
1614 05 Nov 08 peter 243      */
1614 05 Nov 08 peter 244     RNG_state& operator=(const RNG_state&);
1614 05 Nov 08 peter 245
674 10 Oct 06 peter 246   private:
674 10 Oct 06 peter 247     gsl_rng* rng_;
674 10 Oct 06 peter 248
1614 05 Nov 08 peter 249     void clone(const gsl_rng&);
674 10 Oct 06 peter 250   };
674 10 Oct 06 peter 251
2889 07 Dec 12 peter 252
424 07 Dec 05 jari 253   // --------------------- Discrete distribtuions ---------------------
361 04 Aug 05 jari 254
366 05 Aug 05 peter 255   ///
424 07 Dec 05 jari 256   /// @brief Discrete random number distributions.
366 05 Aug 05 peter 257   ///
424 07 Dec 05 jari 258   /// Abstract base class for discrete random number
424 07 Dec 05 jari 259   /// distributions. Given K discrete events with different
424 07 Dec 05 jari 260   /// probabilities \f$ P[k] \f$, produce a random value k consistent
424 07 Dec 05 jari 261   /// with its probability.
367 05 Aug 05 peter 262   ///
425 07 Dec 05 jari 263   class Discrete
362 04 Aug 05 jari 264   {
362 04 Aug 05 jari 265   public:
2901 13 Dec 12 peter 266     /**
2901 13 Dec 12 peter 267        type returned by operator()
2901 13 Dec 12 peter 268
2901 13 Dec 12 peter 269        \since New in yat 0.10
2901 13 Dec 12 peter 270      */
2901 13 Dec 12 peter 271     typedef unsigned long int result_type;
2901 13 Dec 12 peter 272
366 05 Aug 05 peter 273     ///
366 05 Aug 05 peter 274     /// @brief Constructor
366 05 Aug 05 peter 275     ///
706 19 Dec 06 jari 276     Discrete(void);
362 04 Aug 05 jari 277
706 19 Dec 06 jari 278     ///
706 19 Dec 06 jari 279     /// @brief The destructor
706 19 Dec 06 jari 280     ///
706 19 Dec 06 jari 281     virtual ~Discrete(void);
443 15 Dec 05 jari 282
367 05 Aug 05 peter 283     ///
561 15 Mar 06 jari 284     /// @brief Set the seed to \a s.
561 15 Mar 06 jari 285     ///
425 07 Dec 05 jari 286     /// Set the seed to \a s in the underlying rng. If \a s is zero, a
425 07 Dec 05 jari 287     /// default value from the rng's original implementation is used
425 07 Dec 05 jari 288     /// (cf. GSL documentation).
425 07 Dec 05 jari 289     ///
2578 04 Oct 11 peter 290     /// \deprecated Provided for backward compatibility with the 0.7
2578 04 Oct 11 peter 291     /// API. Use RNG::instance()->seed(s) instead.
425 07 Dec 05 jari 292     ///
2578 04 Oct 11 peter 293     void seed(unsigned long s) const YAT_DEPRECATE;
425 07 Dec 05 jari 294
425 07 Dec 05 jari 295     ///
425 07 Dec 05 jari 296     /// @brief Set the seed using the /dev/urandom device.
425 07 Dec 05 jari 297     ///
425 07 Dec 05 jari 298     /// @return The seed acquired from /dev/urandom.
425 07 Dec 05 jari 299     ///
2578 04 Oct 11 peter 300     /// \deprecated Provided for backward compatibility with the 0.7
2580 05 Oct 11 peter 301     /// API. Use RNG::instance()->seed_from_devurandom() instead.
425 07 Dec 05 jari 302     ///
2578 04 Oct 11 peter 303     unsigned long seed_from_devurandom(void) YAT_DEPRECATE;
425 07 Dec 05 jari 304
425 07 Dec 05 jari 305     ///
424 07 Dec 05 jari 306     /// @return A random number.
367 05 Aug 05 peter 307     ///
2901 13 Dec 12 peter 308     virtual result_type operator()(void) const = 0;
2889 07 Dec 12 peter 309
362 04 Aug 05 jari 310   protected:
648 14 Sep 06 peter 311     /// GSL random gererator
362 04 Aug 05 jari 312     RNG* rng_;
361 04 Aug 05 jari 313   };
361 04 Aug 05 jari 314
2889 07 Dec 12 peter 315   /**
2889 07 Dec 12 peter 316      \brief Binomial distribution
2889 07 Dec 12 peter 317
2889 07 Dec 12 peter 318      \see gsl_ran_binomial
2889 07 Dec 12 peter 319
2889 07 Dec 12 peter 320      \since New in yat 0.10
2889 07 Dec 12 peter 321    */
3465 10 Feb 16 peter 322   class Binomial : public Discrete
2889 07 Dec 12 peter 323   {
2889 07 Dec 12 peter 324   public:
2889 07 Dec 12 peter 325     /**
2889 07 Dec 12 peter 326        \brief Constructor
2889 07 Dec 12 peter 327
2889 07 Dec 12 peter 328        Create an object that generates random numbers from binomial
2889 07 Dec 12 peter 329        distribution, the number of successes in \n trials each with
2889 07 Dec 12 peter 330        success probability \a p.
2889 07 Dec 12 peter 331      */
2889 07 Dec 12 peter 332     Binomial(double p, unsigned int n);
2889 07 Dec 12 peter 333
2889 07 Dec 12 peter 334     /**
2889 07 Dec 12 peter 335        \return number from binomial distrubtion as parametrized in constructor.
2889 07 Dec 12 peter 336      */
2889 07 Dec 12 peter 337     unsigned long operator()(void) const;
2889 07 Dec 12 peter 338   private:
2889 07 Dec 12 peter 339     double p_;
2889 07 Dec 12 peter 340     unsigned int n_;
2889 07 Dec 12 peter 341   };
2889 07 Dec 12 peter 342
366 05 Aug 05 peter 343   ///
424 07 Dec 05 jari 344   /// @brief General
366 05 Aug 05 peter 345   ///
2889 07 Dec 12 peter 346   class DiscreteGeneral : public Discrete
362 04 Aug 05 jari 347   {
362 04 Aug 05 jari 348   public:
4009 19 Oct 20 peter 349     /**
4009 19 Oct 20 peter 350        Default constructor leaves the object uninitialized and
4009 19 Oct 20 peter 351        behaviour of operator() is undefined.
4009 19 Oct 20 peter 352
4009 19 Oct 20 peter 353        \since New in yat 0.19
4009 19 Oct 20 peter 354      */
4009 19 Oct 20 peter 355     DiscreteGeneral(void);
4009 19 Oct 20 peter 356
367 05 Aug 05 peter 357     ///
367 05 Aug 05 peter 358     /// @brief Constructor
367 05 Aug 05 peter 359     ///
424 07 Dec 05 jari 360     /// @param hist is a Histogram defining the probability distribution
367 05 Aug 05 peter 361     ///
2089 21 Oct 09 peter 362     explicit DiscreteGeneral(const statistics::Histogram& hist);
2889 07 Dec 12 peter 363
1616 05 Nov 08 peter 364     /**
4384 13 Oct 23 peter 365        Construct a DiscreteGeneral with probability distribution such
4384 13 Oct 23 peter 366        that P(0) is proportional to *first, P(1) is proportional to
4384 13 Oct 23 peter 367        *(++first) etc. The distribution is normalized such that \f$
4384 13 Oct 23 peter 368        \sum P = 1.0 \f$.
4384 13 Oct 23 peter 369
4384 13 Oct 23 peter 370        \since New in yat 0.22
4384 13 Oct 23 peter 371      */
4384 13 Oct 23 peter 372     template<typename Iterator>
4384 13 Oct 23 peter 373     DiscreteGeneral(Iterator first, Iterator last)
4384 13 Oct 23 peter 374       : p_(first, last)
4384 13 Oct 23 peter 375     {
4384 13 Oct 23 peter 376       preproc();
4384 13 Oct 23 peter 377     }
4384 13 Oct 23 peter 378
4384 13 Oct 23 peter 379     /**
1616 05 Nov 08 peter 380        \brief Copy constructor
1616 05 Nov 08 peter 381
1616 05 Nov 08 peter 382        \since Explicitely implemented in yat 0.5
1616 05 Nov 08 peter 383      */
1616 05 Nov 08 peter 384     DiscreteGeneral(const DiscreteGeneral&);
1616 05 Nov 08 peter 385
4010 19 Oct 20 peter 386     /**
4010 19 Oct 20 peter 387        \brief Move constructor
4010 19 Oct 20 peter 388
4010 19 Oct 20 peter 389        \since New in yat 0.19
4010 19 Oct 20 peter 390      */
4010 19 Oct 20 peter 391     DiscreteGeneral(DiscreteGeneral&&);
4010 19 Oct 20 peter 392
367 05 Aug 05 peter 393     ///
424 07 Dec 05 jari 394     /// @brief Destructor
367 05 Aug 05 peter 395     ///
443 15 Dec 05 jari 396     ~DiscreteGeneral(void);
362 04 Aug 05 jari 397
367 05 Aug 05 peter 398     ///
2747 18 Jun 12 peter 399     /// The generated number is an integer and proportional to the
424 07 Dec 05 jari 400     /// frequency in the corresponding histogram bin. In other words,
424 07 Dec 05 jari 401     /// the probability that 0 is returned is proportinal to the size
424 07 Dec 05 jari 402     /// of the first bin.
367 05 Aug 05 peter 403     ///
424 07 Dec 05 jari 404     /// @return A random number.
364 05 Aug 05 peter 405     ///
1271 09 Apr 08 peter 406     unsigned long operator()(void) const;
362 04 Aug 05 jari 407
1616 05 Nov 08 peter 408     /**
1616 05 Nov 08 peter 409        \brief Assignment operator
1616 05 Nov 08 peter 410
1616 05 Nov 08 peter 411        \since Explicitely implemented in yat 0.5
1616 05 Nov 08 peter 412      */
1616 05 Nov 08 peter 413     DiscreteGeneral& operator=(const DiscreteGeneral&);
1616 05 Nov 08 peter 414
4010 19 Oct 20 peter 415     /**
4010 19 Oct 20 peter 416        \brief Move assignment operator
4010 19 Oct 20 peter 417
4010 19 Oct 20 peter 418        \since New in yat 0.19
4010 19 Oct 20 peter 419      */
4010 19 Oct 20 peter 420     DiscreteGeneral& operator=(DiscreteGeneral&&);
4010 19 Oct 20 peter 421
364 05 Aug 05 peter 422   private:
1616 05 Nov 08 peter 423     void free(void);
1616 05 Nov 08 peter 424     void preproc(void);
1616 05 Nov 08 peter 425
1616 05 Nov 08 peter 426     gsl_ran_discrete_t* gen_;
1616 05 Nov 08 peter 427     std::vector<double> p_;
4010 19 Oct 20 peter 428     friend void swap(DiscreteGeneral& lhs, DiscreteGeneral& rhs);
361 04 Aug 05 jari 429   };
361 04 Aug 05 jari 430
1342 18 Jun 08 peter 431   /**
4010 19 Oct 20 peter 432      \relates DiscreteGeneral
4010 19 Oct 20 peter 433
4010 19 Oct 20 peter 434      \since Specialisation new in yat 0.19
4010 19 Oct 20 peter 435    */
4010 19 Oct 20 peter 436   void swap(DiscreteGeneral& lhs, DiscreteGeneral& rhs);
4010 19 Oct 20 peter 437
4010 19 Oct 20 peter 438   /**
1342 18 Jun 08 peter 439      @brief Discrete uniform distribution
2889 07 Dec 12 peter 440
1342 18 Jun 08 peter 441      Discrete uniform distribution also known as the "equally likely
1342 18 Jun 08 peter 442      outcomes" distribution. Each outcome, in this case an integer
1342 18 Jun 08 peter 443      from [0,n-1] , have equal probability to occur.
2889 07 Dec 12 peter 444
1342 18 Jun 08 peter 445      Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le
2889 07 Dec 12 peter 446      k < n \f$ \n
2889 07 Dec 12 peter 447      Expectation value: \f$ \frac{n-1}{2} \f$ \n
1342 18 Jun 08 peter 448      Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$
1342 18 Jun 08 peter 449   */
3893 27 Mar 20 peter 450   class DiscreteUniform : public Discrete
364 05 Aug 05 peter 451   {
364 05 Aug 05 peter 452   public:
3893 27 Mar 20 peter 453     /// argument type is unsigned long int
3893 27 Mar 20 peter 454     typedef unsigned long argument_type;
3893 27 Mar 20 peter 455
738 07 Jan 07 jari 456     /**
738 07 Jan 07 jari 457        \brief Constructor.
377 07 Aug 05 jari 458
738 07 Jan 07 jari 459        The generator will generate integers within the range \f$
738 07 Jan 07 jari 460        [0,n-1] \f$. If \a n is zero, then the whole range of the
738 07 Jan 07 jari 461        underlying RNG will be used \f$ [min,max] \f$. Setting \a n to
738 07 Jan 07 jari 462        zero is the preferred way to sample the whole range of the
738 07 Jan 07 jari 463        underlying RNG, i.e. not setting \n to RNG.max.
424 07 Dec 05 jari 464
738 07 Jan 07 jari 465        \throw If \a n is larger than the maximum number the underlying
738 07 Jan 07 jari 466        random number generator can return, then a GSL_error exception
738 07 Jan 07 jari 467        is thrown.
738 07 Jan 07 jari 468     */
2089 21 Oct 09 peter 469     explicit DiscreteUniform(unsigned long n=0);
738 07 Jan 07 jari 470
738 07 Jan 07 jari 471     /**
738 07 Jan 07 jari 472        \brief Get a random number
738 07 Jan 07 jari 473
738 07 Jan 07 jari 474        The returned integer is either in the range [RNG.min,RNG.max]
738 07 Jan 07 jari 475        or [0,n-1] depending on how the random number generator was
738 07 Jan 07 jari 476        created.
738 07 Jan 07 jari 477
1271 09 Apr 08 peter 478        \see DiscreteUniform(const unsigned long n=0)
738 07 Jan 07 jari 479     */
1271 09 Apr 08 peter 480     unsigned long operator()(void) const;
364 05 Aug 05 peter 481
738 07 Jan 07 jari 482     /**
738 07 Jan 07 jari 483        \brief Get a random integer in the range \f$ [0,n-1] \f$.
364 05 Aug 05 peter 484
738 07 Jan 07 jari 485        All integers in the range [0,n-1] are equally likely. This
738 07 Jan 07 jari 486        function should be avoided for sampling the whole range of the
738 07 Jan 07 jari 487        underlying RNG.
738 07 Jan 07 jari 488
966 11 Oct 07 peter 489        \throw GSL_error if \a n is larger than the range of the
738 07 Jan 07 jari 490        underlying generator.
738 07 Jan 07 jari 491     */
1271 09 Apr 08 peter 492     unsigned long operator()(unsigned long n) const;
738 07 Jan 07 jari 493
3894 27 Mar 20 peter 494     /**
3894 27 Mar 20 peter 495        For n>0 the min value is 0. For n=0 (default) the min value
3894 27 Mar 20 peter 496        is RNG::min(), typically zero but depending on generator in use.
3894 27 Mar 20 peter 497
3894 27 Mar 20 peter 498        \return smallest possible value
3894 27 Mar 20 peter 499
3894 27 Mar 20 peter 500        \since New in yat 0.18
3894 27 Mar 20 peter 501      */
3894 27 Mar 20 peter 502     unsigned long int min(void) const;
3894 27 Mar 20 peter 503
3894 27 Mar 20 peter 504     /**
3894 27 Mar 20 peter 505        For n>0 the max value is n-1. For n=0 (default) the max value
3894 27 Mar 20 peter 506        is RNG::max().
3894 27 Mar 20 peter 507
3894 27 Mar 20 peter 508        \return maximal value that class returns
3894 27 Mar 20 peter 509
3894 27 Mar 20 peter 510        \since New in yat 0.18
3894 27 Mar 20 peter 511      */
3894 27 Mar 20 peter 512     unsigned long int max(void) const;
3894 27 Mar 20 peter 513
364 05 Aug 05 peter 514   private:
1271 09 Apr 08 peter 515     unsigned long range_;
364 05 Aug 05 peter 516   };
364 05 Aug 05 peter 517
3446 02 Dec 15 peter 518
1342 18 Jun 08 peter 519   /**
3446 02 Dec 15 peter 520      @brief Geomtric Distribution
3446 02 Dec 15 peter 521
3446 02 Dec 15 peter 522      The number of independent trials with probability \em p until the
3446 02 Dec 15 peter 523      first success.
3446 02 Dec 15 peter 524
3446 02 Dec 15 peter 525      Probability function \f$ p(k) = p (1-p)^(k-1) \f$
3446 02 Dec 15 peter 526
3446 02 Dec 15 peter 527      \since New in yat 0.14
3446 02 Dec 15 peter 528   */
3446 02 Dec 15 peter 529   class Geometric : public Discrete
3446 02 Dec 15 peter 530   {
3446 02 Dec 15 peter 531   public:
3446 02 Dec 15 peter 532     /**
3446 02 Dec 15 peter 533        \brief Constructor
3446 02 Dec 15 peter 534
3446 02 Dec 15 peter 535        \param p is probability for success in one trial
3446 02 Dec 15 peter 536     */
3446 02 Dec 15 peter 537     Geometric(double p);
3446 02 Dec 15 peter 538
3446 02 Dec 15 peter 539     /*
3446 02 Dec 15 peter 540       \return a number from Geomtric distribution
3446 02 Dec 15 peter 541     */
3446 02 Dec 15 peter 542     unsigned long operator()(void) const;
3446 02 Dec 15 peter 543
3446 02 Dec 15 peter 544     /**
3446 02 Dec 15 peter 545        \return a number from Geomtric distribution with success rate \a p
3446 02 Dec 15 peter 546
3446 02 Dec 15 peter 547        \note this operator ignores parameters set in Constructor
3446 02 Dec 15 peter 548     */
3446 02 Dec 15 peter 549     unsigned long operator()(double p) const;
3446 02 Dec 15 peter 550
3446 02 Dec 15 peter 551   private:
3446 02 Dec 15 peter 552     double p_;
3446 02 Dec 15 peter 553   };
3446 02 Dec 15 peter 554
3446 02 Dec 15 peter 555
3446 02 Dec 15 peter 556   /**
3465 10 Feb 16 peter 557      If we have \a n1 samples of type 1 and \a n2 samples of type 2
3465 10 Feb 16 peter 558      and draw \a t samples with replacement, number of drawn samples
3465 10 Feb 16 peter 559      of type 1 will follow the hyper geometric distribution.
3465 10 Feb 16 peter 560
3465 10 Feb 16 peter 561      \since New in yat 0.14
3465 10 Feb 16 peter 562    */
3465 10 Feb 16 peter 563   class HyperGeometric : public Discrete
3465 10 Feb 16 peter 564   {
3465 10 Feb 16 peter 565   public:
3465 10 Feb 16 peter 566     /**
3465 10 Feb 16 peter 567        \brief Defaul constructor
3465 10 Feb 16 peter 568      */
3465 10 Feb 16 peter 569     HyperGeometric(void);
3465 10 Feb 16 peter 570
3465 10 Feb 16 peter 571     /**
3465 10 Feb 16 peter 572        \brief Constructor
3465 10 Feb 16 peter 573        \param n1 number of samples of type 1
3465 10 Feb 16 peter 574        \param n2 number of samples of type 2
3465 10 Feb 16 peter 575        \param t number of samples to draw
3465 10 Feb 16 peter 576      */
3465 10 Feb 16 peter 577     HyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);
3465 10 Feb 16 peter 578
3465 10 Feb 16 peter 579     /**
3465 10 Feb 16 peter 580        \return random number from hypergeometric distribution using
3465 10 Feb 16 peter 581        parameters set in constructor.
3465 10 Feb 16 peter 582      */
3465 10 Feb 16 peter 583     unsigned long int operator()(void) const;
3465 10 Feb 16 peter 584
3465 10 Feb 16 peter 585     /**
3465 10 Feb 16 peter 586        \return random number from hypergeometric distribution using
3465 10 Feb 16 peter 587        parameters passed.
3465 10 Feb 16 peter 588      */
3465 10 Feb 16 peter 589     unsigned long int operator()(unsigned int n1, unsigned int n2,
3465 10 Feb 16 peter 590                                  unsigned int t) const;
3465 10 Feb 16 peter 591   private:
3465 10 Feb 16 peter 592     unsigned int n1_;
3465 10 Feb 16 peter 593     unsigned int n2_;
3465 10 Feb 16 peter 594     unsigned int t_;
3465 10 Feb 16 peter 595   };
3465 10 Feb 16 peter 596
3465 10 Feb 16 peter 597   /**
3469 29 Feb 16 peter 598      We have \a n1 samples of type 1 and \a n2 samples of type 2.
3469 29 Feb 16 peter 599      Samples are drawn with replacement until \a t samles of type 2
3469 29 Feb 16 peter 600      are drawn. Then \a k, number of drawn samples of type 1, follows
3469 29 Feb 16 peter 601      negative hypergeometric distribution.
3469 29 Feb 16 peter 602
3469 29 Feb 16 peter 603      \since New in yat 0.14
3469 29 Feb 16 peter 604
3469 29 Feb 16 peter 605      \see HyperGeometric
3469 29 Feb 16 peter 606    */
3469 29 Feb 16 peter 607   class NegativeHyperGeometric : public Discrete
3469 29 Feb 16 peter 608   {
3469 29 Feb 16 peter 609   public:
3469 29 Feb 16 peter 610     /**
3469 29 Feb 16 peter 611        \brief Default constructor
3469 29 Feb 16 peter 612      */
3469 29 Feb 16 peter 613     NegativeHyperGeometric(void);
3469 29 Feb 16 peter 614
3469 29 Feb 16 peter 615     /**
3469 29 Feb 16 peter 616        \brief Constructor
3469 29 Feb 16 peter 617        \param n1 number of samples of type 1
3469 29 Feb 16 peter 618        \param n2 number of samples of type 2
3469 29 Feb 16 peter 619        \param t number of samples of type 2 to draw
3469 29 Feb 16 peter 620      */
3469 29 Feb 16 peter 621     NegativeHyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);
3469 29 Feb 16 peter 622
3469 29 Feb 16 peter 623     /**
3469 29 Feb 16 peter 624        \return random number from negative hypergeometric distribution
3469 29 Feb 16 peter 625        using parameters set in constructor.
3469 29 Feb 16 peter 626      */
3469 29 Feb 16 peter 627     unsigned long int operator()(void) const;
3469 29 Feb 16 peter 628
3469 29 Feb 16 peter 629     /**
3469 29 Feb 16 peter 630        \return random number from negative hypergeometric distribution
3469 29 Feb 16 peter 631        using parameters passed.
3469 29 Feb 16 peter 632      */
3469 29 Feb 16 peter 633     unsigned long int operator()(unsigned int n1, unsigned int n2,
3469 29 Feb 16 peter 634                                  unsigned int t) const;
3469 29 Feb 16 peter 635   private:
3469 29 Feb 16 peter 636     unsigned int n1_;
3469 29 Feb 16 peter 637     unsigned int n2_;
3469 29 Feb 16 peter 638     unsigned int t_;
3469 29 Feb 16 peter 639   };
3469 29 Feb 16 peter 640
3469 29 Feb 16 peter 641
3469 29 Feb 16 peter 642   /**
1342 18 Jun 08 peter 643      @brief Poisson Distribution
2889 07 Dec 12 peter 644
1342 18 Jun 08 peter 645      Having a Poisson process (i.e. no memory), number of occurences
1342 18 Jun 08 peter 646      within a given time window is Poisson distributed. This
1342 18 Jun 08 peter 647      distribution is the limit of a Binomial distribution when number
1342 18 Jun 08 peter 648      of attempts is large, and the probability for one attempt to be
1342 18 Jun 08 peter 649      succesful is small (in such a way that the expected number of
1342 18 Jun 08 peter 650      succesful attempts is \f$ m \f$.
2889 07 Dec 12 peter 651
1342 18 Jun 08 peter 652      Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
2889 07 Dec 12 peter 653      k  \f$ \n
2889 07 Dec 12 peter 654      Expectation value: \f$ m \f$ \n
1342 18 Jun 08 peter 655      Variance: \f$ m \f$
1342 18 Jun 08 peter 656   */
368 05 Aug 05 peter 657   class Poisson : public Discrete
364 05 Aug 05 peter 658   {
364 05 Aug 05 peter 659   public:
364 05 Aug 05 peter 660     ///
367 05 Aug 05 peter 661     /// @brief Constructor
364 05 Aug 05 peter 662     ///
367 05 Aug 05 peter 663     /// @param m is expectation value
706 19 Dec 06 jari 664     ///
2089 21 Oct 09 peter 665     explicit Poisson(const double m=1);
364 05 Aug 05 peter 666
364 05 Aug 05 peter 667     ///
368 05 Aug 05 peter 668     /// @return A Poisson distributed number.
364 05 Aug 05 peter 669     ///
1271 09 Apr 08 peter 670     unsigned long operator()(void) const;
364 05 Aug 05 peter 671
364 05 Aug 05 peter 672     ///
718 26 Dec 06 jari 673     /// @return A Poisson distributed number with expectation value
718 26 Dec 06 jari 674     /// \a m
424 07 Dec 05 jari 675     ///
367 05 Aug 05 peter 676     /// @note this operator ignores parameters set in Constructor
424 07 Dec 05 jari 677     ///
1271 09 Apr 08 peter 678     unsigned long operator()(const double m) const;
364 05 Aug 05 peter 679
364 05 Aug 05 peter 680   private:
364 05 Aug 05 peter 681     double m_;
364 05 Aug 05 peter 682   };
364 05 Aug 05 peter 683
4147 28 Feb 22 peter 684   /**
4147 28 Feb 22 peter 685      @brief Negative Binomial Distribution
4147 28 Feb 22 peter 686
4147 28 Feb 22 peter 687      Probability function \f$ p(k) =
4148 01 Mar 22 peter 688      \frac{\Gamma(r+k)}{\Gamma(k+1)\Gamma(r)} p^r (1-p)^k\f$
4147 28 Feb 22 peter 689
4148 01 Mar 22 peter 690      Mean: \f$ m = (1-p)*r/p \f$
4148 01 Mar 22 peter 691      Variance: \f$ V = (1-p) (r-1)/p \f$
4147 28 Feb 22 peter 692
4148 01 Mar 22 peter 693      \note The parameterisation used here is different from wikipedia.
4148 01 Mar 22 peter 694
4147 28 Feb 22 peter 695      \since New in yat 0.20
4147 28 Feb 22 peter 696   */
4147 28 Feb 22 peter 697   class NegativeBinomial : public Discrete
4147 28 Feb 22 peter 698   {
4147 28 Feb 22 peter 699   public:
4147 28 Feb 22 peter 700     /**
4147 28 Feb 22 peter 701        \brief Constructor
4147 28 Feb 22 peter 702
4147 28 Feb 22 peter 703        \param p probability of success [0-1]
4147 28 Feb 22 peter 704        \param r number of successes, r>=1
4147 28 Feb 22 peter 705     */
4147 28 Feb 22 peter 706     NegativeBinomial(double p, double r);
4147 28 Feb 22 peter 707
4147 28 Feb 22 peter 708     /**
4147 28 Feb 22 peter 709        \return A number from a Negative Binomial Distribution
4147 28 Feb 22 peter 710     */
4147 28 Feb 22 peter 711     unsigned long operator()(void) const;
4147 28 Feb 22 peter 712
4147 28 Feb 22 peter 713     /**
4147 28 Feb 22 peter 714         \return A number from a Negative Binomial Distribution
4147 28 Feb 22 peter 715     */
4147 28 Feb 22 peter 716     unsigned long operator()(double p, double r) const;
4147 28 Feb 22 peter 717
4147 28 Feb 22 peter 718   private:
4147 28 Feb 22 peter 719     double p_;
4147 28 Feb 22 peter 720     double r_;
4147 28 Feb 22 peter 721   };
4147 28 Feb 22 peter 722
424 07 Dec 05 jari 723   // --------------------- Continuous distribtuions ---------------------
424 07 Dec 05 jari 724
365 05 Aug 05 peter 725   ///
424 07 Dec 05 jari 726   /// @brief Continuous random number distributions.
365 05 Aug 05 peter 727   ///
424 07 Dec 05 jari 728   /// Abstract base class for continuous random number distributions.
424 07 Dec 05 jari 729   ///
424 07 Dec 05 jari 730   class Continuous
371 05 Aug 05 peter 731   {
371 05 Aug 05 peter 732   public:
2901 13 Dec 12 peter 733     /**
2901 13 Dec 12 peter 734        type returned by operator()
424 07 Dec 05 jari 735
2901 13 Dec 12 peter 736        \since New in yat 0.10
2901 13 Dec 12 peter 737      */
2901 13 Dec 12 peter 738     typedef double result_type;
2901 13 Dec 12 peter 739
371 05 Aug 05 peter 740     ///
371 05 Aug 05 peter 741     /// @brief Constructor
371 05 Aug 05 peter 742     ///
706 19 Dec 06 jari 743     Continuous(void);
371 05 Aug 05 peter 744
706 19 Dec 06 jari 745     ///
706 19 Dec 06 jari 746     /// @brief The destructor
706 19 Dec 06 jari 747     ///
706 19 Dec 06 jari 748     virtual ~Continuous(void);
443 15 Dec 05 jari 749
371 05 Aug 05 peter 750     ///
561 15 Mar 06 jari 751     /// @brief Set the seed to \a s.
561 15 Mar 06 jari 752     ///
425 07 Dec 05 jari 753     /// Set the seed to \a s in the underlying rng. If \a s is zero, a
425 07 Dec 05 jari 754     /// default value from the rng's original implementation is used
425 07 Dec 05 jari 755     /// (cf. GSL documentation).
425 07 Dec 05 jari 756     ///
2578 04 Oct 11 peter 757     /// \deprecated Provided for backward compatibility with the 0.7
2578 04 Oct 11 peter 758     /// API. Use RNG::instance()->seed(s) instead.
425 07 Dec 05 jari 759     ///
2578 04 Oct 11 peter 760     void seed(unsigned long s) const YAT_DEPRECATE;
425 07 Dec 05 jari 761
425 07 Dec 05 jari 762     ///
425 07 Dec 05 jari 763     /// @brief Set the seed using the /dev/urandom device.
425 07 Dec 05 jari 764     ///
425 07 Dec 05 jari 765     /// @return The seed acquired from /dev/urandom.
425 07 Dec 05 jari 766     ///
2578 04 Oct 11 peter 767     /// \deprecated Provided for backward compatibility with the 0.7
2580 05 Oct 11 peter 768     /// API. Use RNG::instance()->seed_from_devurandom() instead.
425 07 Dec 05 jari 769     ///
2578 04 Oct 11 peter 770     unsigned long seed_from_devurandom(void) YAT_DEPRECATE;
425 07 Dec 05 jari 771
425 07 Dec 05 jari 772     ///
424 07 Dec 05 jari 773     /// @return A random number
371 05 Aug 05 peter 774     ///
2901 13 Dec 12 peter 775     virtual result_type operator()(void) const = 0;
371 05 Aug 05 peter 776
424 07 Dec 05 jari 777   protected:
648 14 Sep 06 peter 778     /// pointer to GSL random generator
424 07 Dec 05 jari 779     RNG* rng_;
361 04 Aug 05 jari 780   };
361 04 Aug 05 jari 781
718 26 Dec 06 jari 782   // ContinuousUniform is declared before ContinuousGeneral to avoid
718 26 Dec 06 jari 783   // forward declaration
424 07 Dec 05 jari 784   ///
424 07 Dec 05 jari 785   /// @brief Uniform distribution
424 07 Dec 05 jari 786   ///
424 07 Dec 05 jari 787   /// Class for generating a random number from a uniform distribution
424 07 Dec 05 jari 788   /// in the range [0,1), i.e. zero is included but not 1.
424 07 Dec 05 jari 789   ///
424 07 Dec 05 jari 790   /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
424 07 Dec 05 jari 791   /// Expectation value: 0.5 \n
424 07 Dec 05 jari 792   /// Variance: \f$ \frac{1}{12} \f$
424 07 Dec 05 jari 793   ///
424 07 Dec 05 jari 794   class ContinuousUniform : public Continuous
424 07 Dec 05 jari 795   {
424 07 Dec 05 jari 796   public:
718 26 Dec 06 jari 797     double operator()(void) const;
424 07 Dec 05 jari 798   };
361 04 Aug 05 jari 799
371 05 Aug 05 peter 800   ///
767 22 Feb 07 peter 801   /// @brief Generates numbers from a histogram in a continuous manner.
371 05 Aug 05 peter 802   ///
371 05 Aug 05 peter 803   class ContinuousGeneral : public Continuous
371 05 Aug 05 peter 804   {
371 05 Aug 05 peter 805   public:
371 05 Aug 05 peter 806     ///
371 05 Aug 05 peter 807     /// @brief Constructor
371 05 Aug 05 peter 808     ///
371 05 Aug 05 peter 809     /// @param hist is a Histogram defining the probability distribution
371 05 Aug 05 peter 810     ///
2089 21 Oct 09 peter 811     explicit ContinuousGeneral(const statistics::Histogram& hist);
424 07 Dec 05 jari 812
371 05 Aug 05 peter 813     ///
371 05 Aug 05 peter 814     /// The number is generated in a two step process. First the bin
371 05 Aug 05 peter 815     /// in the histogram is randomly selected (see
371 05 Aug 05 peter 816     /// DiscreteGeneral). Then a number is generated uniformly from
371 05 Aug 05 peter 817     /// the interval defined by the bin.
371 05 Aug 05 peter 818     ///
371 05 Aug 05 peter 819     /// @return A random number.
371 05 Aug 05 peter 820     ///
718 26 Dec 06 jari 821     double operator()(void) const;
371 05 Aug 05 peter 822
371 05 Aug 05 peter 823   private:
371 05 Aug 05 peter 824     const DiscreteGeneral discrete_;
818 17 Mar 07 peter 825     const statistics::Histogram hist_;
371 05 Aug 05 peter 826     ContinuousUniform u_;
371 05 Aug 05 peter 827   };
371 05 Aug 05 peter 828
1342 18 Jun 08 peter 829   /**
1342 18 Jun 08 peter 830      \brief Generator of random numbers from an exponential
1342 18 Jun 08 peter 831      distribution.
2889 07 Dec 12 peter 832
1342 18 Jun 08 peter 833      The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
1342 18 Jun 08 peter 834      \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
1342 18 Jun 08 peter 835      variance \f$ m^2 \f$
1342 18 Jun 08 peter 836   */
424 07 Dec 05 jari 837   class Exponential : public Continuous
424 07 Dec 05 jari 838   {
424 07 Dec 05 jari 839   public:
424 07 Dec 05 jari 840     ///
424 07 Dec 05 jari 841     /// @brief Constructor
424 07 Dec 05 jari 842     ///
648 14 Sep 06 peter 843     /// @param m is the expectation value of the distribution.
424 07 Dec 05 jari 844     ///
2089 21 Oct 09 peter 845     explicit Exponential(const double m=1);
371 05 Aug 05 peter 846
424 07 Dec 05 jari 847     ///
424 07 Dec 05 jari 848     /// @return A random number from exponential distribution.
424 07 Dec 05 jari 849     ///
718 26 Dec 06 jari 850     double operator()(void) const;
424 07 Dec 05 jari 851
424 07 Dec 05 jari 852     ///
424 07 Dec 05 jari 853     /// @return A random number from exponential distribution, with
424 07 Dec 05 jari 854     /// expectation value \a m
424 07 Dec 05 jari 855     ///
424 07 Dec 05 jari 856     /// @note This operator ignores parameters given in constructor.
424 07 Dec 05 jari 857     ///
718 26 Dec 06 jari 858     double operator()(const double m) const;
424 07 Dec 05 jari 859
424 07 Dec 05 jari 860   private:
424 07 Dec 05 jari 861     double m_;
424 07 Dec 05 jari 862   };
424 07 Dec 05 jari 863
1342 18 Jun 08 peter 864   /**
1342 18 Jun 08 peter 865      @brief Gaussian distribution
2889 07 Dec 12 peter 866
1342 18 Jun 08 peter 867      Class for generating a random number from a Gaussian distribution
1342 18 Jun 08 peter 868      between zero and unity. Utilizes the Box-Muller algorithm, which
1342 18 Jun 08 peter 869      needs two calls to random generator.
2889 07 Dec 12 peter 870
1342 18 Jun 08 peter 871      Distribution function \f$ f(x) =
1342 18 Jun 08 peter 872      \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
1342 18 Jun 08 peter 873      \f$ \n
1342 18 Jun 08 peter 874      Expectation value: \f$ \mu \f$ \n
1342 18 Jun 08 peter 875      Variance: \f$ \sigma^2 \f$
1342 18 Jun 08 peter 876   */
424 07 Dec 05 jari 877   class Gaussian : public Continuous
424 07 Dec 05 jari 878   {
424 07 Dec 05 jari 879   public:
424 07 Dec 05 jari 880     ///
424 07 Dec 05 jari 881     /// @brief Constructor
706 19 Dec 06 jari 882     ///
424 07 Dec 05 jari 883     /// @param s is the standard deviation \f$ \sigma \f$ of distribution
648 14 Sep 06 peter 884     /// @param m is the expectation value \f$ \mu \f$ of the distribution
424 07 Dec 05 jari 885     ///
2089 21 Oct 09 peter 886     explicit Gaussian(const double s=1, const double m=0);
424 07 Dec 05 jari 887
424 07 Dec 05 jari 888     ///
424 07 Dec 05 jari 889     /// @return A random Gaussian number
424 07 Dec 05 jari 890     ///
718 26 Dec 06 jari 891     double operator()(void) const;
424 07 Dec 05 jari 892
424 07 Dec 05 jari 893     ///
424 07 Dec 05 jari 894     /// @return A random Gaussian number with standard deviation \a s
424 07 Dec 05 jari 895     /// and expectation value 0.
706 19 Dec 06 jari 896     ///
424 07 Dec 05 jari 897     /// @note this operator ignores parameters given in Constructor
424 07 Dec 05 jari 898     ///
718 26 Dec 06 jari 899     double operator()(const double s) const;
424 07 Dec 05 jari 900
424 07 Dec 05 jari 901     ///
424 07 Dec 05 jari 902     /// @return A random Gaussian number with standard deviation \a s
424 07 Dec 05 jari 903     /// and expectation value \a m.
706 19 Dec 06 jari 904     ///
424 07 Dec 05 jari 905     /// @note this operator ignores parameters given in Constructor
424 07 Dec 05 jari 906     ///
718 26 Dec 06 jari 907     double operator()(const double s, const double m) const;
424 07 Dec 05 jari 908
424 07 Dec 05 jari 909   private:
424 07 Dec 05 jari 910     double m_;
424 07 Dec 05 jari 911     double s_;
424 07 Dec 05 jari 912   };
424 07 Dec 05 jari 913
4281 29 Jan 23 peter 914
1004 23 Jan 08 peter 915   /**
4281 29 Jan 23 peter 916      \brief Multivariate Gaussian Distribution.
4281 29 Jan 23 peter 917
4281 29 Jan 23 peter 918      The distribution is defined by the mean vector and covariance
4281 29 Jan 23 peter 919      matrix.
4281 29 Jan 23 peter 920
4281 29 Jan 23 peter 921      \see <a
4281 29 Jan 23 peter 922      href="https://www.gnu.org/software/gsl/doc/html/randist.html#the-multivariate-gaussian-distribution">The
4281 29 Jan 23 peter 923      Multivariate Gaussian Distribution</a> in the GSL Manual.
4281 29 Jan 23 peter 924
4281 29 Jan 23 peter 925      \since New in yat 0.21
4281 29 Jan 23 peter 926    */
4281 29 Jan 23 peter 927   class MultivariateGaussian
4281 29 Jan 23 peter 928   {
4281 29 Jan 23 peter 929   public:
4281 29 Jan 23 peter 930     /**
4281 29 Jan 23 peter 931        \param m the mean vector
4281 29 Jan 23 peter 932        \param C the covariance matrix that has been processed by a
4281 29 Jan 23 peter 933        CholeskyDecomposer
4281 29 Jan 23 peter 934      */
4281 29 Jan 23 peter 935     MultivariateGaussian(const utility::VectorBase& m,
4281 29 Jan 23 peter 936                          const utility::CholeskyDecomposer& C);
4281 29 Jan 23 peter 937
4281 29 Jan 23 peter 938     /**
4281 29 Jan 23 peter 939        \param m the mean vector
4281 29 Jan 23 peter 940        \param C the covariance matrix that has been processed by a
4281 29 Jan 23 peter 941        CholeskyDecomposer
4281 29 Jan 23 peter 942      */
4281 29 Jan 23 peter 943     MultivariateGaussian(const utility::VectorBase& m,
4281 29 Jan 23 peter 944                          utility::CholeskyDecomposer&& C);
4281 29 Jan 23 peter 945
4281 29 Jan 23 peter 946     /**
4281 29 Jan 23 peter 947        Generate a random vector given the parameters provided in
4281 29 Jan 23 peter 948        constructor.
4281 29 Jan 23 peter 949      */
4281 29 Jan 23 peter 950     void operator()(utility::VectorMutable& x);
4281 29 Jan 23 peter 951
4281 29 Jan 23 peter 952     /**
4281 29 Jan 23 peter 953        Generate a random vector given the passed parameters \c m and
4281 29 Jan 23 peter 954        \c S, ignoring the parameters provided in constructor.
4281 29 Jan 23 peter 955     */
4281 29 Jan 23 peter 956     void operator()(const utility::VectorBase& m,
4281 29 Jan 23 peter 957                     const utility::CholeskyDecomposer& S,
4281 29 Jan 23 peter 958                     utility::VectorMutable& x);
4281 29 Jan 23 peter 959   private:
4281 29 Jan 23 peter 960     utility::Vector mean_;
4281 29 Jan 23 peter 961     utility::CholeskyDecomposer cholesky_decomposer_;
4281 29 Jan 23 peter 962   };
4281 29 Jan 23 peter 963
4281 29 Jan 23 peter 964
4281 29 Jan 23 peter 965   /**
1004 23 Jan 08 peter 966      \brief Convenience function to shuffle a range with singleton RNG.
1004 23 Jan 08 peter 967
1004 23 Jan 08 peter 968      Wrapper around std::random_shuffle using DiscreteUniform as
1004 23 Jan 08 peter 969      random generator and thereby using the underlying RNG class,
1004 23 Jan 08 peter 970      which is singleton.
2145 15 Jan 10 peter 971
3278 05 Jul 14 peter 972      Type Requirements:
3518 05 Oct 16 peter 973      - RandomAccessIterator is \random_access_iterator
1004 23 Jan 08 peter 974    */
2142 15 Jan 10 peter 975   template<typename RandomAccessIterator>
2142 15 Jan 10 peter 976   void random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
1004 23 Jan 08 peter 977   {
1004 23 Jan 08 peter 978     DiscreteUniform rnd;
4252 18 Nov 22 peter 979     std::shuffle(first, last, std::default_random_engine(rnd()));
1004 23 Jan 08 peter 980   }
1004 23 Jan 08 peter 981
680 11 Oct 06 jari 982 }}} // of namespace random, yat, and theplu
424 07 Dec 05 jari 983
361 04 Aug 05 jari 984 #endif

Generated on Sat Oct 14 02:23:03 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)