Pseudorandom numbers in Gaussian distribution

最新推荐文章于 2022-04-22 15:34:56 发布

Memories off

最新推荐文章于 2022-04-22 15:34:56 发布

阅读量261

点赞数

分类专栏： C++数据结构

本文链接：https://blog.csdn.net/happy_single/article/details/89175664

版权

C++数据结构专栏收录该内容

30 篇文章 0 订阅

订阅专栏

Question:

This project is a continuation of the previous project. Many applications require pseudorandom number sequences that are not uniformly distributed. For example, a program that simulates the birth of babies can use random numbers for the birth weights of the newborns. But these birth weights should have a Gaussian distribution. In a Gaussian distribution, numbers are more likely to fall in intervals near the center of the overall distribution. The exact probabilities of falling in a particular interval can be computed from knowing two numbers: (1) the center of the overall distribution (called the median), and (2) a number called the $standard\enspace{deviation}$ , which indicates how widely spread the distribution appears.
Generating a pseudorandom number sequence with an exact Gaussian distribution can be difficult, but there is a good way to approximate a Gaussian distribution using uniformly distributed random numbers in the range [0, 1). The approach is to generate 12 uniformly distributed pseudorandom numbers, each in the range [0, 1). These numbers are then combined to produce the next number in the Gaussian sequence. The formula to combine the numbers is given here, where sum is the sum of 12 numbers and sd is the desired standard deviation:
$Next\,number\,in\,the\,Gaussian\,sequence$ $=median+(sum-6)\times{sd}$ Add a member function to the random number class, which produces a sequence of pseudorandom numbers with approximate Gaussian distribution.

My answer:

random.h

#ifndef RANDOM_H
#define RANDOM_H

class pseudorandom_generator{
public:
    pseudorandom_generator(int a, int b, int c, int d){
        multiplier = a;
        increment = b;
        modulus = c;
        seed = d;
    };
    void set_parameters(int a, int b, int c, int d);
    int generate_number();
    double generate_double_number();
    double generate_gaussian_sequence();
private:
    int multiplier, seed, increment, modulus;
};

#endif

random.cpp

#include "random.h"
#include <iostream>
#include <cmath>

using namespace std;
int ne = 0;
void pseudorandom_generator::set_parameters(int a, int b, int c, int d){
    multiplier = a;
    increment = b;
    modulus = c;
    seed = d;
}

int pseudorandom_generator::generate_number(){
    seed = (multiplier*seed + increment)%modulus;
    // cout << seed << " ";
    return seed;
}

double pseudorandom_generator::generate_double_number(){
    seed = (multiplier*seed + increment)%modulus;
    double t1 = (double)seed;
    double t2 = (double)modulus;
    // cout << t1/t2 << " ";
    return t1/t2;
}

double pseudorandom_generator::generate_gaussian_sequence(){
    // I think that median is 0.5
    double a[10] = {0};
    int i = 0;
    for(; i < 10; i++){
        a[i] = generate_double_number();
    }
    double avg = 0, sd = 0, s = 0;
    for(i = 0; i < 10; i++){
        avg += a[i];
    }
    avg /= 10;
    for(i = 0; i < 10; i++){
        s += pow(a[i]-avg, 2);
    }
    s /= 10;
    sd = sqrt(s);
    if((0.5 + (10 * avg - 6) * sd) < 0)
        ne++;
    return (0.5 + (10 * avg - 6) * sd);
    // return (avg + (10 * avg - 6) * sd);
}

main.cpp

#include "random.h"
#include <iostream>
#include <iomanip>
#include <cassert>

using namespace std;

const int NUM = 10;

enum Range{
    R1, R2, R3, R4, R5, R6, R7, R8, R9, R10
};

unsigned int range[NUM] = {0};

void update_list(double r);
void print_list();

extern int ne;

int main(){
    pseudorandom_generator a(40, 725, 729, 1);
    cout << "Test 1: multiplier = 40, increment = 725, modulus = 729, seed = 1" << endl;
    int i = 0;
    while(i < 1000000){
        // a.generate_number();
        // double r = a.generate_double_number();
        double r = a.generate_gaussian_sequence();
        update_list(r);
        i++;
    }
    print_list();
    cout << ne << endl;
    return 0;
}

void update_list(double r){
    if(r >= 0 && r < 0.1){
        range[R1] += 1;
    } else if(r >= 0.1 && r < 0.2){
        range[R2] += 1;
    } else if(r >= 0.2 && r < 0.3){
        range[R3] += 1;
    } else if(r >= 0.3 && r < 0.4){
        range[R4] += 1;
    } else if(r >= 0.4 && r < 0.5){
        range[R5] += 1;
    } else if(r >= 0.5 && r < 0.6){
        range[R6] += 1;
    } else if(r >= 0.6 && r < 0.7){
        range[R7] += 1;
    } else if(r >= 0.7 && r < 0.8){
        range[R8] += 1;
    } else if(r >= 0.8 && r < 0.9){
        range[R9] += 1;
    } else if(r >= 0.9 && r < 1.0){
        range[R10] += 1;
    }
}

void print_list(){
    cout << "------------------------------------" << endl;
    cout << setiosflags(ios::left) << setw(14) << "Range" << resetiosflags(ios::left)
         << setiosflags(ios::right) << setw(9) << "Number of Occurrences" << resetiosflags(ios::right) << endl;
    cout << "------------------------------------" << endl;
    for(int i = R1; i <= R10; i++){
        switch(i){
            case R1:
                cout << setiosflags(ios::left) << setw(14) << "[0.0, 0.1)" << resetiosflags(ios::left)
                     << setiosflags(ios::right) << setw(9) << range[R1] << resetiosflags(ios::right) << endl;
                break;
            case R2:
                cout << setiosflags(ios::left) << setw(14) << "[0.1, 0.2)" << resetiosflags(ios::left)
                     << setiosflags(ios::right) << setw(9) << range[R2] << resetiosflags(ios::right) << endl;
                break;
            case R3:
                cout << setiosflags(ios::left) << setw(14) << "[0.2, 0.3)" << resetiosflags(ios::left)
                     << setiosflags(ios::right) << setw(9) << range[R3] << resetiosflags(ios::right) << endl;
                break;
            case R4:
                cout << setiosflags(ios::left) << setw(14) << "[0.3, 0.4)" << resetiosflags(ios::left)
                     << setiosflags(ios::right) << setw(9) << range[R4] << resetiosflags(ios::right) << endl;
                break;
            case R5:
                cout << setiosflags(ios::left) << setw(14) << "[0.4, 0.5)" << resetiosflags(ios::left)
                     << setiosflags(ios::right) << setw(9) << range[R5] << resetiosflags(ios::right) << endl;
                break;
            case R6:
                cout << setiosflags(ios::left) << setw(14) << "[0.5, 0.6)" << resetiosflags(ios::left)
                     << setiosflags(ios::right) << setw(9) << range[R6] << resetiosflags(ios::right) << endl;
                break;
            case R7:
                cout << setiosflags(ios::left) << setw(14) << "[0.6, 0.7)" << resetiosflags(ios::left)
                     << setiosflags(ios::right) << setw(9) << range[R7] << resetiosflags(ios::right) << endl;
                break;
            case R8:
                cout << setiosflags(ios::left) << setw(14) << "[0.7, 0.8)" << resetiosflags(ios::left)
                     << setiosflags(ios::right) << setw(9) << range[R8] << resetiosflags(ios::right) << endl;
                break;
            case R9:
                cout << setiosflags(ios::left) << setw(14) << "[0.8, 0.9)" << resetiosflags(ios::left)
                     << setiosflags(ios::right) << setw(9) << range[R9] << resetiosflags(ios::right) << endl;
                break;
            case R10:
                cout << setiosflags(ios::left) << setw(14) << "[0.9, 1.0)" << resetiosflags(ios::left)
                     << setiosflags(ios::right) << setw(9) << range[R10] << resetiosflags(ios::right) << endl;
                break;
        }
    }
}

结果：

Test 1: multiplier = 40, increment = 725, modulus = 729, seed = 1
------------------------------------
Range         Number of Occurrences
------------------------------------
[0.0, 0.1)       127571
[0.1, 0.2)       156383
[0.2, 0.3)       148148
[0.3, 0.4)       131686
[0.4, 0.5)       102883
[0.5, 0.6)        74072
[0.6, 0.7)        50754
[0.7, 0.8)        17831
[0.8, 0.9)         2744
[0.9, 1.0)         1372
186556

可见经过题中给的公式计算，会得到负的伪随机数，此例中有186556个伪随机数。
将公式中的0.5改为0.8，再次执行得到：

Test 1: multiplier = 40, increment = 725, modulus = 729, seed = 1
------------------------------------
Range         Number of Occurrences
------------------------------------
[0.0, 0.1)        26063
[0.1, 0.2)        50755
[0.2, 0.3)       101508
[0.3, 0.4)       127571
[0.4, 0.5)       156383
[0.5, 0.6)       148148
[0.6, 0.7)       131686
[0.7, 0.8)       102883
[0.8, 0.9)        74072
[0.9, 1.0)        50754
8230

看起来大致是服从正态分布的，仅有8230个负数。