CDF of probability distribution with replacement [duplicate]

I want to get every color of gumball in a gumball machine (where there are 16 types of gumballs, each with a 1/16 chance of obtaining a particular color [assume there are an infinite amount of gumballs]).

I'm interested in knowing what the PMF (as a function of # of attempts) of getting every single gumball. ie. what is the probability of getting all 16 colors of gumballs in n trials?

I'd rather not have to brute-force the answer and am hoping for a somewhat elegant solution (maybe we can use the negative binomial distribution to find this?).

Let $c=16$ be the number of colours.

Imagine the position immediately before the $n$th gumball. The probability of having exactly $c-1$ colours covered is the ratio of the number of ways of choosing the $c-1$ colours, namely ${c \choose c-1}=c$, multiplied by the number of arrangements giving a positive number of each of the $c-1$ colours from the $n-1$ gumballs, namely $S_2(n-1,c-1)$ using Stirling numbers of the second kind since the order of the chosen colurs matters, multiplied by the possible orders of the chosen colours since they are distinguishable, namely $(c-1)!$, all divided by the full number of possibilities for $n-1$ balls and $c$ colours, namely $c^{n-1}$, which gives $\frac{c \times S_2(n-1,c-1) \times (c-1)! }{c^{n-1}} = \frac{c!}{c^{n-1}} S_2(n-1,c-1).$

The conditional probability that the $n$th ball is then the final colour is $\frac1c$, making the overall probability that the $n$th gumball completes the set $$\dfrac{c!}{c^n}S_2(n-1,c-1).$$

For $c=16$, this suggests the following approximate probabilities. Using unrounded figures confirms that this calculation gives an expected value of $16H_{16}\approx 54.09$ as is known for the coupon collector's problem based on the harmonic numbers, and that the distribution is slightly right-skewed.

n   Probability
1   0.000000
2   0.000000
3   0.000000
4   0.000000
5   0.000000
6   0.000000
7   0.000000
8   0.000000
9   0.000000
10  0.000000
11  0.000000
12  0.000000
13  0.000000
14  0.000000
15  0.000000
16  0.000001
17  0.000009
18  0.000035
19  0.000102
20  0.000241
21  0.000489
22  0.000885
23  0.001460
24  0.002239
25  0.003232
26  0.004435
27  0.005832
28  0.007393
29  0.009082
30  0.010856
31  0.012671
32  0.014482
33  0.016249
34  0.017933
35  0.019505
36  0.020939
37  0.022216
38  0.023324
39  0.024256
40  0.025009
41  0.025586
42  0.025992
43  0.026236
44  0.026328
45  0.026280
46  0.026105
47  0.025816
48  0.025427
49  0.024950
50  0.024398
51  0.023784
52  0.023118
53  0.022411
54  0.021672
55  0.020910
56  0.020133
57  0.019348
58  0.018560
59  0.017775
60  0.016998
61  0.016231
62  0.015479
63  0.014744
64  0.014029
65  0.013334
66  0.012661
67  0.012012
68  0.011386
69  0.010785
70  0.010209
71  0.009656
72  0.009128
73  0.008624
74  0.008144
75  0.007686
76  0.007251
77  0.006837
78  0.006445
79  0.006073
80  0.005720
81  0.005386
82  0.005070
83  0.004772
84  0.004489
85  0.004223
86  0.003971
87  0.003734
88  0.003510
89  0.003299
90  0.003100
91  0.002912
92  0.002736
93  0.002570
94  0.002413
95  0.002266
96  0.002128
97  0.001998
98  0.001875
99  0.001760
100 0.001652
101 0.001551
102 0.001455
103 0.001366
104 0.001281
105 0.001202
106 0.001128
107 0.001058
108 0.000993
109 0.000931
110 0.000874
111 0.000819
112 0.000769
113 0.000721
114 0.000676
115 0.000634
116 0.000595
117 0.000558
118 0.000523
119 0.000491
120 0.000460
121 0.000431
122 0.000405
123 0.000379
124 0.000356
125 0.000334
126 0.000313
127 0.000293
128 0.000275
129 0.000258
130 0.000242
131 0.000227
132 0.000213
133 0.000199
134 0.000187
135 0.000175
136 0.000164
137 0.000154
138 0.000144
139 0.000135
140 0.000127
141 0.000119
142 0.000112
143 0.000105
144 0.000098
145 0.000092
146 0.000086
147 0.000081
148 0.000076
149 0.000071
150 0.000067
151 0.000062
152 0.000059
153 0.000055
154 0.000051
155 0.000048
156 0.000045
157 0.000042
158 0.000040
159 0.000037
160 0.000035
161 0.000033
162 0.000031
163 0.000029
164 0.000027
165 0.000025
166 0.000024
167 0.000022
168 0.000021
169 0.000020
170 0.000018
171 0.000017
172 0.000016
173 0.000015
174 0.000014
175 0.000013
176 0.000012
177 0.000012
178 0.000011
179 0.000010
180 0.000010
181 0.000009
182 0.000008
183 0.000008
184 0.000007
185 0.000007
186 0.000007
187 0.000006
188 0.000006
189 0.000005
190 0.000005
191 0.000005
192 0.000004
193 0.000004
194 0.000004
195 0.000004
196 0.000003
197 0.000003
198 0.000003
199 0.000003
200 0.000003
201 0.000002
202 0.000002
203 0.000002
204 0.000002
205 0.000002
206 0.000002
207 0.000002
208 0.000002
209 0.000001
210 0.000001
211 0.000001
212 0.000001
213 0.000001
214 0.000001
215 0.000001
216 0.000001
217 0.000001
218 0.000001
219 0.000001
220 0.000001
221 0.000001
222 0.000001
223 0.000001
224 0.000001
225 0.000001

Here's a much simpler approach that I just thought of:

Note that $P(A)=1-P(\neg A)$ In our case, lets define $\neg A$ to be event where at least one color $C_i$is missing $\neg A :=\{ \neg C_1 \cup \neg C_2....\cup \neg C_N\}$

What is $P(\neg A)$?

Lets say you draw $N$ gumballs $G_1...G_N$. The the probability that you dont get one of a particular color in N tosses can be modeled as Binomial:

$$P(G_1 \neq C_i,...,G_N \neq C_i)=\left(\frac{15}{16}\right)^{N}$$

However, there are $16$ colors, so multiply this by $16$ and then correct for the fact that more than one color can be missing to get:

$$P(\neg A)=16\left(\frac{15}{16}\right)^{N}-P(I)$$ where $I$ is calculated as a sum of intersection probabilities as per the inclusion-exclusion principle. For 16 colors, this will be a long formula, but mathematically it is rather simple: $P(\neg C_1 \cap \neg C_2)=\left(\frac{14}{16}\right)^{N}$ since there are now two colors that must be avoided.