# Generating a list of integers that roughly satisfy a distribution

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).

-
Hello, I don't know why you wanted to delete this, but I got a flag from a user that spent a lot of time writing up an answer and wanted to post it, so I undeleted it. If you deleted because it is unclear or you wanted different specifications, then please feel free to edit your post and clarify it. If you deleted it for other reasons, please let me know. – R. M. Jan 6 '14 at 0:28
I didn't think I had described the problem well, and now, thanks to everyone, I have an answer! Thanks for undeleting it. – jtth Jan 6 '14 at 19:48
@jtth As noted your question is tricky to formulate (!). I would write something like this: "I am looking for a list of $n$ integers ($n$ is given) such that samples drawn from that list appear, as much as possible, to be drawn from a given normal distribution". There is certainly a name for this interesting problem, maybe we should ask on a different forum. If that sounds right to you feel free to edit you question. – A.G. Jan 6 '14 at 21:27
@jtth This is an interesting problem, may I ask where you came upon it and if there are areas where you are aware of its use/relevance? -- Thks. – A.G. Jan 8 '14 at 22:47
@A.G. If I want to generate stimuli that vary in particular ways, I can use this to make lists of parameters that follow the kind of distribution we need. This way I have at least some claim about the parametric space of the stimuli, but not the psychological space (scaling, red vs green as different than red vs blue), which is still something. It also lets me generate sets of stimuli that have peaks and valleys, places where I can ask participants to infer what's left based on the properties of the distribution used to generate the stimuli. – jtth Jan 10 '14 at 14:38

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.


-
I did a terrible job describing what I wanted, but you knew. Thank you. – jtth Jan 6 '14 at 19:44

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]}]


-
These are random, though, not deterministic. There are none greater than 41 and it seems heavier on the left than the right. – jtth Jan 5 '14 at 21:25
You need to describe what you are looking for then, because this is essentially the same as your code, only more compact. The process is deterministic as are almost all "random number" generators. – bill s Jan 5 '14 at 21:36

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:

-

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]}]


-
btw, I didn't read the other answers carefully, so I apologise for any repetition. – Aky Jan 6 '14 at 10:13

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}

-