# Why is this probability distribution biased towards even numbers?

I am working on the probabilities of a game where

1. they pick 20 unique number(non-repeating) from range 1 to 80

2. sort those 20 numbers in ascending order

3. sum every 4 numbers and take the last digit as result (so end with 5 numbers).

Normally I could find these probabilities with rules based on combination and permutation but for this one because of the sorting I am not sure how to work on it.

I tried to brute force all combinations but with some test it would take my old laptop about few thousand years. So I did some simulations, sampling 100 million samples. Below is the final result.

So there seems to be a small edge towards even numbers, but I don't understand where this edge comes from (even numbers appear about .5095925 of the time).

I went further down about the first 4 numbers out of the sorted 20 numbers and try based on the fact that the 1st sorted number must be in the range of 1 to 61, 2nd sorted must be in the range of 2 to 62, etc. So odd positions will have 31 out of 61 chances to be odd and even positions have 31 out of 61 to be even. And I did some checks below, but the resulting bias towards even numbers is very small (.50000004 vs the stimulation .5095925).

I wonder where and how the game tends to favor even numbers quite a bit. Or is my simulation code just wrong?

• +1 For a quick insight into this, think about the results you would get if you were to select 76 numbers rather than 20. It will help to reason in terms of parity: how many odd numbers do you expect to see in each group of four numbers?
– whuber
May 20 '19 at 15:14
• @whuber you mean because there are 40 odd and 40 even, the chance of getting 10 odd and 10 even for the 20 number is higher, then the chance of getting 2 odds and 2 evens (result in even) is higher in a 4 number group base? May 21 '19 at 9:14
• No, that's not it. You need to look at the chances, within any given group, of observing an even number of odd values, because that characterizes the cases where the total will end in 0, 2, 4, 6, or 8. For instance, any group of four numbers in a row is guaranteed to have two odd and two even values and therefore its sum cannot end in an odd digit.
– whuber
May 21 '19 at 12:00
• Your results look good. In particular, the first-number percentages for digits 0 - 5 are typically correct in the first four or five significant digits.
– whuber
May 23 '19 at 17:16

This is an interesting question because the result is a little surprising; it is subtle; and it has been well researched with a large simulation.

Here is a graphical representation of ten iterations of the game. Each row is an iteration. It shows the 20 numbers selected, $$k_1,k_2, \ldots, k_{20},$$ and uses colors to depict their partition into five groups of four. The game sums the numbers in each group, thereby producing five sums.

Nevertheless, within the large number of simulations reported in the question, there appears to be an imbalance among odd and even sums, although this procedure seemingly would not favor one parity. Indeed, if we were only to take a subset of four numbers from $$\{1,2,3,\ldots, 80\}$$ and sum them modulo ten, each digit would have exactly $$1/10$$ chance of appearing.

To appreciate what's happening, consider extreme versions of this procedure, where instead of partitioning a sample of $$20$$ into five groups of four, let's take a sample of $$76$$ and partition it into $$19$$ groups of four:

What strikes us first about this image are the prominent gray "bubbles" depicting the numbers that were not sampled. With only $$80-76=4$$ bubbles in each iteration, we're unlikely to find any bubbles in the first group of four. Moreover, when there is a bubble in that group, as in the fourth row from the top, it's just as likely to be an even number as an odd number. That means the first group (in this extreme case) is highly likely to include two odd numbers and two even numbers, making its sum even.

This consideration of bubbles suggests analyzing the game in terms of gaps between the selected numbers. The first gap is $$\Delta_1 = k_1-0,$$ the second gap is $$\Delta_2 = k_2-k_1,$$ and so on. All gaps are $$1$$ or larger. Evidently the last number chosen is $$k_{20} = \Delta_1+\Delta_2 + \cdots + \Delta_{20},$$ which therefore must be $$80$$ or smaller.

Let's begin with the first gap, $$\Delta_1=k_1.$$ What is the chance this gap is only $$1,$$ the smallest possible value? That's easy to compute, because

1. There are $$\binom{80}{20}$$ distinct, equally-probable subsamples of size $$20$$ that can be obtained from all $$80$$ numbers.

2. Of these, $$\binom{79}{20}$$ are from the set $$\{2,3,\ldots, 80\}.$$

Subtracting, we find the number of subsamples with $$k_1=1$$ and divide that by the total number of subsamples to obtain the probability:

$$\Pr(\Delta_1 = 1) = \frac{\binom{80}{20} - \binom{79}{20}}{\binom{80}{20}} = \frac{\binom{79}{19}}{\binom{80}{20}} = \frac{1}{4}.$$

(Although only one of the ten iterations in the first figure has $$\Delta_1=1$$--the fifth row from the top--a longer set of simulations confirms that $$\Delta_1=1$$ in about a quarter of them.)

Analogous reasoning establishes the remaining chances:

$$\Pr(\Delta_1 = j_1) = \frac{\binom{80+1-j_1}{20} - \binom{80-j_1}{20}}{\binom{80}{20}} = \frac{\binom{80-j_1}{19}}{\binom{80}{20}},\ j_1=1, 2, \ldots, 61.$$

Now suppose we have observed $$\Delta_1.$$ What could $$\Delta_2$$ be? If it is equal to $$j_2,$$ that means the last $$20-1$$ numbers are a subset of the values from $$\Delta_1+j_2$$ through $$80$$ and the remaining $$20-2$$ numbers all lie between $$\Delta_1 + j_2 +1$$ and $$80.$$ In effect, the role of $$80$$ is now played by $$80-j_1,$$ $$20$$ is reduced to $$20-1=1,$$ and $$j_2$$ plays the role of $$j_1.$$ Thus, writing $$k_2=j_1+j_2$$ for the second lowest value selected,

\eqalign{ \Pr(k_2\text{ selected next}\mid k_1\text{ selected first}) &= \Pr(\Delta_2 = j_2 \mid \Delta_1=j_1) \\ &= \frac{\binom{80-k_2}{20-2}}{\binom{80-k_1}{20-1}},\ k_2=j_1+1, 2, \ldots, 80+1-(20-1).}

The pattern is clear. In particular, to study the first group, we can calculate the probability distribution of the sum of just the lowest four numbers in the subsample by summing over all possible values of $$\Delta_1,\Delta_2,\Delta_3,\Delta_4.$$ There are sufficiently few of these that this computation is feasible.

By the laws of conditional probability, the chance of such a combination $$(j_1,j_2,j_3,j_4)$$ is obtained by the products of the component probabilities

\eqalign{ q_{j_4,j_3,j_2,j_1} &= \Pr(\Delta_4=j_4\mid \Delta_3=j_3,\Delta_2=j_2,\Delta_1=j_1) \\ &\times \Pr(\Delta_3=j_3\mid \Delta_2=j_2,\Delta_1=j_1) \\ &\times \Pr(\Delta_2=j_2\mid \Delta_1=j_1) \\ &\times \Pr(\Delta_1=j_1)}

as

$$\Pr(j_1,j_2,j_3,j_4) = \sum_{j_1=1}^{81-20}\ \sum_{j_2=1}^{81-j_1-19}\ \sum_{j_3=1}^{81-j_1-j_2-18}\ \sum_{j_4=1}^{81-j_1-j_2-j_3-17} q_{j_4,j_3,j_2,j_1} .$$

For the game described in the question, there are "only" 635,376 possible combinations in this quadruple sum: the calculations can be carried out in seconds. Here is the interesting part of the probability distribution of the group sum $$k_1+k_2+k_3+k_4 = 4j_1+3j_2+2j_3+j_1:$$

(The full distribution extends to the largest possible sum that could appear in this first group of four numbers, equal to $$61+62+63+64=250.$$)

The solid red dots are probabilities a little larger than suggested by their neighbors: the alternation between even and odd is clear.

When we further reduce this distribution by taking just the last digits of the sums, the alternation between even and odd persists:

This plot shows the relative differences between the probabilities and $$1/10.$$ The actual values (to eight digits) are

Digit Probability
0  0.10235369
1  0.09779473
2  0.10159497
3  0.09859690
4  0.10239713
5  0.09838777
6  0.10201990
7  0.09821967
8  0.10121774
9  0.09741751


Many of them agree closely (to four or five digits) with the simulation results reported in the question.

Notice (from the equations or the simulations) that the first group of four smallest numbers tends to lie to the left of $$20$$ or so, leaving about $$60$$ values from which to select the next four groups. Thus, the situation with the second group is typically just like the situation with five groups, but with $$80$$ reduced to around $$60$$ and only $$20-4=16$$ numbers to select. Similar reasoning suggests the situation with the third and fourth groups also is similar. The simulations in the question bear this out: they still favor even sums over odd.

Analogous discrepancies in probabilities occur when representing the sums in other (non-decimal) bases. Different patterns of discrepancies occur when taking groups of other sizes, especially odd sizes: groups of size $$4$$ are special because they have a pronounced tendency to include either zero, two, or four odd numbers, which guarantees their sum is even.

Here is the R used to do the calculations. It is coded so you can experiment with versions of this game simply by changing the parameters n, m, and g.

library(data.table)
n <- 80 # Maximum number available (100 or more may be problematic)
m <- 4  # Group sizes >= 2 (8 or more is problematic)
g <- 5  # Number of groups >= 1 (more is easier!)
#
# Generate all possible vectors of gaps.
#
delta <- 1:(n - m*(g-1) + 1)
X <- do.call("expand.grid", lapply(1:m, function(i) delta))
names(X) <- paste0("k", 1:m)
X <- as.data.table(X)
i <- rowSums(X) < max(delta)
X <- X[i, ]
#
# Generate the numbers corresponding to each gap vector.
#
Y <- apply(X, 1, cumsum)
rownames(Y) <- paste0("k", 1:m)
Y <- as.data.table(t(Y))
#
# Compute log conditional probabilities for each selection.
#
f <- function(k, n, g) {
m <- length(k)
k0 <- c(0, k[-m])
j <- (m*g-1) : (m*(g-1))
lchoose(n - k, j) - lchoose(n - k0, j+1)
}
Z <- apply(Y, 1, function(k) f(k, n, g))
rownames(Z) <- paste0("p", 1:m)
Z <- as.data.table(t(Z))
#
# Compute the probability of each vector.
#
Z[, pi := exp(rowSums(Z))]
Y[, Sum := rowSums(Y)]
X <- cbind(X, Y, Z)
1-sum(X\$pi) # Should be near zero, within floating point roundoff error
#
# Summarize by group sum.
#
H <- X[, .(Probability = sum(pi)), keyby=Sum]
H.mod <- rbindlist(lapply(2:40, function(b)
H[, .(Probability = sum(Probability), Base=b),
keyby=(Sum - choose(m+1,2)) %% b]))

with(H[Probability > 0.00025],  {
above <- sapply(1:length(Sum), function(i) {
g <- splinefun(Sum[-i], Probability[-i])
Probability[i] > g(Sum[i])
})
colors <- ifelse(above, "Red", "Gray")
plot(Sum, Probability, type="l")
abline(h=0)
abline(v = seq(0, max(Sum), by=5), col="Gray", lty=3)
abline(v = seq(0, max(Sum), by=10), col="Gray")
points(Sum, Probability, pch=21, bg=colors)
})
#
# Explore the last digit in various bases.
#
base <- 10
with(H.mod[Base==base], {
plot(Sum, base*Probability - 1, type="b", ylim=base*range(Probability)-1,
pch=21, bg=ifelse(base*Probability-1 > 0, "Red", "Gray"),
xlab=paste("Sum modulo", base, "offset by", choose(m+1,2)))
abline(h=0, col="Gray")
})

• I am moved by the details and visualization that you offer. There is still much for me to dig into the step by step forming of the formula and the code but first and formal let me thank you for this amazing answer. May 24 '19 at 21:23