Generating a list of integers that roughly satisfy a distribution

Question

Say I have a normal distribution,

NormalDistribution[25, 7]

I'd like integers from 1 to 50, not sampled at random from that distribution, but sampled in a way that the mean is more common than the bounds, following that distribution.

What is that called, and is there a way of doing it concisely within Mathematica?

I suppose one way to do it would be selecting a number from 1-50 with the probabilities given by

Table[NormalDistribution[25, 7], {x, 1, 50, 1}]
{0.00015967, 0.000257934, 0.000408253, 0.000633121, 0.000962014, \
0.00143223, 0.00208921, 0.00298598, 0.00418147, 0.0057373, 0.007713, \
0.0101596, 0.0131119, 0.0165803, 0.0205426, 0.0249376, 0.0296614, \
0.0345672, 0.0394707, 0.0441593, 0.0484068, 0.051991, 0.0547124, \
0.0564132, 0.0569918, 0.0564132, 0.0547124, 0.051991, 0.0484068, \
0.0441593, 0.0394707, 0.0345672, 0.0296614, 0.0249376, 0.0205426, \
0.0165803, 0.0131119, 0.0101596, 0.007713, 0.0057373, 0.00418147, \
0.00298598, 0.00208921, 0.00143223, 0.000962014, 0.000633121, \
0.000408253, 0.000257934, 0.00015967, 0.0000968449}

Where each position is the probability that that number would be selected: 1 would show up with p = .00015967, etc. Is that right? I'm not really sure how to map and select the integers 1-50 along these probabilities, or even if that's on track. I'm not sure because I don't want it to be left to probability: the function that makes these should be deterministic. I don't think I'm describing this very well.

The idea is to have parameters vary normally around the mean, 25, but are integers and not randomly sampled, as I only need a few (though in this case 50) though that number may change.

I hope this makes some semblance of sense. The end result would be a list

{1, 5, 10, 13, 17, 20, 22, 23, 24, 25, 25, 26, 27, 28, 30, 33, 37, 40, 45, 50}

Or something like that, but with whatever distribution and number of numbers (I'd like to do it for something like MixtureDistribution[{1, 1}, {NormalDistribution[25, 7], NormalDistribution[75, 7]}] as well).

Answer 1

One way to approach your question is to generate a list of integers whose cumulative distribution function (CDF) is as close as possible to that of the normal distribution. For a discrete statistical distribution the CDF is a staircaise function whose steps correspond to the numbers in the list. Ideally the numbers would correspond to where the discrete CDF is just 1/2 under the continuous CDF on the left and 1/2 above on the right. I suggest getting those points and then rounding them to get integers, which should be good enough as long as you have a reasonable number of integers. Here the continuous and discrete CDF's are F and G.

(* your distribution is from 1 to a (a=50), you are looking for n numbers *) 
a = 50; n = 20;
X = CensoredDistribution[{0, a}, NormalDistribution[25, 7]];
F := CDF[X]
cuttingPoints = Range[n]/n - 1/(2 n);
numbers = Round[InverseCDF[X, cuttingPoints]]
Print["Mean = ", N@Mean[numbers], " Standard deviation = ", N@StandardDeviation[numbers]]
(* the discrete CDF *)
G[x_] := Length@Select[numbers, # <= x &]/n
Plot[{F[x], G[x]}, {x, 0, a + 1}, PlotRange -> {0, 1}, PlotPoints -> 200]

Output (note that there are two "25"'s):

{11, 15, 17, 18, 20, 21, 22, 23, 24, 25, 25, 26, 27, 28, 29, 30, 32, 33, 35, 39}
Mean = 25. Standard deviation = 6.98871

Check the results

ShapiroWilkTest[numbers, {"PValue", "TestConclusion"}] // TableForm
QuantilePlot[numbers]
Show[
 Plot[PDF[NormalDistribution[25, 7]][x], {x, 0, a}],
 SmoothHistogram[numbers, PlotStyle -> Red]]

Output:

 1.
 The null hypothesis that the data is distributed according to the 
 NormalDistribution[\[FormalX],\[FormalY]] is not rejected at the 5 percent level
 based on the Shapiro-Wilk test.

Answer 2

I don't really understand what you are looking for, but here is one way to get 50 arbitrary integers with a normal-like distribution centered around 25 with a variance of approximately 7.

list = Round[Abs[RandomVariate[NormalDistribution[25, 7], 50]]]

Here's one way to visualize the result:

Histogram[list, {Range[50]}]

Answer 3

IMHO your question is a bit unclear. I hope the following addresses what you want.

You want a sample from a normal distribution, an unbounded random distribution, and the sample must be in the range of 1..50 but not random, no, it must be "sampled in a way that the mean is more common than the bounds".

Restricting the output from random samplings from a given distribution to a certain range can be done in two ways: using CensoredDistribution or using TruncatedDistribution. The former draws its samples from the specified distribution and any output outside of the specified bounds is set equal to the closest bound. The latter redraws a sample if it lies outside of the bounds and repeats until the it lies within the bounds.

They look like this:

Histogram[#, {1}, "PDF", 
    Epilog -> Line[Table[{x+.5, PDF[NormalDistribution[25, 15], x]}, {x, -20, 70}]],
    PlotRange -> {{-20, 70}, {0, Automatic}}, ImageSize -> 300, 
    Frame -> True, Axes -> None] & /@
  {
   RandomVariate[CensoredDistribution[{1, 50}, NormalDistribution[25, 15]], {100000}],
   RandomVariate[TruncatedDistribution[{1, 50},NormalDistribution[25, 15]], {100000}]
  } // GraphicsRow

The line in the figures depicts the unbounded normal distribution. I used a slightly wider normal distribution than in your question (SD = 15 instead of 7) to better demonstrate the difference. With SD = 7 the probability of a sample to get outside of the bounds is rather low:

Probability[x < 1 || x > 50, x \[Distributed] NormalDistribution[25, 7]] // N

0.0004809031132

So, the unbounded and the bounded distributions are pretty much the same:

The code you displayed in your question doesn't actually generate the output you provided. It should be

Table[PDF[NormalDistribution[25, 7], x], {x, 1, 50, 1}] // N

{0.0001596702666, 0.0002579337301, 0.0004082527082, 0.0006331212017, 0.0009620142107, 0.001432230655, 0.002089206122, 0.002985977025, 0.004181465147, 0.005737297206, 0.007712995216, 0.01015957693, 0.01311188199, 0.01658025814, 0.02054255182, 0.02493758204, 0.02966136545, 0.03456724636, 0.03947074079, 0.0441593444, 0.04840684797, 0.05199096025, 0.05471239428, 0.05641316285, 0.05699175434, 0.05641316285, 0.05471239428, 0.05199096025, 0.04840684797, 0.0441593444, 0.03947074079, 0.03456724636, 0.02966136545, 0.02493758204, 0.02054255182, 0.01658025814, 0.01311188199, 0.01015957693, 0.007712995216, 0.005737297206, 0.004181465147, 0.002985977025, 0.002089206122, 0.001432230655, 0.0009620142107, 0.0006331212017, 0.0004082527082, 0.0002579337301, 0.0001596702666, 0.00009684491216}

Of course, these probabilities don't add to zero because they're not from the distribution's full range, which is $(-\infty,\infty)$. So, you can't have a distributions that produces this. Distributions must sum up to one.

You could use these probabilities to make an EmpiricalDistribution. This distribution can use weights for every value present in the sample space and scales them so they add up to 1. Using the wide normal distribution for demonstration:

Histogram[#, {1}, "PDF", 
   Epilog -> Line[Table[{x+.5, PDF[NormalDistribution[25, 15], x]}, {x, 1, 50}]],
   PlotRange -> {0, Automatic}, ImageSize -> 300] &@
 RandomVariate[
  EmpiricalDistribution[Table[PDF[NormalDistribution[25, 15], x], {x, 1, 50}] -> Range[50]], 
  {100000}
 ]

I suppose this may be the same as the TruncatedDistribution.

The SD=7 normal distribution:

Answer 4

Here's an approach that maps from a continuous probability distribution to a discrete one, using the formula:

$$\int_{i-1/2}^{i+1/2} p(x) dx = P[i]$$

(Here i is from set of all integers. I suspect this has all sorts of problems, rigor-wise, but anyway...)

Let's generate the probabilities of the numbers 1 through 50 using this approach, by doing:

w = Table[NProbability[i - 0.5 <= x <= i + 0.5, 
  Distributed[x, NormalDistribution[25, 7]]], {i, 1, 50}]

(* {0.000161133, 0.000260085, 0.000411338, 0.000637432, 0.000967878, \
0.00143999, 0.00209918, 0.00299841, 0.00419649, 0.00575482, \
0.00773266, 0.0101807, 0.0131335, 0.016601, 0.0205607, 0.0249514, \
0.0296691, 0.0345672, 0.0394618, 0.044141, 0.0483791, 0.0519549, \
0.0546697, 0.0563662, 0.0569433, 0.0563662, 0.0546697, 0.0519549, \
0.0483791, 0.044141, 0.0394618, 0.0345672, 0.0296691, 0.0249514, \
0.0205607, 0.016601, 0.0131335, 0.0101807, 0.00773266, 0.00575482, \
0.00419649, 0.00299841, 0.00209918, 0.00143999, 0.000967878, \
0.000637432, 0.000411338, 0.000260085, 0.000161133, 0.0000978148} *)

By the way, if we do Total[w] we get around .999633, which would mean that the numbers 1 to 50 account for almost 99.97% of the probability mass. (With the remaining integers - a countably infinite set! - accounting for the remaining 0.03%, according to our "model".)

Now generate 10000 (say) random samples from Range[50] using w as a weight:

samples = RandomChoice[w -> Range[50], 10000];

Here's the histogram:

Histogram[samples, {Range[50]}]

Answer 5

Round[Quantile[NormalDistribution[25, 7], Range[20]/21]]

{13, 16, 18, 19, 20, 21, 22, 23, 24, 25, 25, 26, 27, 28, 29, 30, 31, 32, 34, 37}