[C++]Random库--正态分布

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