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;

It is also possible to put several of these graphs together in a 3-D plot.

> n := 'n':
t := 't':
eta:=(t,n) -> product(k/n*t/(1-(1-k/n)*t), k=1..n):
for n from 2 to 12 do
series(eta(t,n), t, 51):
Prob[n] := [ seq(coeff(%, t, i), i = 1..50)]:
pmf[n] := zip((x,y)->(x,y),domain,Prob[n]):
od:
yellowsides := proc(x, y, z, u)
POLYGONS([[x,y,0], [x+u, y, 0], [x+u, y, z], [x,y,z]],
[[x,y,0], [x,y+u,0], [x,y+u,z], [x,y,z]],
[[x + u, y, 0], [x+u, y+u, 0], [x+u, y+u, z], [x+u,y,z]],
[[x+u,y+u,0], [x,y+u,0], [x,y+u,z], [x+u,y+u,z]],
COLOR(RGB,1,1,0));
end:
greentop := proc(x,y,z,u)
POLYGONS([[x,y,z], [x+u,y,z], [x+u, y+u,z], [x,y+u,z]],
COLOR(RGB,0,1,0));
end:
sides := seq(seq(yellowsides(i-1/2, j - 1/2, Prob[i][j],1), i = 3 .. 12), j = 1 .. 50):
tops := seq(seq(greentop(i-1/2, j - 1/2, Prob[i][j],1), i = 3 .. 12), j = 1 .. 50):

> P3d := PLOT3D(sides, tops, STYLE(PATCH), AXESSTYLE(BOXED), ORIENTATION(24, 66)):
display3d({P3d}, title=`Number of rolls needed to observe each face of\nan n-sided die at least once, n = 3, 4, ..., 12`, titlefont=[TIMES,BOLD,14]);