This solution was written by Joshua Levy during the summer of 1994 when he was a Hope College student.

First find the probability generating function for the geometric distribution and then enter the probability generating function of X , the number of rolls needed for an n -sided die.

> f := (1-p)^(x-1)*p;
eta := sum(t^x*f, x = 1 .. infinity);

Note that for a 6-sided die, for example, that the respective probabilities for observing a new face, given that k faces, k = 0, 1, ... , 5, have been observed, are 6/6, 5/6, 4/6, 3/6, 2/6, 1/6. Thus, in general, we can write the probability generating function, when rolling an n -sided die as follows:

> n := 'n':
t := 't':
eta:=(t,n) -> product(k/n*t/(1-(1-k/n)*t), k=1..n); #Written reverse order

Maple can simplify this expression .

> k := 'k':
eta(t,n);
eta(t,n) := simplify(%);

To find the mean of X symbolically , evaluate the first derivative at t = 1.

> etap := diff(%, t):
mu := subs(t=1, %);
mu := simplify(%);

To find particular values of the mean, e.g., when n = 6 or n = 20 , substitute these values for n .

> simplify(subs(n=6, mu)), simplify(subs(n=20, mu)); evalf(%);

To find the variance of X symbolically, use the second derivative of the probability-generating function .

> diff(etap, t):
subs(t=1, %) + mu - mu^2:
var := simplify(%);

Now find particular values of the variance, for example, when n = 6 or n = 20 .

> simplify(subs(n=6, var)), simplify(subs(n=20, var)); evalf(%);

We now look at some expansions of probability-generating functions. By "stripping" off the respective coefficients of , we can display the respective probabilities. Recall that the second argument gives the value of n , the number of sides on the die .

> eta(t,1);

> eta(t,2);
series(eta(t,2), t, 10);

> eta(t,3);
series(eta(t,3), t, 16);

> Prob[3] := [ seq(coeff(%, t, i), i = 1..15) ];

In order to " see " what this p.m.f. looks like, we shall plot a probability histogram.

> domain := [k$k=1..15];
pmf[3] := zip((x,y)->(x,y),domain,Prob[3]);

> k:='k':
xtics := [k$k=1..15]:
plot3 := ProbHistFill(pmf[3]):
ppp := plot([[0,0],[0,0]], x = 0 .. 15.7, xtickmarks=xtics, labels=[``,``]):
t3 := textplot([10, 0.2, `n = 3`],font=[TIMES,BOLD,14]):
display({plot3,t3,ppp}, title=`Number of Rolls Needed to Observe\nEach Face at Least Once, 3-sided Dice`, titlefont=[TIMES,BOLD,14]);

Here is an animation for n = 4, 5, ..., 12 .

> k := 'k':
xtics := [2*k$k=1..35]:
domain := [k$k=1..70]:
ytics := [k*0.01$k=1..15]:

> for n from 4 to 12 do
series(eta(t,n), t, 71):
Prob[n] := [ seq(coeff(%, t, i), i = 1..70)]:
pmf[n] := zip((x,y)->(x,y),domain,Prob[n]):
od:

> for k from 4 to 12 do
plotpmf[k] := ProbHistFill(pmf[k]):
plotscale := plot([[0,0], [0,0]], x = 0 .. 70.2, y = 0 .. 0.152, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]):
n := k:
txt[k] := textplot([50, 0.12, cat(`n = `,convert(n,string))], font=[TIMES,BOLD,12]):
mean := evalf(simplify(subs(n=j, mu))):
txtm[k] := textplot([50, 0.105, cat(`m = `,convert(mean,string))], font=[SYMBOL,12]):
SD := evalf(sqrt(simplify(subs(n=j, var)))):
txtsd[k] := textplot([50, 0.09, cat(`s = `,convert(SD,string))], font=[SYMBOL,12]):
H[k] := display({plotpmf[k], plotscale, txt[k], txtm[k], txtsd[k]}, title=`Number of rolls needed, n-sided die`, titlefont=[TIMES,BOLD,14]):
od:

> display([seq(H[k], k = 4 .. 12)],insequence=true, title = `The number of rolls needed to observe each face\nof an n-sided die at least once`, titlefont = [TIMES,BOLD,14]);

> k := 'k':
xtics := [5*k$k=1..30]:
domain := [k$k=1..150]:
ytics := [k*0.005$k=1..7]:

> series(eta(t,20), t, 151):
Prob[20] := [ seq(coeff(%, t, i), i = 1..150)]:
pmf[20] := zip((x,y)->(x,y),domain,Prob[20]):

> plotpmf[20] := ProbHistFill(pmf[20]):
plotscale := plot([[0,0], [0,0]], x = 0 .. 150.2, y = 0 .. 0.022, xtickmarks=xtics, ytickmarks=ytics, labels=[``,``]):
n := 20:
txt[20] := textplot([120, 0.02, cat(`n = `,convert(n,string))], font=[TIMES,BOLD,12]):
mean := evalf(simplify(subs(n=20, mu))):
txtm[20] := textplot([120, 0.018, cat(`m = `,convert(mean,string))], font=[SYMBOL,12]):
SD := evalf(sqrt(simplify(subs(n=20, var)))):
txtsd[20] := textplot([120, 0.016, cat(`s = `,convert(SD,string))], font=[SYMBOL,12]):
display({plotpmf[20], plotscale, txt[20], txtm[20], txtsd[20]}, title=`Number of rolls needed to observe each face\nof a 20-sided die at least once`, titlefont=[TIMES,BOLD,14]);

Suppose you want to collect each of 20 different English porcelain miniatures that are put into a 100-bag box of Red Rose tea and you use one tea bag per day. On the average, approximately how long will it take to complete the collection, in years?

> years := (mean*100)/365;