I have done a last-digit test on a draw of numbers between 1 og 370. Team A drew 64 numbers from the pool, and the last digit of those numbers spread out like this:

Digit: 0 1 2 3 4 5 6 7 8 9

Frequency: 9 (14,06 %) 7 (10,94%) 13 (20,31%) 3 (4,69 %) 5 (7,81%) 4 (6,25%) 4 (6,25 %) 8 (12,50 %) 7 (10,94 %) 4 (6,25 %).

Normally one would expect the last digit's to be uniformly distributed between 0 and 9, but in this case the last-digit is 3 in 20,31% of the time.

My question is: What are the chance of this occurring at random? I apologize if this is a to simple question.

  • $\begingroup$ Before anyone can answer, you have to make your question more precise. What is the chance of what happening? That one number gets 20% or more? That exactly this distribution occurs? That a distribution this far from uniform occurs? That the number "3" occurs most often? That the number 3 occurs 20% of the time? Something else? $\endgroup$
    – Peter Flom
    Mar 24 '13 at 12:07
  • $\begingroup$ Sorry for being imprecise. I am asking how to figure out what the probability for one number to get 20% or more is. $\endgroup$
    – Matsuhui
    Mar 24 '13 at 12:23
  • 1
    $\begingroup$ btw also have a look at the en.wikipedia.org/wiki/Benford%27s_law if you are interested in the occurencies of certain numbers. $\endgroup$
    – Le Max
    Mar 24 '13 at 13:00
  • $\begingroup$ Is this sampling without replacement? $\endgroup$ Mar 24 '13 at 13:20

Restatement of the question

Within the set $\{1,2,\ldots, 370\}$ there are $37$ copies of each terminal digit, information we can represent as a vector $x = (x_0, x_1, \ldots, x_9) = (37, 37, \ldots, 37)$ whose subscripts denote the digits and whose values give their counts. Consider sampling $k = 64$ of those numbers without replacement, as suggested in the question.

Evidently $20$% of the time is $20$% of $64$, or $12.8$, lying between $12$ and $13$. I will take the question to ask for the chance that at least one of the nine digits occurs $13$ or more times in the set of $k$ values. To answer this, it's more convenient to compute the complementary chance that no digit appears more than $m=12$ times.


The zero can appear anywhere from $0$ through $12$ times. The chance that it appears exactly $i$ times is the chance that $i$ of the $x_0=37$ zeros were drawn and, independently, that the remaining $k-i = 64-i$ values drawn were among the $370 - 37 = 333$ non-zeros. Because all $\binom{370}{64}$ draws are equally likely, this chance equals


This must be multiplied by the chance that, when $64-i$ numbers are drawn from the $333$ ones, twos, ..., and nines, none appears more than $12$ times. Summing these chances gives a formula. To write it down, all we need do is give a name to the chances we have been discussing. To this end,

$$p(m, k, x)$$

be the chance that when $k$ numbers are drawn (randomly and without replacement) from a set of values as counted by the entries in a vector $x$, none of those numbers will appear more than $m$ times in the sample. We have seen that

$$p(m, k, x) = \sum_{i=0}^{\min(x_0, m, k)}\frac{\binom{x_0}{i}\binom{n-x_0}{k-i}}{\binom{n}{k}}p(m, k-i, (x_1, x_2, \ldots, x_s))$$

where $x = (x_0, x_1, \ldots, x_s)$ ($s+1 = 10$ is the number of distinct possible digits) and $n = x_0 + x_1 + \cdots + x_s = 370$ is the total number of digits to draw from.


An efficient solution is obtained with a dynamic program. It goes through at most $(s+1)(m+1)k$ steps and so can be very fast even for large problems. (A dynamic program effectively performs a recursive calculation, storing intermediate results along the way to avoid recalculating them whenever they might be needed again.)

Here is a Mathematica implementation:

p[m_, k_, x_] := 
  p[m, k, x] =  Module[{n = Total[x], x0 = First[x], y = Rest[x]},
    Sum[(Binomial[x0, i] Binomial[n - x0, k - i]/Binomial[n, k]) p[m, k - i, y],
     {i, 0, Min[x0, m, k]}]];
p[m_, k_, {}] := p[m, k, {}] = 1;
p[m_, k_, x_] /; Total[x] < k := 0


The answer to the question is obtained by the command

p[12, 64, ConstantArray[37, 10]]

with the resulting output (after less than $0.1$ second)


(I have written it out in full to demonstrate that it's unlikely we will find any simple closed formula for the answer.) This is approximately $95.49000$%, indicating the chance of observing $13$ or more of some digit is approximately $4.51000$%. That answers the question.

In the course of this calculation, we will have obtained values of $p(m,k,x)$ for all smaller values of $k$ and all suffixes of $x$, which are conveniently arranged in a two-way table. I use rows for $x$, running from $(x_9)$ to $(x_8,x_9)$ down to the original $(x_0,x_1,\ldots, x_9)$, and columns for $k$ running from $12$ (the smallest value for which any conceivable calculation is needed) through $64$. Here is a plot of the logarithms of $1-p(m,k,x)$ in such an arrangement, with lighter values corresponding to the largest logs:

Array plot

A clearer picture can be obtained by using a third dimension and restricting the heights to the largest values (corresponding to all probabilities of $0.0001$ or greater):

3D plot

These are base-$10$ logs. The surface has been reflected so that its front boundary corresponds to $s=10$ (all $10$ digits) and the leftmost point on that boundary corresponds to our answer (the common log of $0.0451$ is about $-1.3$). Values along the frontmost row of this surface were computed in terms of values along the next row (where only nine digits have to be considered), which in turn were computed in terms of values along the row behind it (eight digits), and so on until all digits were used up.


A small simulation would help verify whether the preceding analysis and implementation are correct. Let's repeat the experiment a million times and tally the results:

x = Mod[Range[370], 10] (* All 370 last digits *);
sim = Table[(Tally@RandomSample[x, 64])[[All, 2]] // Max, {i, 10^6}]; 

This result--which can be expected to be accurate only to about $3.5$ decimal places (as we will see below)--takes considerably longer than computing the exact result ($6.5$ seconds vs $0.08$ seconds), but gives more information. For instance, it estimates the sampling distribution of the maximum digit count:


For instance, it is now evident that observing $16$ or more of any digit in this experiment is quite a rare event and that generally we should expect to see the maximum count range from $8$ through $13$ or so. (In order to compute this distribution exactly, we would have to compute $p(m,k,x)$ separately for $m=0$ through $37 = \max(x)$; there is no overlap among the intermediate results. That's a $5$ second calculation, taking about the same amount of time as this simulation.)

The proportion of trials with counts of $13$ or greater is $95.4785$%. The standard error of this estimate is $\sqrt{0.954785(1-0.954785)/10^6}\approx 0.000208$. The difference between this proportion and the probability computed by $p(m,k,x)$ is only $0.55$ standard errors, indicating that the difference can be attributed to chance: this is evidence that the Mathematica implementation of $p(m,k,x)$ is correct for these arguments.

A simpler version of the problem

The corresponding chances when sampling with replacement are more easily computed. However, as another test of the Mathematica implementation, let's estimate those chances by making the population larger and larger:

Table[p[12, 64, ConstantArray[10^i, 10]] // N, {i, 2, 9}]

Instead of using just $37$ of each final digit, this calculation uses $100$, then $1000$, ..., and finally $10^9$ of each final digit. The output is


Evidently a limiting value of $90.2078$% is being approached. Its complement, $9.792\ldots$%, is the answer we should get in an analysis that assumes sampling with replacement.

Because the answer differs quite a bit from the previous one--it is twice as great--this is a situation where it would be important to be clear about whether the sampling is with or without replacement. Intuitively, $13$ is a sufficiently large fraction of the $37$ copies of each digit available that by the time $12$ copies have been sampled, it is relatively more difficult to get the $13$th into the sample, thereby lowering the probability that $13$ or more copies will be observed (of any digit).


One way to do this would be by simulation. I think this is correct (but it's early in the morning!)

n <- 1000
maxes <- vector("numeric", n)
for (i in 1:n)
x <- round(runif(64, 0, 9))
maxes[i] <- max(table(x)/64)
sum(maxes > 0.20)/n

I got that this happens about 20% of the time, with that seed. You could increase N to get a better estimate.

Please see the comment below.... My code is not quite right.

  • 3
    $\begingroup$ @Peter Flom: I think that there is mistake in Your code, which leads to bad result. After applying round function You don't sample from uniform distribution. You will get digit 0 only when sampled number is lower than 0.5, I have the same concern about 9 - round will return 9 only when number is greater that 8.5. It would be better to use floor(ruinf(64, 0 , 10)) or sample(0:9, 64, replace=TRUE). After the change probability is a bit lower than 0.1. $\endgroup$ Mar 24 '13 at 16:15
  • $\begingroup$ Good catch @TomekTarczynski ! You are right. $\endgroup$
    – Peter Flom
    Mar 24 '13 at 19:44
  • 2
    $\begingroup$ -1. In addition to the unfixed coding error already mentioned, this code also neglects the fact that the you're sampling without replacement from a pre-specified finite collection of items - the numbers 1 up to 370 - which seems clear from the OP's first two sentences. This gives rise to a different distribution from the one you're sampling from - for example, in your code, it is possible to get more than $37$ of a particular digit whereas that is not possible in the experiment described. These things probably explain why your numbers disagree with whuber's (+1) analytic expression. $\endgroup$
    – Macro
    Mar 25 '13 at 12:24

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.