Random庫–正態分佈

背景介紹

正態分佈（Normal distribution）又名高斯分佈（Gaussian distribution），是一個在數學、物理及工程等領域都非常重要的概率分佈，在統計學的許多方面有著重大的影響力。若隨機變數X服從一個數學期望為μ、方差為σ^2 的高斯分佈，記為N(μ，σ^2 )。其概率密度函式為正態分佈的期望值μ決定了其位置，其標準差σ決定了分佈的幅度。因其曲線呈鐘形，因此人們又經常稱之為鐘形曲線。我們通常所說的標準正態分佈是μ = 0,σ = 1的正態分佈。

Normal distribution

Random number distribution that produces floating-point values according to a normal distribution, which is described by the following probability density function:

This distribution produces random numbers around the distribution mean (μ) with a specific standard deviation (σ).

The normal distribution is a common distribution used for many kind of processes, since it is the distribution that the aggregation of a large number of independent random variables approximates to, when all follow the same distribution (no matter which distribution).

The distribution parameters, mean (μ) and stddev (σ), are set on construction.

To produce a random value following this distribution, call its member function operator().

實現要求

———- adaptor ———-

about: normalize the random numbers to [0, 1].

adaptor() // default constructor

double min() const // return the minimun that the adaptor
// can express.
double max() const // return the maximun that the adaptor
// can express.
double operator()() // normalize the random numbers to [0, 1].

Tips(pseudo code):
remember, you need to realize adaptor as the following way!

# b = the digits of numeric_limits<double>
# r = generator.max - generator.min + 1
# log2r = log(r) / log(2)
# k = max(1, (b + log2r - 1) / log2r)
# sum = 0 and tmp = 1
# loop k -> 0
# sum += (generator.randomNumber - generator.min) * tmp
# tmp = tmp * r
# end loop
# return sum / tmp

———- normal_distribution_my ———-
about: adapt the random number with adaptor and simulate
the normal distribution.

double mean, stddev; // two parameter: mean and standard deviation.

adaptor aurng; // adaptor to normalize the random number

normal_distribution_my(); // default constructor, mean = 0, stddev = 1.

double min() const; // return the minimun real number it can express.

double max() const; // return the maximun real number it can express.

double operator()(); // use its adaptor to simulate normal distribution.

Tips(pseudo code):
remember, you need to realize normal_distribution as the following way!

# loop until r2 -> (0, 1]
# x = 2 * adaptor - 1
# y = 2 * adaptor - 1
# r2 = x * x + y * y
# end loop
# mult = sqrt(-2 * log(r2) / r2)
# ret = y * mult
# ret = ret * standard deviation + mean
# return ret

程式碼實現

測試檔案

#include <iostream>
#include <cstring>
#include <string>
#include "distribution.h"
using namespace std;
using namespace RAND_GEN;
using namespace DISTRIBUTION;

void test_adaptor() {
    cout << "---------- test_adaptor ----------\n";
    linear_congruential_engine_my lce;
    adaptor adp(lce);
    cout << "min = " << adp.min() << endl;
    cout << "max = " << adp.max() << endl;
    for (int i = 0; i < 19; i++)
        cout << adp() << ", ";
    cout << adp() << endl;
}

void test_normal_distribution() {
    cout << "---------- test_normal_distribution ----------\n";
    const int nrolls = 10000;  // number of experiments
    const int nstars = 100;  // maximum number of stars to distribute
    linear_congruential_engine_my lce;
    normal_distribution_my dist(5.0, 2.0, lce);
    cout << "min = " << dist.min() << endl;
    cout << "max = " << dist.max() << endl;
    int p[10] = { 0 };
    for (int i = 0; i < nrolls; ++i) {
        double number = dist();
        ++p[static_cast<int>(number)];
    }
    cout << "normal_distribution(5.0 , 2.0):\n";
    for (int i = 0; i < 10; ++i) {
        cout << i << "-" << (i + 1) << ": ";
        cout << string(p[i]*nstars/nrolls, '*') << endl;
    }
}

void hard_test() {
    int nrolls, nstars;
    cin >> nrolls >> nstars;
    linear_congruential_engine_my lce;
    double mean, stddev;
    cin >> mean >> stddev;
    normal_distribution_my dist(mean, stddev, lce);
    int p[10] = { 0 };
    for (int i = 0; i < nrolls; ++i) {
        double number = dist();
        ++p[static_cast<int>(number)];
    }
    cout << "normal_distribution(5.0 , 2.0):\n";
    for (int i = 0; i < 10; ++i) {
        cout << i << "-" << (i + 1) << ": ";
        cout << string(p[i]*nstars/nrolls, '*') << endl;
    }
}

int main() {
    int t;
    cin >> t;
    if (t == 0) {
        test_adaptor();
        test_normal_distribution();
    } else {
        hard_test();
    }
    return 0;
}

隨機數生成檔案

#ifndef RANDOM_MY_H
#define RANDOM_MY_H

namespace RAND_GEN {
class mod_my {
  public:
    int m, a, c;
    mod_my(int _m, int _a, int _c) : m(_m), a(_a), c(_c) {}

    //  General case for x = (ax + c) mod m -- use Schrage's algorithm
    //  to avoid integer overflow.
    //  (ax + c) mod m can be rewritten as:
    //    a(x mod q) - r(x / q)  if >= 0
    //    a(x mod q) - r(x / q)  otherwise
    //  where: q = m / a , r = m mod a
    //
    //  Preconditions:  a > 0, m > 0.
    //
    //  Note: only works correctly for m % a < m / a.
    int calc(int x) {
        if (a == 1) {
            x %= m;
        } else {
            int q = m / a;
            int r = m % a;
            int t1 = a * (x % q);
            int t2 = r * (x / q);
            if (t1 >= t2) x = t1 - t2;
            else x = m - t2 + t1;
        }
        if (c != 0) {
            const int d = m - x;
            if (d > c) x += c;
            else x = c - d;
        }
        return x;
    }
};

class linear_congruential_engine_my {
  public:
    int multiplier, increment, modulus;
    unsigned int default_seed_my, seed_my;
    mod_my mod_temp;

    linear_congruential_engine_my()
    : multiplier(16807), increment(1), modulus(2147483647)
    , default_seed_my(1u), mod_temp(modulus, multiplier, increment)
    { seed(default_seed_my); }

    linear_congruential_engine_my(int _m, int _a, int _c, int _s)
    : multiplier(_a), increment(_c), modulus(_m)
    , default_seed_my(_s), mod_temp(modulus, multiplier, increment)
    { seed(default_seed_my); }

    void seed(unsigned int s)
    { seed_my = s; }

    int min()
    { return  increment == 0u ? 1u : 0u; }

    int max()
    { return modulus - 1u; }

    void discard(unsigned long long z)
    { for (; z != 0ULL; --z) (*this)(); }

    int operator()() {
        seed_my = mod_temp.calc(seed_my);
        return seed_my;
    }
};

}  // namespace RAND_GEN

#endif

正態分佈實現：

#ifndef DISTRIBUTION_H
#define DISTRIBUTION_H
#include <limits>
#include <cmath>
#include <algorithm>
#include <cstddef>
#include "random_my.h"
using namespace RAND_GEN;

namespace DISTRIBUTION {

struct adaptor {
  private:
      linear_congruential_engine_my my_gen;
  public:
    adaptor() { }

    explicit adaptor(const linear_congruential_engine_my& gen)
    : my_gen(gen) { }

    double min() const
    { return static_cast<double>(0); }

    double max() const
    { return static_cast<double>(1); }

    double operator()() {
        const size_t b =
        static_cast<size_t>(std::numeric_limits<double>::digits);
        const long double r = static_cast<long double>(my_gen.max())
            - static_cast<long double>(my_gen.min()) + 1.0L;
        const size_t log2r = std::log(r) / std::log(2.0L);
        size_t k = std::max<size_t>(1UL, (b + log2r - 1UL) / log2r);
        double sum = static_cast<double>(0);
        double tmp = static_cast<double>(1);
        for (; k != 0; --k) {
            sum += static_cast<double>(my_gen() - my_gen.min()) * tmp;
            tmp *= r;
        }
        return sum / tmp;
    }
};

class normal_distribution_my {
    double mean, stddev;
    adaptor aurng;
  public:
    normal_distribution_my() : mean(0), stddev(1)
    { }

    normal_distribution_my(double _mean, double _stddev,
      linear_congruential_engine_my &gen)
      : mean(_mean), stddev(_stddev), aurng(gen)
    { }

    double getmean() const
    { return mean; }

    double getstddev() const
    { return stddev; }

    void setparam(const double _mean, const double _stddev) {
        mean = _mean;
        stddev = _stddev;
    }

    double min() const
    { return std::numeric_limits<double>::min(); }

    double max() const
    { return std::numeric_limits<double>::max(); }

    double operator()() {
        double ret, x, y, r2;
        do {  // r2 = (0, 1];
            x = static_cast<double>(2) * aurng() -1.0;
            y = static_cast<double>(2) * aurng() -1.0;
            r2 = x * x + y * y;
        } while (r2 > 1.0 || r2 == 0.0);
        const double mult = std::sqrt(-2 * std::log(r2) / r2);
        ret = y * mult;
        ret = ret * stddev + mean;
        return ret;
    }
};

}  // namespace DISTRIBUTION

#endif