Box Muller Transform在实现正态分布PRNG时的问题



使用线性一致生成器我能够创建两个均匀分布的独立伪随机数序列。我正在尝试更改我的程序,以允许其使用这些序列并执行盒装变换以将其更改为正态分布的集合。

我遇到的问题是,无论两个均匀序列的输入种子值如何,我的新"正态分布随机数"(z)总是等于零。

任何技巧都将不胜感激。非常感谢

#define _USE_MATH_DEFINES
#include <iostream>
#include <cmath>
#include <math.h>
using namespace std;
#define M 4294967295
unsigned long get_rand(unsigned long x)                                                                 //establishing function to generate random numbers
{
    unsigned long a = 2269477;
    unsigned long b = 1;                                                                                //Values taken from wikipedia for Linear Congruence Method
    unsigned long m = M;
    unsigned long y;
    y = (a * x + b) % m;
    return y;
}
unsigned long get_normal(unsigned long x1, unsigned long x2)
{
    unsigned long R;
    unsigned long phi;
    unsigned long u;
    R = sqrt(-2 * log(x1));                                                                             //Box-Muller Transform
    phi = (2 * M_PI*x2);
    u = R*cos(phi);
    return u;
}
double u1, u2, Z;
double bin0 = 0;
double bin1 = 0;
double bin2 = 0;                                                                                        //Variables used to store frequency of each number range
double bin3 = 0;
double bin4 = 0;
double bin5 = 0;
double bin6 = 0;
double bin7 = 0;
double bin8 = 0;
double bin9 = 0;
int  main() {
    double seed1,seed2;
    cout << "Please enter seed values " << endl;        
    cin >> seed1;
    cout << "n";
    cin >> seed2;
    double x;
    cout << "Please enter how many random numbers you want " << endl;       
    cin >> x;
    cout << endl;
    cout << "Random Numbers generated shown below: " << endl;

    for (int i = 0; i < x; i++)                                                                         //generate as many random numbers as the user has asked for 
    {
        seed1 = get_rand(seed1);
        seed2 = get_rand(seed2);
        u1 = (double(seed1) / M);                                                                       //changing to double and dividing by 'M' gets all values between 0 and 1
        cout <<"U1 = " << u1 << endl;                                                                   //type conversion to prevent integer rounding in division
        u2 = (double(seed2) / M);
        cout << "U2 = " << u2 << endl;
        Z = get_normal(u1, u2);
        cout <<"Z = " << Z << endl;
        if (Z >= 0.0 && Z <= 0.1)
        {                                                                                               //checking for which intervals each random number falls into
            bin0++;                                                                                     //if a number falls into this interval, increment the counter by 1 each time
        }
        else if (Z > 0.1 && Z <= 0.2)                                                                   //if it doesnt fall into first interval, it will check the next interval, and so on...
        {
            bin1++;
        }
        else if (Z > 0.2 && Z <= 0.3)
        {
            bin2++;
        }
        else if (Z > 0.3 && Z <= 0.4)
        {
            bin3++;
        }
        else if (Z > 0.4 && Z <= 0.5)
        {
            bin4++;
        }
        else if (Z > 0.5 && Z <= 0.6)
        {
            bin5++;
        }
        else if (Z > 0.6 && Z <= 0.7)
        {
            bin6++;
        }
        else if (Z > 0.7 && Z <= 0.8)
        {
            bin7++;
        }
        else if (Z > 0.8 && Z <= 0.9)
        {
            bin8++;
        }
        else if (Z > 0.9 && Z <= 1.0)
        {
            bin9++;
        }
    }
    double binTotal = bin0 + bin1 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 + bin8 + bin9;
    cout << endl;
    int bin0Percent = (bin0 / binTotal) * 100;                                                          //working out a percentage 
    cout << " Number of values in range 0.0-0.1:      " << bin0 << endl;                                //output screen for each interval
    cout << " Percentage of values in this interval:   " << bin0Percent << "%" << endl;
    cout << endl;
    int bin1Percent = (bin1 / binTotal) * 100;
    cout << " Number of values in range 0.1-0.2:      " << bin1 << endl;
    cout << " Percentage of values in this interval:   " << bin1Percent << "%" << endl;
    cout << endl;
    int bin2Percent = (bin2 / binTotal) * 100;
    cout << " Number of values in range 0.2-0.3:      " << bin2 << endl;
    cout << " Percentage of values in this interval:   " << bin2Percent << "%" << endl;
    cout << endl;
    int bin3Percent = (bin3 / binTotal) * 100;
    cout << " Number of values in range 0.3-0.4:      " << bin3 << endl;
    cout << " Percentage of values in this interval:   " << bin3Percent << "%" << endl;
    cout << endl;
    int bin4Percent = (bin4 / binTotal) * 100;
    cout << " Number of values in range 0.4-0.5:      " << bin4 << endl;
    cout << " Percentage of values in this interval:   " << bin4Percent << "%" << endl;
    cout << endl;
    int bin5Percent = (bin5 / binTotal) * 100;
    cout << " Number of values in range 0.5-0.6:      " << bin5 << endl;
    cout << " Percentage of values in this interval:   " << bin5Percent << "%" << endl;
    cout << endl;
    int bin6Percent = (bin6 / binTotal) * 100;
    cout << " Number of values in range 0.6-0.7:      " << bin6 << endl;
    cout << " Percentage of values in this interval:   " << bin6Percent << "%" << endl;
    cout << endl;
    int bin7Percent = (bin7 / binTotal) * 100;
    cout << " Number of values in range 0.7-0.8:      " << bin7 << endl;
    cout << " Percentage of values in this interval:   " << bin7Percent << "%" << endl;
    cout << endl;
    int bin8Percent = (bin8 / binTotal) * 100;
    cout << " Number of values in range 0.8-0.9:      " << bin8 << endl;
    cout << " Percentage of values in this interval:   " << bin8Percent << "%" << endl;
    cout << endl;
    int bin9Percent = (bin9 / binTotal) * 100;
    cout << " Number of values in range 0.9-1.0:      " << bin9 << endl;
    cout << " Percentage of values in this interval:   " << bin9Percent << "%" << endl;
    cout << endl;
}

get_normal返回一个不能在0到1之间的 long,因为它是一个整数。将函数返回的整数存储为doubleZ)不会神奇地恢复废弃的分数部分。

我认为您应该在get_normal中使用浮点算术(即double s),还可以更改返回类型。

顺便说一句,C 标准库具有许多随机数分布。您可能想使用它而不是尝试编写自己的。

M太大,在 long的极限中。因此,任何long除以该M的划分(或模量)将导致0。perphaps您应该使用unsigned long long

也:

而不是R = sqrt(-2 * log(x1));尝试R= sqrt(fabs(2 * log(x1)));

phi = (2 * M_PI*x2);
u = R*cos(phi);

phi始终是2*pi的倍数,因此cos(phi)=1。

最新更新