Yes, this has been discussed before, but I have calculated the probability at 0.513858194042236 and the calculations are on my website: Occupancy Calculator Part 4
(I’m certain the calculations are accurate.)
However, is there an easier way to do this? (I’ve read that it can be done using Stirling Numbers of the second kind but I have not seen any explanation for this.)
Any explanation would be greatly appreciated.
Sure, using a Stirling number of the second kind could be expressed in an “easier” fashion, but a full explanation would require more in-depth knowledge and would ultimately be as mechanically difficult to compute.
The same calculation you want to do could be simplified in this case to:
{13,6}*P(6,6)/(6^13)
where 6^13 is the number of outcomes of 13 rolls of a single 6-sided die and P(6,6) is the number of permutations of 6 things taken 6 at a time (the 6 being the number of faces of the die).
The {13,6} is the Stirling number of the second kind which represents the number of ways to partition the 13 rolls into 6 non-empty subsets. That is, the 6 subsets are the 6 possible outcomes and you need them all nonempty, i.e. one of each face shows up.
How to come up with this number? Well, there’s an explicit formula for Stirling numbers of the second kind, though it’s about as ‘complicated’ as the work you’ve already done. There’s also a recurrence relation, but that can also be mechanically complicated to compute for numbers as large as {13,6}. But the fast way is to look it up. It’s the 84-th entry of A008277, which happens to be 9321312. You can also probably get it via Wolfram alpha or some other computer based calculator.
{13,6} = 9,321,312
P(6,6) = 720
6^13 = 13,060,694,016
9,321,312 *720 / 13,060,694,016 = 0.513858…
Note that the lookup (or computer version) is “easier” in one sense but explaining Stirling numbers may be conceptually more difficult than the calculation you already have.
I would do basically the same calculation but phrase it somewhat differently.
1st approximation would be the one minus the sum of the six probabilities that each number was not rolled: 1 - 6 * (5/6)^13. (this corresponds to line 2 of your step 4)
But that computation double subtracts the results in which two numbers are both not rolled, so the next correction is to add back the probability corresponding to the 15 pairs that were double-subtracted in the first step, or + 15 *(2/3)^13. (this corresponds to line 3 of your step 4)
Next, the results with exactly three unrolled numbers have been completely canceled out by the first two steps, but as they should be subtracted once we’re now overcounting and need to subtract those back out. There are 20 of these, so - 20 * (1/2) ^13.
etc.
Great Antibob - thanks for the reply.
I guess I never found an explanation for the “Stirling numbers” solution because it would have been just as involved as the solution on my website.
It’s funny. I’ve always wondered about that probability problem ever since I was in school. I’d sometimes ask a teacher or professor about it and they’d say they’d give me an answer after some time. (They NEVER did). 
You would think the answer would be much simpler than that but it sure as heck isn’t.
To help in searching, this problem— finding the distribution of how many of n items chosen uniformly at random with replacement have been seen after t draws— is often called the ‘coupon collector problem’.
Some Call Me… Tim and Itself
Thanks for your replies also.
I guess those 3 replies have given me some hints for further study.
For the sake of completeness, I should probably elaborate a bit more.
The partitions counted by the Stirling number of the second kind don’t differentiate between results. It’s really counting how many ways to partition without caring about the particular partitions at all.
So, for example, let’s say the single die came out “1” eight times, and each of the other faces came up only once. If “4” came up on the 5th roll and “6” came up on the 9th roll, it counts the same for the Stirling number as if “6” came up on the 5th roll and “4” came up on the 9th roll with the other rolls being the same.
But for this particular computation we want to count these as separate results, not as the same result. Hence we multiply by P(6,6), i.e. the number of permutations of the 6 faces.
Repeating what others have said …
Here is the calculation shown at the website, with the division by 6^13 omitted to reduce clutter:
+1 * 6^13 -6 * 5^13 +15 * 4^13 -20 * 3^13 +15 * 2^13 -6 * 1^13
The numbers colored red are the 6th row of Pascal’s triangle, but with alternating signs. (Like most in the thread, I think I memorized the first 6 rows of Pascal’s triangle decades ago, but not the 7th row. Maybe if Backgammon were played with seven-sided dice … 
)
As so often with these problems, the alternating signs of “Inclusion/Exclusion” are key.
If someone tells you that the Stirling partition number {13,6} is 9321312 then you can replace the above calculation with 9321312 * 6!
If, OTOH, we needed to first calculate this Stirling number, we’d end up doing the exact same calculation as above … with one exception: If we didn’t notice the redundancy, we’d divide by 6! per the Stirling formula, and then, for reasons Great Antibob just mentioned, multiply by 6! to solve OP’s problem.
The principle of inclusion/exclusion is covered pretty well above so I’ll throw something else out here even though I don’t necessarily think it’s easier.
Let A be the matrix…
0 1 0 0 0 0 0
0 1/6 5/6 0 0 0 0
0 0 1/3 2/3 0 0 0
0 0 0 1/2 1/2 0 0
0 0 0 0 2/3 1/3 0
0 0 0 0 0 5/6 5/6
0 0 0 0 0 0 1
And let v be the vector…
1 0 0 0 0 0 0
Then v A^n is a row vector that contains the probabilities of having seen k distinct die faces after n rolls of a six sided die.
e.g. v A^13 is…
0 4.59394E-10 9.40609E-06 0.002403777 0.069805693 0.413922929 0.513858194
Here’s one more just for fun…
Let:
A[sub]1[/sub] = 0
A[sub]2[/sub] = 0
A[sub]3[/sub] = 0
A[sub]4[/sub] = 0
A[sub]5[/sub] = 0
A[sub]6[/sub] = 5/324
and
A[sub]n[/sub] = 7/2 A[sub]n-1[/sub] - 175/36 A[sub]n-2[/sub] + 245/72 A[sub]n-3[/sub] - 203/162 A[sub]n-4[/sub] + 49/216 A[sub]n-5[/sub] -5/324 A[sub]n-6[/sub] for n > 6
Then A[sub]n[/sub] gives the probability of having seen all six faces of a die after n rolls.
Thus A[sub]13[/sub] is approximately 0.513858194.
The derivation of this linear recurrence is left as an exercise for the reader.
Great, Sunday morning and I have homework already. Sheesh.
Here’s one more. It’s a little light on why it works, but it is a completely generalizable for ‘dice’ with other than six sides.
Consider the polynomial
p(x) = (x - 1) (x - 2) (x - 3) (x - 4) (x - 5) (x - 6)
= x[sup]6[/sup] - 21 x[sup]5[/sup] + 175 x[sup]4[/sup] - 735 x[sup]3[/sup] + 1624 x[sup]2[/sup] - 1764 x + 720.
Now define a sequence A[sub]n[/sub] recursively so that p(x) is its characteristic polynomial.
A[sub]n[/sub] = 21 A[sub]n-1[/sub] - 175 A[sub]n-2[/sub] + 735 A[sub]n-3[/sub] - 1624 A[sub]n-4[/sub] + 1764 A[sub]n-5[/sub] - 720 A[sub]n-6[/sub]
Using seeds:
A[sub]1[/sub] = 0
A[sub]2[/sub] = 0
A[sub]3[/sub] = 0
A[sub]4[/sub] = 0
A[sub]5[/sub] = 0
A[sub]6[/sub] = 6! = 720
Then A[sub]n[/sub] counts the number of distinct sequences of n numbers chosen from 1 to 6 inclusive where each number from 1 to 6 appears at least once and A[sub]n[/sub] / 6[sup]n[/sup] is the probability of randomly producing a sequence of n numbers chosen from 1 to 6 inclusive where each number from 1 to 6 appears at least once.
Thus A[sub]13[/sub] = 6,711,344,640 and A[sub]13[/sub] / 6[sup]13[/sup]
= 6,711,344,640 / 13,060,694,016
~= 0.513858…