jeudi 15 décembre 2016

Issue with Box Muller Transform when implementing a Normal Distribution PRNG

Using a Linear Congruent Generator I am able to create two independent Pseudo-Random number sequences which are uniformly distributed. I am trying to alter my program to allow it to use these sequences and perform a Box-Muller transform to change them into a normally distributed set.

The issue I am having however is that my new "normally distributed random number" (Z) is always equal to zero regardless of the input seed values for the two uniform sequences.

Any tips would be gratefully appreciated. Many Thanks

#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;
}




Aucun commentaire:

Enregistrer un commentaire