CODEWITHSUNDEEP
c++normal-distribution
user391838
user391838

Reputation: 61

Need a normally-distributed random number generator

I am new to C++ but have a background in mathematics. I am trying to create a random number generator that will spit out a decimal (dollar amount) between 1.00 and 2.00 with an approximately normal distribution (so that the mean is 1.50).

Now, I have never tried anything like this and have not found a similar question that specifically spits out an element from a set of numbers with normally distributed probabilities. In this model, 1.50 would be most likely to come up, 1.00 or 2.00 would have almost no chance of appearing.

It's easy enough to write a p.d.f. for a normal distribution with mean=1.50 and 3*sigma = 0.5 --> sigma = 1/6 (so that nearly all data is between 1.00 and 2.00). However, more than just not knowing how to integrate 101 regions under this curve (which I dont think can be solved analytically) with C++, it doesn't sound efficient to me. I know there is a normal distribution function in C++.

Can somebody please help me write this out, with comments? Thank you

Upvotes: 2

Views: 1409

Answers (3)

einpoklum
einpoklum

Reputation: 131970

The C++ standard library has a normal distribution class - just what you've asked for. Use it - and clip the value so it's between your minimum and maximum:

#include <algorithm> // for std::clamp()
#include <iostream>
#include <random>

int main() {
    std::random_device randomness_device{};
    std::mt19937 pseudorandom_generator{randomness_device()};

    auto mean = 1.5;
    auto std_dev = 0.5;
    auto min_allowed = 1.0;
    auto max_allowed = 2.0;
    std::normal_distribution<> distribution{mean, std_dev};
    auto sample = distribution(pseudorandom_generator);
    auto clamped = 
        // C++17 and later
        std::clamp(sample, min_allowed, max_allowed);
        // C++14 or earlier:
        // std::max(min_allowed,std::min(sample, max_allowed));
        // 

    std::cout 
        << "A value from a normal distribution with mean " << mean
        << " and standard deviation " << std_dev << ": "   << sample
        << "; when clamped to [" << min_allowed << ", " 
        << max_allowed << "], we get: " << clamped << "\n";
}

In terms of the distribution - this changes the measure so that the entire range (-infinity,1) is concentrated at 1, and similarly (2,infinity) is concentrated at 2. As commenters, suggest, there are other ways to interpret your request for an "approximately normal" distribution, such as re-sampling until you hit a value in the required range; or applying a continuous transformation which maps (infinity, infinity) to (1,2), e.g. x -> arctan(x). But you didn't specify what exactly you're after.

Upvotes: 5

Sid S
Sid S

Reputation: 6125

You can use std::normal_distribution<> for this.

It generates samples from the full range of the double data type, but you can throw away samples outside the interval you're interested in and you will still be close to a normal distribution with three sigma.

#include <iostream>
#include <random>

using namespace std;

int main()
{
    auto mean   = 1.5,
         stddev = 1.0 / 6;

    // Create a normal distribution to pull samples from
    // The distribution has mean 1.5 and ~1/6 std dev
    random_device rd;
    mt19937_64 generator(rd());
    normal_distribution<> distribution(mean, stddev);

    cout << "some samples:" << endl;
    for (int i = 0; i < 10; ++i)
    {
        // Generate a sample.  The sample will be in (-infinity, infinity), so 
        // we throw away values that are outside of 3 std devs.
        // The distribution will no longer be normal, but close enough.
        double v;
        do
        {
            v = distribution(generator);
        } while (v < mean - 3 * stddev || v >= mean + 3 * stddev);

        cout << v << endl;
    }
    return 0;
}

some samples:
1.70539
1.49569
1.53731
1.42872
1.34029
1.54886
1.66154
1.54685
1.60833
1.36282

Upvotes: 2

phoxis
phoxis

Reputation: 61940

Adding to the other answers, you can use the (pseudo) unform random number generator to make a normal random number generator by using the Polar method

Sample code:

#include <math.h>
#include <stdlib.h>

double
randn (double mu, double sigma)
{
  double U1, U2, W, mult;
  static double X1, X2;
  static int call = 0;

  if (call == 1)
    {
      call = !call;
      return (mu + sigma * (double) X2);
    }

  do
    {
      U1 = -1 + ((double) rand () / RAND_MAX) * 2;
      U2 = -1 + ((double) rand () / RAND_MAX) * 2;
      W = pow (U1, 2) + pow (U2, 2);
    }
  while (W >= 1 || W == 0);

  mult = sqrt ((-2 * log (W)) / W);
  X1 = U1 * mult;
  X2 = U2 * mult;

  call = !call;

  return (mu + sigma * (double) X1);
}

Next, you can use this to get the normally distributed random numbers for your application.

Upvotes: 3

Related Questions