Discrete distribution not symmetrical

Question

I'm trying to sample a discrete distribution using the std::discrete_distribution function. Here is a mwe:

// discrete_distribution
#include <iostream>
#include <random>

int main()
{
  const int nrolls = 10000; // number of experiments
  const int nstars = 100;   // maximum number of stars to distribute
  std::vector<double> weights;
  weights = {1.28503e-22, 1.67881e-17, 8.99861e-13, 1.70418e-08, 9.27031e-05,
    0.106935, 16.1967, 140.325, 16.1967, 0.106935, 9.27031e-05, 1.70418e-08,
    8.99861e-13, 1.67881e-17, 1.28503e-22};

  std::default_random_engine generator;
  std::discrete_distribution<int> distribution(weights.begin(), weights.end());

  for (double x:distribution.probabilities()) std::cout << x << " ";
  std::cout << std::endl;

  int p[15]={};

  for (int i=0; i<nrolls; ++i) {
    int number = distribution(generator);
    ++p[number];
  }

  std::cout << "a discrete_distribution:" << std::endl;
  for (int i=0; i<15; ++i)
    std::cout << i << ": " << std::string(p[i]*nstars/nrolls,'*') << std::endl;

  return 0;
}

this gives:

7.43082e-25 9.70789e-20 5.20354e-15 9.8546e-11 5.36065e-07 
0.000618363 0.0936591 0.811444 0.0936591 0.000618363 5.36065e-07
9.85459e-11 5.10703e-15 0 0

a discrete_distribution:
0: 
1: 
2: 
3: 
4: 
5: 
6: *********
7: ********************************************************************************
8: *********
9: 
10: 
11: 
12: 
13: 
14:

Note the asymmetry, especially the zeros at the end. I can't see what I've done wrong. Is there something wrong with the code, or is some rounding taking place that I can't see. Thanks.

Answer 1

The problem appears to be floating point math. In particular, it seems likely that the distribution keeps a running total while normalizing its weights, leading it to lose the tiny probabilities at the end. In a simple example, assume that doubles could only store 2 significant digits (the reality is closer to 16) and you weights were 0.001, 1.0, 1.0 and 0.001:

It would sum the weights to 2.002 (which it can only represent as 2.00), then went ahead and normalized the weights. The first one becomes 0.001/2.00 = 0.0005. Then next one is 0.5, running total 0.5005 (i.e. 5.00). The third weight is also 0.5, so the total is 1.00 already. The difference to the allowed sum is 0.00, so it cannot give positive weight to the last event.

I am aware that this is not a perfect example (because the weights still don't fully add up), but I hope you get the point - your standard library implementation and/or your floating point settings are messing up your results here due to cancellation. Not that getting your event of probability 1e-20 to ever happen is within realms of reason, but you are right that it should theoretically keep the symmetry.

For those who say "there was not enough precision to print it": I disagree, because the values should ideally still be symmetric and the first value is not printed as 0, unlike the last. Seeing as only those values smaller than one ULP around 1 are printed as zero, I stand by cancellation as the problem.