Identical packs of Skittles

Posted on January 9, 2019 by possiblywrong

Introduction

“No two rainbows are the same. Neither are two packs of Skittles. Enjoy an odd mix.” – Skittles label

Analyzing packs of Skittles (or sometimes M&Ms) seems to be a very common exercise in introductory statistics. How many candies are in each pack? How many candies of each individual color? Are the colors uniformly distributed?

The motivation for this post is to ask some questions raised by the claim in the above quote:

How many different possible packs of Skittles are there? Here we consider two packs of Skittles to be distinguishable only by the number of candies of each color.
What is the probability that two randomly purchased packs of Skittles are the same?
What is the expected number of packs that must be purchased until first encountering a duplicate?

Number of possible packs

If there are $n$ candies in a 2.17-ounce pack, and each candy is one of $d=5$ colors, then the number of possible distinguishable packs is

${n+d-1 \choose d-1}$

Interestingly, a Skittles commercial from the 1990s suggests that there are 371,292 different possible packs (about 27 seconds into the video). It’s not clear whether this number is based on any mathematics or just marketing… but it’s actually reasonably close to the value 367,290 that we get with the above formula assuming that there are $n=52$ candies in a pack.

However, the total number of candies varies from pack to pack. Most studies suggest an average of about 60 candies per pack– with 635,376 possible packs of exactly that size– and this study includes a couple of outlier packs with as few as 42 and as many as 85 candies. We can sum the binomial coefficients over this range of possible pack sizes, to get 42,578,514 possible distinguishable packs… but fortunately, we can more conservatively allow for all pack sizes from empty up to, say, 100 candies, and still remain within an order of magnitude, and confidently assert that there are at most 100 million different possible packs of Skittles.

Total number of possible packs of Skittles, from “empty” to given maximum number of candies per pack.

Mars Wrigley Confectionery advertises 200 million individual candies produced every day; if we make the extremely conservative assumption that these are distributed among equal numbers of packs of various sizes (i.e., for very 2.17-ounce pack, there is, for example, a corresponding 3.4-pound party pack), so that the 2.17-ounce packs make up only about 1.6% of total production, this works out to approximately 50,000 packs per day. By the pigeonhole principle, after only about 1800 days of production, there must be two identical packs of Skittles out there somewhere… but identical packs are likely much more common than this worst-case analysis suggests, as we’ll see shortly.

Probability of identical packs

The second question requires a bit more work. Let’s temporarily assume that there are always exactly $n=60$ candies in every pack, so that there are ${64 \choose 4}$ , or 635,376 possible packs. We might naively assume that the probability that two randomly purchased packs are identical is 1/635376. But this is only true if every possible pack is equally likely, e.g., a pack of 60 all-red candies is just as likely as a pack with 12 candies of each color. I can’t envision a method of production and packaging that would yield this uniform distribution.

Instead, let’s assume that each individual candy’s color is independently and identically distributed, with each of the $d=5$ colors equally likely. For example, imagine an equally large number of candies of each color mixed together in one giant urn (I’m a mathematician, so I have to call it an urn), and dispensed roughly 60 at a time into each pack.

This tracks very well with the actual observed distribution of colors of candies within and across packs. Some packs have more yellows, some have more reds, etc., but there is almost never a color missing, and almost never exactly the same number of every color… but with more and more observed packs, the overall distribution of colors is indeed very nearly uniform. In other words, although it can be unintuitive, giving the impression of “designed” non-uniformity of the distribution of colors, we should expect variability using this uniform model, even for a seemingly “large” sample from, say, a party-size bag of over a thousand candies.

So, given two randomly purchased packs each with $n$ candies of $d$ equally likely colors, what is the probability $p(n,d)$ that they are identical? We can compute this probability exactly as the coefficient of an appropriate generating function:

$p(n,d) = \frac{1}{d^{2n}}[\frac{x^{2n}}{(n!)^2}](\sum\limits_{k=0}^n (\frac{x^k}{k!})^2)^d$

For completeness, following is example C++ source code for computing these probabilities using arbitrary-precision rational arithmetic:

This yields, for example, $p(60,5)$ equals approximately 1/10254. In other words, buy two 60-candy packs, and observe whether they are identical. Then buy another pair of 60-candy packs, etc., until you find an identical pair. You should expect to buy about 10,254 pairs of packs on average.

(Aside: This is effectively the same problem as discussed here a couple of years ago.)

However, we must again account for the variability in total number $n$ of candies per pack. That is, $p(n,5)$ is the probability that two packs are identical, conditioned on them both containing exactly $n$ candies. In question (1), we only needed the maximum range of possible values of $n$ , but here we need the actual probability density $f(n)$ , so that we can integrate $f(n)^2 p(n,5)$ over all possible $n$ .

Fortunately, the variance of the approximately normally distributed pack size (we can be reasonably confident in the mean of 60) doesn’t change the answer much: the probability that two randomly purchased packs– of possibly different sizes– will be identical is somewhere between about 1/100,000 and 1/200,000.

Expected number of packs until first duplicate

Finally, question (3) is essentially a messy birthday problem: instead of insisting that a particular pair of Skittles packs be identical, how many packs must we buy on average for some pair to be identical? The numbers from question (2) in the hundred thousands may seem large, but a back-of-the-envelope square-root estimate of “only” about 400-500 packs to find a duplicate was, to me, surprisingly small. This is confirmed by simulation: let’s assume that every pack of Skittles independently contains a (100,0.6)-binomially distributed number of candies– for a mean of 60 candies per pack– with each individual candy’s color independently uniformly distributed. We repeatedly buy packs until we first encounter one that is identical to one we have previously purchased. The figure below shows the distribution of the number of packs we need to buy, based on one million Monte Carlo simulations, where the mean of 524 packs is shown in red.

Distribution of number of packs needed until encountering a first duplicate.

A box of 36 packs of Skittles costs about 21 dollars… so an experiment to search for two identical packs should only cost about $300… on average.

[Edit: See this follow-up post describing the results of this experiment.]

This entry was posted in Uncategorized. Bookmark the permalink.

14 Responses to Identical packs of Skittles

Pingback: Follow-up: I found two identical packs of Skittles, among 468 packs with a total of 27,740 Skittles | Possibly Wrong
nyc says:

April 18, 2019 at 3:24 pm

The number of distinguishable packs is listed as n+d-1 choose d-1. However, the linked entry about multisets of n elements seems to suggest it should be n+d-1 choose d. Typo or subtle difference in the skittles case?

Reply
- possiblywrong says:
  
  April 18, 2019 at 3:43 pm
  
  This is a good question. The variable “n” is playing the same role here and in the Wikipedia page (the total number of Skittles in a pack), and the “d” here corresponds to “x” in the Wikipedia page (the number of possible colors). Translating to the notation here, the Wikipedia formula is C(d+n-1,n). To see that this is equivalent to the formula in this article, note the identity C(h,k)=C(h,h-k).
  
  Reply
wilmington says:

May 6, 2019 at 3:26 pm

Nice article. Perhaps you can tackle how many licks does it take to get to the center of a Tootsie Pop next.

Three questions:
1. In the expression for p(d, n), are you sure that n! belongs in the denominator?
2. What is the function of X**2n in p(d, n)?
3. Can you elaborate on the “back-of-the-envelope square-root estimate”?

Reply
- possiblywrong says:
  
  May 6, 2019 at 4:13 pm
  
  1. Yes, in fact there are two of them, so to speak (the n! is squared), one for each of the two identical n-candy packs we are constructing. Note, however, that in implementation we extract the coefficient of x^(2n), and *multiply* by (n!)^2. (In other words, the coefficient of (x^k/a) is the same as the coefficient of x^k multiplied by a.)
  
  2. Imagine arranging the candies in two equal-size packs in a row, with the first half of the row containing the n candies in the first pack, and the second half containing those in the second pack. In the generating function (i.e., polynomial in formal variable x), each power of x corresponds to a particular *total* number 2n of candies in the row. (And the coefficient indicates the number of possible such arrangements of two n-candy packs with the same color composition.)
  
  3. Suppose, temporarily, that there are exactly M different compositions of Skittles packs, and that they are all equally likely (i.e., the probability that two randomly selected packs are identical is 1/M). Then a very good approximation for the expected number of packs until we first encounter a duplicate is sqrt(pi/2 x M) (see here). In this case, different types of packs are not equally likely, but we can extend this approximation by letting M=1/q, where q is the probability of two identical packs, so that our estimated average is sqrt(pi/(2q)).
  
  Reply
  - wilmington says:
    
    May 7, 2019 at 3:14 pm
    
    Thanks for the reference on sqr(pi/2 x M). Between probability generating functions and that, I’m learning.
    
    But, I still don’t understand how you arrived at p(n,d). I do get that you want to take the coefficient of the X**2n term, but think that the (sum(k=0 to n of X**k / k!)**d would have already given that and don’t see why another factor of x**2n sits outside the sum.
    
    I’m also still not clear on why n! squared is in the denominator rather than numerator. I thought you were essentially treating it as multinomial and would have expected n! squared to be in the numerator (as it ends up in your code).
    
    Suppose we have a cheap knock off of Skittles that has only d=2 colors red and green and you only get n=3 Skittles in a pack. There are 4 possibilities for what you will get in a randomly selected pack. 1/8 of the time, you’ll get 3 red. 1/8 of the time, you’ll get 3 green. 3/8 of the time you’ll get two red and one green and 3/8 of the time you’ll get two green and one red. Square each of these and sum for the probability that two randomly selected bags are the same. 1/64 + 9/64 + 9/64 + 1/64 = 20/64.
    
    Now, let’s use the probability generating function. we have 1/2**(2*3) * (1/3!)**2 * [1 + 2x**2 + x**4 + (8/6)x**6 + (7/12)x**8 + (1/12)x**10 + (1/36)x**12] if I did the math right and ignoring the leading x**2n term that I doubt.
    
    Taking the coefficient of X**6 should give us the probability of randomly selecting two packs with 3 skittles of two possible colors. So that is (1/64) * (1/36) * (8/6) = 1/1728. That’s pretty far off from 20/64.
    
    If I un-ignore the leading x**2n term, and multiply by it, then well be taking the constant term from the power of the sum factor. This will always be 1. In that case, the probability would be (1/64) * (1/36) * 1 = 1/2304. This is also pretty far from 20/64.
    
    So, maybe you could comment on my interpretation and even better, explain how you derived the probability generating function.
possiblywrong says:

May 7, 2019 at 7:25 pm

@wilmington, Ah, I think I see the issue now. There isn’t “another factor of x^(2n),” that’s just notation for coefficient extraction. That is, the expression [x^n]g(x) means, “write out the polynomial g(x), and extract the coefficient of x^n in that polynomial.” For example, [x^2](1+2x+3x^2) equals 3.

See here for a great introduction to generating functions, in particular the definition at the bottom of page 7 for this notation. Note the continuing discussion on the following page 8, which gets to your next question: we could write the “generating function portion” of the formula for p(n,d) in two different ways: (n!)^2[x^(2n)]blah^d, meaning, “write out blah^d, extract the coefficient of x^(2n), and multiply the result by (n!)^2”; or equivalently, [x^(2n)/(n!)^2]blah^d, meaning, “write out blah^d, and extract the coefficient of x^(2n)/(n!)^2.” See the equation (1.2.8) and subsequent example in the link.

Your by-hand calculation of the probability 20/64 for your example with n=3 and d=2 is correct; to use the generating function approach, if we evaluate just the polynomial sum((x^k/k!)^2 for k=0..3)^2, we get (1+x^2+x^4/4+x^6/36)^2, which expanding yields 1+2*x^2+3/2*x^4+5/9*x^6+17/144*x^8+x^10/72+x^12/1296. The coefficient of x^6 is 5/9, times (3!)^2 is 20, divided by 2^6 is 20/64. Alternatively, the coefficient of x^6/((3!)^2) is 20, divided by 2^6 is 20/64.

One final note: this isn’t really a *probability* generating function. There are such things, where *each* coefficient is a probability of the corresponding value of an integer random variable (thus, the coefficients sum to 1). In this case, only the coefficient of x^(2n) (including the additional factors (n!)^2/d^(2n)) is a probability.

Reply
- wilmington says:
  
  May 7, 2019 at 7:37 pm
  
  Ah! Thank you. I should ask for my money back from my college as I didn’t learn that notation as subscripting. But, I was aware that it wasn’t a proper PGF and was only valid at X^(2n). I might have used the term “squirrelly” with colleagues… Thanks for correcting my algebra; I neglected to square the terms of the sum before raising to the d power (plus the misunderstanding about subscripting).
  
  Reply
wilmington says:

May 7, 2019 at 8:10 pm

p.s. Thanks for the reference to generatingfunctionology. Looks pretty. Your use of them seems to be similar to convolution functions.

Reply
Anon Omous says:

December 19, 2019 at 2:19 am

Hi! So, for a class assignment I implemented my own simulation, assuming 60 skittles to a bag. The average number of bags you need to open before finding a duplicate how I have coded it is roughly 127. I then randomized the number of skittles in each bag to between 42 and 85 for each bag “opened” and am now seeing about 830.
I was wondering if there was any way you can point me to the code for the monte carlo simulation where the number of bags to open before a duplicate was calculated at roughly 540 – I am interested in the code. Thank you very much.
snx1

Reply
- possiblywrong says:
  
  December 19, 2019 at 5:30 pm
  
  Note that the approximate average number of packs mentioned in the article is 524, not 540. Beyond that correction, the simulation results you describe are pretty much as expected– as you can see, the probability q of two packs being identical (and thus the number of packs we need to open to find a duplicate pair) depends significantly on whether/how much the number of candies *per pack* might vary. If the density of this number n of candies in a pack is f(n), then the probability q=sum(f(n)^2*p(n,5)), and a very good approximation for the expected number of packs needed to find a duplicate is sqrt(pi/(2q)).
  
  At the end of this initial article, I assumed that f(n) was the PDF of the binomial(100,0.6) distribution; see here for a Monte Carlo simulation in Python making this assumption; the corresponding sqrt(pi/(2q)) estimate is about 527 packs.
  
  In this Python code, you can replace the number of packs np.random.binomial(100,0.6) with the constant 60 to replicate your first experiment; the corresponding sqrt(pi/(2q)) estimate is about 127 packs. Or you can replace with np.random.randint(42,86) to replicate your second experiment; in this case the sqrt(pi/(2q)) estimate is about 836 packs. (Note that this assumes that the number of candies per pack is *uniformly* distributed between 42 and 85, i.e., each pack size is equally likely. See the follow-on article for a look at the actual sample distribution, which is roughly *normally* distributed with a much tighter variance; if we assume that f(n) matches this sample distribution, then the expected number of packs needed is approximately 293.)
  
  Reply
Pingback: NEEEEEEEEEEEERDS! (Actually, Skittles) | Eigenblogger
Jenny Ward says:

April 4, 2020 at 10:15 pm

Hi, degree level maths student here looking at the same problem! Could you explain to me how you knew to use the coefficient of that term of the polynomial found? Just struggling to get my head round it all. Great article.

Reply
- possiblywrong says:
  
  April 5, 2020 at 2:07 pm
  
  @Jenny, see my comment reply here in the follow-up article that describes the generating function in some more detail. Let me know if this doesn’t help answer your question.
  
  Reply