21
$\begingroup$

I need to generate 100 random integers in R, where each integer is between 1 and 4 (hence 1,2,3,4) and the mean is equal to a specific value.

If I draw random uniform numbers between 1 and 5 and get floor, I have a mean of 2.5.

x = floor(runif(100,min=1, max=5)) 

I need to fix the mean to 1.9 or 2.93 for example.

I guess I can generate random integers that add to 100 * mean but I don't know how to restrict to random integers between 1 and 4.

$\endgroup$
  • 7
    $\begingroup$ I think this is a bit under-determined... One for instance can get a mean of 1.9 with sample(size=n, x= 1:4, prob=c(3.666,1,1,1), replace=TRUE) but also with sample(size=n, x= 1:4, prob=c(3,1,1,0.715), replace=TRUE). $\endgroup$ – usεr11852 Jan 5 at 22:53
  • 2
    $\begingroup$ Are you asking how to constrain the mean of the underlying distribution, or the sample mean? $\endgroup$ – user20160 Jan 5 at 22:53
  • 1
    $\begingroup$ tha sample mean @user20160 $\endgroup$ – Fierce82 Jan 5 at 23:04
  • 2
    $\begingroup$ Integers between 1 and 4 only allows for 2 and 3. You also need to specify the distribution that they are drawn randomly from (or make one up). $\endgroup$ – wolfies Jan 6 at 6:28
  • 5
    $\begingroup$ I voted to leave this open because there's an interesting algorithmic question in here--the R part is incidental; you could just as easily implement this in Python or with a pad and some dice. $\endgroup$ – Matt Krause Jan 6 at 21:14
41
$\begingroup$

I agree with X'ian that the problem is under-specified. However, there is an elegant, scalable, efficient, effective, and versatile solution worth considering.

Because the product of the sample mean and sample size equals the sample sum, the problem concerns generating a random sample of $n$ values in the set $\{1,2,\ldots, k\}$ that sum to $s$ (assuming $n \le s \le kn,$ of course).

To explain the proposed solution and, I hope, justify the claim of elegance, I offer a graphical interpretation of this sampling scheme. Lay out a grid of $k$ rows and $n$ columns. Select every cell in the first row. Randomly (and uniformly) select $s-n$ of the remaining cells in rows $2$ through $k.$ The value of observation $i$ in the sample is the number of cells selected in column $i:$

Figure

This $4\times 100$ grid is represented by black dots at the unselected cells and colored patches at the selected cells. It was generated to produce a mean value of $2,$ so $s=200.$ Thus, $200-100=100$ cells were randomly selected among the top $k-1=3$ rows. The colors represent the numbers of selected cells in each column. There are $28$ ones, $47$ twos, $22$ threes, and $3$ fours. The ordered sample corresponds to the sequence of colors from column $1$ through column $n=100.$

To demonstrate scalability and efficiency, here is an R command to generate a sample according to this scheme. The question concerns the case $k=4, n=100$ and $s$ is $n$ times the desired average of the sample:

tabulate(sample.int((k-1)*n, s-n) %% n + 1, n) + 1

Because sample.int requires $O(s-n)$ time and $O((k-1)n)$ space, and tabulate requires $O(n)$ time and space, this algorithm requires $O(\max(s-n,n))$ time and $O(kn)$ space: that's scalable. With $k=4$ and $n=100$ my workstation takes only 12 microseconds to perform this calculation: that's efficient.

(Here's a brief explanation of the code. Note that integers $x$ in $\{1,2,\ldots, (k-1)n\}$ can be expressed uniquely as $x = nj + i$ where $j \in \{0,1,\ldots, k-2\}$ and $i\in\{1,2,\ldots, n\}.$ The code takes a sample of such $x,$ converts them to their $(i,j)$ grid coordinates, counts how many times each $i$ appears (which will range from $0$ through $k-1$) and adds $1$ to each count.)

Why can this be considered effective? One reason is that the distributional properties of this sampling scheme are straightforward to work out:

  • It is exchangeable: all permutations of any sample are equally likely.

  • The chance that the value $x \in\{1,2,\ldots, k\}$ appears at position $i,$ which I will write as $\pi_i(x),$ is obtained through a basic hypergeometric counting argument as $$\pi_i(x) = \frac{\binom{k-1}{x-1}\binom{(n-1)(k-1)}{s-n-x+1}}{\binom{n(k-1)}{ s-n}}.$$ For example, with $k=4,$ $n=100,$ and a mean of $2.0$ (so that $s=200$) the chances are $\pi = (0.2948, 0.4467, 0.2222, 0.03630),$ closely agreeing with the frequencies in the foregoing sample. Here are graphs of $\pi_1(1), \pi_1(2), \pi_1(3),$ and $\pi_1(4)$ as a function of the sum:

    Figure 2

  • The chance that the value $x$ appears at position $i$ while the value $y$ appears at position $j$ is similarly found as $$\pi_{ij}(x,y) = \frac{\binom{k-1}{x-1}\binom{k-1}{y-1}\binom{(n-1)(k-1)}{s-n-x-y+2}}{\binom{n(k-1)}{ s-n}}.$$

These probabilities $\pi_i$ and $\pi_{ij}$ enable one to apply the Horvitz-Thompson estimator to this probability sampling design as well as to compute the first two moments of the distributions of various statistics.

Finally, this solution is versatile insofar as it permits simple, readily-analyzable variations to control the sampling distribution. For instance, you could select cells on the grid with specified but unequal probabilities in each row, or with an urn-like model to modify the probabilities as sampling proceeds, thereby controlling the frequencies of the column counts.

| cite | improve this answer | |
$\endgroup$
  • 7
    $\begingroup$ (+1) Ultimate elegance, indeed. $\endgroup$ – Xi'an Jan 6 at 18:56
  • 3
    $\begingroup$ The answer is too difficult for me to follow, appreciate it nonetheless $\endgroup$ – Fierce82 Jan 7 at 11:33
  • 2
    $\begingroup$ What an elegant and beautifully presented answer. If you don't mind my humble suggestion as a reader, you might consider presenting the solution first (the counting patches and the great diagram), and then talking about the implementation and how your argument about how it fits the intuition, and finally why it's efficient. It might make it a bit easier to follow. $\endgroup$ – Neil G Jan 8 at 8:49
  • 3
    $\begingroup$ @Neil Thank you for your suggestion. I think it's a good one and will consider it carefully. $\endgroup$ – whuber Jan 8 at 15:02
  • $\begingroup$ This is a lovely and satisfying answer. I did want to note that the numbers are small enough in this case (100 numbers summing to 190) that we can calculate the uniform distribution of all values that satisfy. I ran some calculations to compare your distribution against this and found that yours is much much more likely (billions in some cases) to select small non-1 values. For example, your model will almost never give distributions with >45 "ones" (~0.002% chance for 46, vanishing for more), but that comprises ~58% of the uniform model values. $\endgroup$ – Cireo Jan 31 at 5:29
15
$\begingroup$

The question is under-specified in that the constraints on the frequencies \begin{align}n_1+2n_2+3n_3+4n_4&=100M\\n_1+n_2+n_3+n_4&=100\end{align} do not determine a distribution: "random" is not associated with a particular distribution, unless the OP means "uniform". For instance, if there exists one solution $(n_1^0,n_2^0,n_3^0,n_4^0)$ to the above system, then the distribution degenerated at this solution is producing a random draw that is always $(n_1^0,n_2^0,n_3^0,n_4^0)$.

In the case the question is about simulating a Uniform distribution over the grid\begin{align}n_1+2n_2+3n_3+4n_4&=100M\\n_1+n_2+n_3+n_4&=100\end{align}one can always use a Metropolis-Hastings algorithm. Starting from $(n_1^0,n_2^0,n_3^0,n_4^0)$, create a Markov chain by proposing symmetric random perturbations of the vector $(n_1^t,n_2^t,n_3^t,n_4^t)$ and accept if the result is within $\{1,2,3,4\}^4$ and satisfies the constraints.

For instance, here is a crude R rendering:

cenM=293
#starting point (n¹,n³,n⁴)
n<-sample(1:100,3,rep=TRUE)
while((sum(n)>100)|(n[2]-n[1]+2*n[3]!=cenM-200))
    n<-sample(1:100,3,rep=TRUE)
#Markov chain
for (t in 1:1e6){
  prop<-n+sample(-10:10,3,rep=TRUE)
  if ((sum(prop)<101)&
      (prop[2]-prop[1]+2*prop[3]==cenM-200)&
      (min(prop)>0)) 
        n=prop}
c(n[1],100-sum(n),n[-1])

with the distribution of $(n_1,n_3,n_4)$ over the 10⁶ iterations:

enter image description here

In case you want draws of the integers themselves,

 sample(c(rep(1,n[1]),rep(2,100-sum(n)),rep(3,n[2]),rep(4,n[3])))

is a quick & dirty way to produce a sample.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ thanks. but I cannot understand how i can utilize this to get the 4 integers (between 1 and 4) $\endgroup$ – Fierce82 Jan 7 at 11:37
  • $\begingroup$ This generates the numbers of 1,2,3,4's $n_1,n_2,n_3,n_4)$ so that there are 100 of them and the sum is cenM. The integer themselves are a random permutation of $n_1$ 1's,..., $n_4$ 4's. $\endgroup$ – Xi'an Jan 8 at 7:35
10
$\begingroup$

I want to ... uh ... "attenuate" @whuber's amazing answer, which @TomZinger says is too difficult to follow. By that I mean I want to re-describe it in terms that I think Tom Zinger will understand, because it's clearly the best answer here. And as Tom gradually uses the method and finds that he needs, say, to know the distribution of the samples rather than just their mean, whuber's answer will be just what he's looking for.

In short: there are no original ideas here, only a simpler explanation.

You'd like to create $n$ integers from $1$ to $4$ with mean $r$. I'm going to suggest computing $n$ integers from $0$ to $3$ with mean $r-1$, and then adding one to each of them. If you can do that latter thing, you can solve the first problem. For instance, if we want 10 integers between $1$ and $4$ with mean $2.6$, we can write down these $10$ integers between $0$ and $3$...

0,3,2,1,3,1,2,1,3,0

whose mean is $1.6$; if we increase each by $1$, we get

1,4,3,2,4,2,3,2,4,1

whose mean is $2.6$. It's that simple.

Now let's think about the numbers $0$ through $3$. I'm going to think of those as "how many items do I have in a 'small' set?" I might have no items, one item, two items, or three items. So the list

0,3,2,1,3,1,2,1,3,0

represents ten different small sets. The first is empty; the second has three items, and so on. The total number of items in all the sets is the sum of the ten numbers, i.e., $16$. And the average number of items in each set is this total, divided by $10$, hence $1.6$.

whuber's idea is this: suppose you make yourself ten small sets, with the total number of items being $10t$ for some number $t$. Then the average size of the sets will be exactly $t$. In the same way, if you make yourself $n$ sets with a total number of items being $nt$, the average number of items in a set will be $t$. You say you're interested in the case $n = 100$.

Let's make this concrete for your example: you want 100 items between 1 and 4 whose average is $1.9$. Using the idea of my first paragraph, I'm going to change this to "make $100$ ints between $0$ and $3$ whose average is $0.9$". When I'm done, I'll add $1$ to each of my ints to get a solution to your problem. So my target average is $t = 0.9$.

I want to make $100$ sets, each with between $0$ and $3$ items in it, with an average set-size of $0.9$.

As I've observed above, this means that there have to be a total of $100 \cdot 0.9 = 90$ items in the sets. From the numbers $1, 2, \ldots, 300$, I'm going to select exactly $90$. I can indicate the selected ones by making a list of 300 dots and Xs:

..X....X...XX...

where the list above indicates that I selected the numbers 3, 9, 13, 14, and then many others that I haven't shown because I got sick of typing. :) I can take this sequence of 300 dots and Xs and break it into three groups of 100 dots each, which I arrange one atop the other, getting something that looks like this:

...X....X..X.....X...
.X...X.....X...X.....
..X...X.X..X......X..

but goes on for a full 100 items in each row. The number of Xs in each row might differ -- there might be 35 in the first row, 24 in the second, and 31 in the third, for instance, and that's OK. [Thanks to whuber for pointing out that I had this wrong in a first draft!]

Now look at each column: each column can be considered as a set, and that set has between 0 and 3 "X"s in it. I can write the tallies below the rows to get something like this:

...X....X..X.....X...
.X...X.....X...X.....
..X...X.X..X......X..
011101102003000101100

That is to say, I've produced 100 numbers, each between 1 and 3. And the sum of those 100 numbers must be the number of Xs, total, in all three rows, which was 90. So the average must be $90/100 = 0.9$, as desired.

So here are the steps to getting 100 integers between 1 and 4 whose average is exactly $s$.

  1. Let $t = s - 1$.
  2. Compute $k = 100 t$; that's how many Xs we'll place in the rows, total.
  3. Make a list of 300 dots-or-Xs, $k$ of which are Xs.
  4. Split this into three rows of 100 dots-or-Xs, each containing about a third of the Xs, more or less.
  5. Arrange these in an array, and compute column sums, getting 100 integers between $0$ and $3$. Their average will be $t$.
  6. Add one to each column sum to get 100 integers between $1$ and $4$ whose average is $s$.

Now the tricky part of this is really in step 4: how do you pick $300$ items, $k$ of which are "X" and the other $300-k$ of which are "."? Well, it turns out that R has a function that does exactly that.

And then whuber tells you how to use it: you write

tabulate(sample.int((k-1)*n, s-n) %% n + 1, n)

For your particular case, $n = 100$, and $s$, the total number of items in all the small sets, is $100r$, and you want numbers between $1$ and $4$, so $k = 4$, so $k-1$ (the largest size for a 'small set') is 3, so this becomes

tabulate(sample.int(3*100, 100r-100) %% 100 + 1, n)

or

tabulate(sample.int(3*100, 100*(r-1)) %% 100 + 1, 100)

or, using my name $t$ for $r - 1$, it becomes

tabulate(sample.int(3*100, 100*t) %% 100 + 1, 100)

The "+1" at the end of his original formula is exactly the step needed to convert from "numbers between $0$ and $3$" to "numbers between $1$ and $4$".

Let's work from the inside out, and let's simplify to $n = 10$ so that I can show sample outputs:

tabulate(sample.int(3*10, 10*t) %% 10 + 1, 10)

And let's aim for $t = 1.9$, so this becomes

tabulate(sample.int(3*10, 10*1.9) %% 10 + 1, 10)

Starting with sample.int(3*10, 10*1.9): this produces a list of $19$ integers between $1$ and $30$. (i.e., it solved the problem of picking $k$ numbers out of your total -- $300$ in your real problem, $30$ in my smaller example).

As you'll recall, we want to produce three rows of ten dots-and-Xs each, something like

 X.X.XX.XX.
 XXXX.XXX..
 XX.X.XXX..

We can read this left-to-right-top-to-bottom (i.e., normal reading order) to produce a list of locations for Xs: the first item's a dot; the second and third are Xs, and so on, so our list of locations starts out $1, 3, 5, 6, \ldots$. When we get to the end of a row, we just keep counting up, so for the picture above, the X-locations would be $1, 3, 5, 6, 8, 9, 11, 12, 13, 14, 16, 17, 18, 21, 22, 24, 26, 27, 28$. Is that clear?

Well, whubers code produces exactly that list of locations with its innermost section.

The next item is %% 10; that takes a number and produces its remainder on division by ten. So our list becomes $1, 3, 5, 6, 8, 9, 1, 2, 3, 4, 6, 7, 8, 1, 2, 4, 6, 7, 8$. If we break that into three groups --- those that came from numbers between $1$ and $10$, those that came from numbers from $11$ to $20$, and those that came from numbers $21$ to $30$, we get $1, 3, 5, 6, 8, 9$, then $1, 2, 3, 4, 6, 7, 8,$, and finally $1, 2, 4, 6, 7, 8$. Those tell you where the Xs in each of the three rows are. There's a subtle problem here: if there had been an X in position 10 in the first row, the first of our three lists would have been $1, 3, 5, 6, 8, 9, 0$, and the tabulate function doesn't like "0". So whuber adds 1 to each item in the list to get $2, 4, 6, 7, 9, 10, 1$. Let's move on to the overall computation:

tabulate(sample.int(3*10, 10*1.9) %% 10 + 1, 10)

This asks "for those $30$ numbers, each indicating whether there's an X in some column, tell me how many times each column (from $1$ to $10$ --- that's what the final "10" tells you) appears, i.e., tell me how many Xs are in each column. The result is 0 3 2 2 2 1 3 2 3 1 which (because of the shift-by-one thing) you have to read as "there are no Xs in the 10th column; there are 3 Xs in the first column; there are 2 Xs in the second column," and so on up to "there is one X in the 9th column".

That gives you ten integers between $0$ and $3$ whose sum is $19$, hence whose average is $1.9$. If you increase each by 1, you get ten integers between $1$ and $4$ whose sum is $29$, hence an average value of $2.9$.

You can generalize to $n = 100$, I hope.

| cite | improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ +1 Welcome to our site, John. I appreciate your efforts to explain and clarify these ideas. At one point your description departs from what the code does: one does not divide the three rows into groups of 30 each. Instead, 90 cells out of the 300 cells in those rows are selected. Usually, each row will have a different number of cells. $\endgroup$ – whuber Jan 7 at 17:15
  • 1
    $\begingroup$ Thanks...I actually worried about that a little bit as I wrote it, but I was in mid-sentence, and by the time I was finished, the thought had flown. I'll edit to try to fix it up. $\endgroup$ – John Jan 7 at 21:25
4
$\begingroup$

You can use sample() and select specific probabilities for each integer. If you sum the product of the probabilities and the integers, you get the expected value of the distribution. So, if you have a mean value in mind, say $k$, you can solve the following equation: $$k = 1\times P(1) + 2\times P(2) + 3\times P(3) + 4\times P(4)$$ You can arbitrarily choose two of the probabilities and solve for the third, which determines the fourth (because $P(1)=1-(P(2)+P(3)+P(4))$ because the probabilities must sum to $1$). For example, let $k=2.3$, $P(4)=.1$, and $P(3)=.2$. Then we have that $$k = 1 \times [1-(P(2)+P(3)+P(4)] + 2\times P(2) + 3\times P(3) + 4\times P(4)$$ $$2.3 = [1 - (P(2)+.1+.2)] + 2*P(2) + 3\times .2 + 4\times .1$$ $$2.3 = .7 + P(2) + .6 + .4$$ $$P(2)=.6$$ $$P(1)=1-(P(2)+P(3)+P(4)=1 - (.6+.1+.2)=.1$$

So you can run x <- sample(c(1, 2, 3, 4), 1e6, replace = TRUE, prob = c(.1, .6, .2, .1)) and mean(x) is approximately $2.3$

| cite | improve this answer | |
$\endgroup$
  • 6
    $\begingroup$ This explains how to constrain the mean of the distribution. But, the OP specified in the comments that they want to constrain the sample mean (which won't match the mean of the distribution, except in expectation). On the other hand, it seems the OP accepted this answer anyway, so perhaps that's not what they wanted after all. $\endgroup$ – user20160 Jan 6 at 0:03
  • $\begingroup$ are you sure? @user20160 why sample mean is not contrainted? it's equal to target $\endgroup$ – Fierce82 Jan 6 at 0:08
  • 3
    $\begingroup$ This answer does not provide a way to make the sample mean equal the target value: most of the time the mean will not equal the target. $\endgroup$ – whuber Jan 6 at 0:10
  • 2
    $\begingroup$ @TomZinger Yes. This answer nicely describes how to sample from a distribution with the given target mean. But, the mean of a sample drawn from a distribution will not generally equal the mean of the distribution. $\endgroup$ – user20160 Jan 6 at 4:59
  • 1
    $\begingroup$ I wrote my answer before I saw that comment, but I figured this would be useful anyway. I imagined it would require an integer programming optimization problem to get a sample mean exactly equal to some value. $\endgroup$ – Noah Jan 6 at 8:20
1
$\begingroup$

Here is a simple algorithm: Create $n-1$ random integers in the range $[1,4]$ and calculate the $n^{th}$ integer for the mean to be equal to the specified value. If that number is smaller than $1$ or larger than $4$, then one by one distribute the surplus/lacking onto other integers, e.g. if the integer is $5$, we have $1$ surplus; and we may add this to the next integer if it's not $4$, else add to the next etc. Then, shuffle the entire array.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ One big problem with this proposal is that it doesn't come along with any indication of what the expected frequencies of the resulting values are. $\endgroup$ – whuber Jan 5 at 23:32
  • 2
    $\begingroup$ Although interesting, I thought the OP only requires an algorithm to generate the desired array of integers in a non-deterministic manner. $\endgroup$ – gunes Jan 5 at 23:37
  • 3
    $\begingroup$ I think that avoids the essence of the question rather than providing a satisfactory answer. A good answer should be able to characterize the distribution it proposes in a meaningful way, such as by giving a formula for the probabilities or at least giving the first couple of moments. $\endgroup$ – whuber Jan 5 at 23:46
  • $\begingroup$ A minor adjustment of the simulated data is likely 'proper', however, looking at the expertimental design in cases where more significant mean deviation is required, depending on the intended purpose, could be, from a hypothesis testing perspective, 'suspect', in my judgement. Either over or under loading a random design to justify or reject possible non-random effects that have been actually observed can be questionable practice. So, any method that makes a very small adjustment to the last of say a 100 observations is probably keeping in good practice, in my opinion. $\endgroup$ – AJKOER Jan 7 at 17:33
1
$\begingroup$

As a supplement to whuber's answer, I've written a script in Python which goes through each step of the sampling scheme. Note that this is meant for illustrative purposes and is not necessarily performant.

Example output:

n=10, s=20, k=4

Starting grid
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
X X X X X X X X X X

Filled in grid
X X . . X . X . . X
. . X X X . . . . .
. . . . X X . . . .
X X X X X X X X X X

Final grid
X X . . X . X . . X
. . X X X . . . . .
. . . . X X . . . .
X X X X X X X X X X
2 2 2 2 4 2 2 1 1 2

The script:

import numpy as np

# Define the starting parameters
integers = [1, 2, 3, 4]
n = 10
s = 20
k = len(integers)


def print_grid(grid, title):
    print(f'\n{title}')
    for row in grid:
        print(' '.join([str(element) for element in row]))


# Create the starting grid
grid = []
for i in range(1, k + 1):
    if i < k:
        grid.append(['.' for j in range(n)])
    else:
        grid.append(['X' for j in range(n)])

# Print the starting grid
print_grid(grid, 'Starting grid')

# Randomly and uniformly fill in the remaining rows
indexes = np.random.choice(range((k - 1) * n), s - n, replace=False)
for i in indexes:
    row = i // n
    col = i % n
    grid[row][col] = 'X'

# Print the filled in grid
print_grid(grid, 'Filled in grid')

# Compute how many cells were selected in each column
column_counts = []
for col in range(n):
    count = sum(1 for i in range(k) if grid[i][col] == 'X')
    column_counts.append(count)
grid.append(column_counts)

# Print the final grid and check that the column counts sum to s
print_grid(grid, 'Final grid')
print()
print(f'Do the column counts sum to {s}? {sum(column_counts) == s}.')
| cite | improve this answer | |
$\endgroup$
1
$\begingroup$

I've turned whuber's answer into an r function. I hope it helps someone.

  • n is how many integers you want;
  • t is the mean you want; and
  • k is the upper limit you want for your returned values
whubernator<-function(n=NULL, t=NULL, kMax=5){
  z = tabulate(sample.int(kMax*(n), (n)*(t),replace =F) %% (n)+1, (n))
  return(z)
}

It seems to work as expected:

> w = whubernator(n=10,t=4.2)
> mean(w)
[1] 4.2
> length(w)
[1] 10
> w
 [1] 3 5 3 5 5 3 4 5 5 4

It can return 0s, which matches my needs.

> whubernator(n=2,t=0.5)
[1] 1 0
| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.