I will try to give (A) a solution on the lines first suggested by lulu short after the question shown up. It uses only elementary steps. We will get an explicit rational number $E$ as expectation after computing a finite sum obtained by simple combinatorial means:
$$
\color{blue}{
E =
\frac{14127973228249375}{380420285792256}
=\frac{5^4 \cdot 11^3 \cdot 16983288629}{2^{31} \cdot 3^{11}}
\approx
37.137801941415212874\dots
}
$$
The formula for the sum given in the sequel and leading to the above fraction explains without computations why the denominator should have that shape, it is a $\Bbb Z$-linear combinations of expressions like
$
\left(\frac 1p-1\right)^{N_1}
\left(\frac 1p-2\right)^{N_2}
\left(\frac 1p-3\right)^{N_3}
\left(\frac 1p-4\right)^{N_4}
$, where $p=\frac15$ is the probability to get each ball in one / each bin, and the involved powers are (specific random) natural numbers.
It is nice and affordable to have (B) a confirmation of the $E$-value from a solution based on different ideas, awkward's solution is indeed awkward, i like it indeed, it's giving also a rational number after computing an integral on $[0,\infty)$, same number, the code is simpler, because this solution is structural. Finally, some simulation part (C) should convince the experimentally thinking reader that the mean is around $37.1$. Computer support will be needed on the road to have a "short answer", i will use sage with hopefully readable code.
(A)
The used modelling probability space $\Omega$ is the space of "paths" which are infinite words
$\omega=w_1w_2w_3\dots$
with letters in the alphabet $A=\{1,2,3,4,5\}$. At time $n$ we see only
the truncated word
$\omega'=w_1w_2w_3\dots w_n$ from $\omega$. Denote by $|\omega'|$ the length of $\omega'$,
which is $n$ in this last sample. Denote by $|\omega'|_k$ the number of letters $=k$ in $\omega'$. So
$|\omega'|=
|\omega'|_1 +
|\omega'|_2 +
|\omega'|_3 +
|\omega'|_4 +
|\omega'|_5$.
The filtration of $\Omega$ is at time $n$ the $\sigma$-algebra generated by the events $E(\omega')$, $E(\omega')$ being the set of all $\omega$ paths starting with that $\omega'$ word of length $n$.
We consider only the paths having the first bin filled first with $5$ balls, then the second bin, then the third one, then the fourth one.
The symmetric group is acting on the four bins to get at least five in time.
So we have to multiply with $4!$.
So we are passing in order through the following "states" (which are events):
$$
\boxed{0\ 0\ 0\ 0\ |\ 0} \overset{(1)}\longrightarrow
\boxed{5\ a\ b\ c\ |\ \#}\overset{(2)}\longrightarrow
\boxed{*\ 5\ d\ f\ |\ \#}\overset{(3)}\longrightarrow
\boxed{*\ *\ 5\ g\ |\ \#}\overset{(4)}\longrightarrow
\boxed{*\ *\ *\ 5\ |\ \#}
\ .
$$
They are each specific unions of $E(\omega')$ cylinders.
The $k$.th component of a state counts the number of balls in the $k$.th bin,
i.e. matches $|\omega'|_k$.
The $\#$ stays for any natural number,
the $*$ for any natural number $\ge 5$, the $5$ for the five, and
for this five we insist that we first get this fifth ball at last in $\omega'$. Explicitly:
$$
\begin{aligned}[]{}
\boxed{0\ 0\ 0\ 0\ |\ 0} &= E(\text{ empty word })=\Omega\ ,
\\
\boxed{5\ a\ b\ c\ |\ \#} &= \bigsqcup_{\omega'}E(\omega')
,\ &&|\omega'|_1=5
,\ |\omega'|_2=a
,\ |\omega'|_3=b
,\ |\omega'|_4=c
\ ;\ \omega'\text{ ends in }1\ ,
\\
\boxed{*\ 5\ d\ f\ |\ \#} &= \bigsqcup_{\omega'}E(\omega')
,\ &&|\omega'|_2=5
,\ |\omega'|_3=d
,\ |\omega'|_4=f
\ ;\ \omega'\text{ ends in }2\ ,
\\
\boxed{*\ *\ 5\ g\ |\ \#} &= \bigsqcup_{\omega'}E(\omega')
,\ &&|\omega'|_3=5
,\ |\omega'|_4=g
\ ;\ \omega'\text{ ends in }3\ ,
\\
\boxed{*\ *\ *\ 5\ |\ \#} &= \bigsqcup_{\omega'}E(\omega')
,\ &&|\omega'|_4=5
\ ;\ \omega'\text{ ends in }4\ .
\end{aligned}
$$
During the passage $(1)$ there are exactly $5$, $a$, $b$, $c$ balls falling respectively in the bins $1,2,3,4$, but there may be some "wasted" balls falling in bin $5$, let $j_1$ be their number. Let $N_1=5+a+b+c$ be the number of "useful" balls for this step.
During the passage $(2)$ there are exactly $a'=5-a$, $b'=d-b$, $c'=f-c$
balls falling respectively in the bins $2,3,4$, but there may be some "wasted" balls falling in bins $1,5$, let $j_2$ be their number. Let $N_2=a'+b'+c'$ be the number of "useful" balls for this step.
During the passage $(3)$ there are exactly $b''=5-d$, $c''=g-f$
balls falling respectively in the bins $3,4$, but there may be some "wasted" balls falling in bins $1,2,5$, let $j_3$ be their number. Let $N_3=b'+c'$ be the number of "useful" balls for this step.
During the passage $(4)$ there are exactly $c'''=5-g$
balls falling respectively in the last needed bin $4$, but there may be some "wasted" balls falling in bins $1,2,3,5$, let $j_4$ be their number. Let $N_4=c'''$ be an other alias the number of "useful" balls for this step.
Then we can write down the formula for the mean value $M$ of steps to get $5$ balls in each bin, in the bin order $1,2,3,4$:
$$
\begin{aligned}
M
&=
\sum_{
\substack{a,b,c;d,f;g\\j_1,j_2,j_3,j_4\ge 0}
}
\binom{N_1-1+j_1}{5-1,a,b,c,j_1}p^{N_1-1}\cdot p^{j_1}\cdot p
\\
&\qquad\qquad\qquad\qquad
\cdot
\binom{N_2-1+j_2}{a'-1,b',c',j_2}p^{N_2-1}\cdot (2p)^{j_2}\cdot p
\\
&\qquad\qquad\qquad\qquad\qquad\qquad
\cdot
\binom{N_3-1+j_3}{b''-1,c'',j_3}p^{N_3-1}\cdot (3p)^{j_3}\cdot p
\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
\cdot
\binom{N_4-1+j_4}{c'''-1,j_4}p^{N_4-1}\cdot (4p)^{j_4}\cdot p
\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad
\cdot
\Big(\underbrace{N_1+N_2+N_3+N_4}_{=20}+j_1+j_2+j_3+j_4\Big)
\\
%
%
%
&=
\sum_{
\substack{a,b,c;d,f;g\\j_1,j_2,j_3,j_4\ge 0}
}
\binom{N_1-1}{5-1,a,b,c}
\binom{N_2-1}{a'-1,b',c'}
\binom{N_3-1}{b''-1,c''}p^{20}
\\
&\qquad\qquad
\cdot
\binom{N_1-1+j_1}{j_1}
\binom{N_2-1+j_2}{j_2}
\binom{N_3-1+j_3}{j_3}
\binom{N_4-1+j_4}{j_4}
\cdot p^{j_1}(2p)^{j_2}(3p)^{j_3}(4p)^{j_4}
\\
&\qquad\qquad\qquad\qquad
\cdot
\Big(20+j_1+j_2+j_3+j_4\Big)
\\
%
%
%
&=
\sum_{a,b,c;d,f;g}
\binom{N_1-1}{5-1,a,b,c}
\binom{N_2-1}{a'-1,b',c'}
\binom{N_3-1}{b''-1,c''}p^{20}
\\
&\qquad\qquad
\cdot
\frac 1{(1-p)^{N_1}}
\cdot
\frac 1{(1-2p)^{N_2}}
\cdot
\frac 1{(1-3p)^{N_3}}
\cdot
\frac 1{(1-4p)^{N_4}}
\\
&\qquad\qquad\qquad\qquad
\cdot
\left(20
+ \frac{N_1\cdot p}{1-p}
+ \frac{N_2\cdot 2p}{1-2p}
+ \frac{N_3\cdot 3p}{1-3p}
+ \frac{N_4\cdot 4p}{1-4p}
\right)
\ .
\end{aligned}
$$
The sum was splitted in pieces corresponding to the terms
$20$, $j_1$, $j_2$, $j_3$, $j_4$. We have used the formula
$$
\sum_{j\ge 0}
j\cdot \binom{N-1+j}{j}q^j
=
\frac {Nq}{(1-q)^{N+1}} \ .
$$
The last sum is finite and can be computed.
To get the mean $E$ of step, considered
without any restriction on the order of the bins
first getting five balls, recalling the action of the symmetric group, we have $E=4!\; M$.
M = 0
p, q1, q2, q3, q4 = 1/5, 4/5, 3/5, 2/5, 1/5
for a, b, c in cartesian_product([[0..4], [0..4], [0..4]]):
N1, C1 = (5+a+b+c), multinomial(5-1, a, b, c)
for d, f in cartesian_product([[b..4], [c..4]]):
N2, C2 = (5+d+f)-(a+b+c), multinomial(5-a-1, d-b, f-c)
for g in [f..4]:
N3, N4, C3 = (5+g)-(d+f), 5-g, multinomial(5-d-1, g-f)
M += C1 * C2 * C3 * p^20 / q1^N1 / q2^N2 / q3^N3 / q4^N4 \
* ( 20 + N1*p/q1 + N2*2*p/q2 + N3*3*p/q3 + N4*4*p/q4 )
E = QQ(24*M)
print(f'Computed value of M is:\nM = {M}')
print(f'Answer to the question is E = 24 M:\nE = {E} = {E.factor()}')
print(f'E ~ {E.n(200)}')
And the above code delivers:
Computed value of M is:
M = 14127973228249375/9130086859014144
Answer to the question is E = 24 M:
E = 14127973228249375/380420285792256 = 2^-31 * 3^-11 * 5^4 * 11^3 * 16983288629
E ~ 37.137801941415212874629304030929071265672012509384861160505
Explicitly:
$$
\color{blue}{
E =
\frac{14127973228249375}{380420285792256}
=\frac{5^4 \cdot 11^3 \cdot 16983288629}{2^{31} \cdot 3^{11}}
\approx
37.137801941415212874\dots
}
$$
(B)
As in the answer of awkward,
and as in the book Analytic Combinatorics, Philippe Flajolet, Robert Sedgewick, page 113 (out of more than 800), II.3 Surjections, Set Partitions, and Sets, Example II.9, formulas (21) and (22),
the usage of exponential generating functions (EGF) is a natural approach.
An early form of the book is / may have been Analytic Combinatorics, Symbolic Combinatorics, Philippe Flajolet, Robert Sedgewick, 2002, page 78, Example 8, Random allocations (balls-in-bins model).
For the convenience of the reader i will cite from either book.
Example. Random allocations (balls-in-bins model). Throw at random $n$ distinguishable
balls into $m$ distinguishable bins. A particular realization is described by a word of length $n$
(balls are distinguishable, say, as numbers from $1$ to $n$) over an alphabet of cardinality m (representing the bins chosen).
Let
$\operatorname{Min}$ and
$\operatorname{Max}$ represent the size of the least filled and most filled
bins, respectively. Then,
$$
\tag{$21$}
$$
$$
\begin{aligned}
\mathbb P\{\ \operatorname{Max} \le b\ \}
&= n! \ [z^n]\ e_b\left(\frac zm\right)^m
\\
\mathbb P\{\ \operatorname{Min} > b\ \}
&= n! \ [z^n]\ \left(\exp\frac zm -e_b\left(\frac zm\right)\right)^m
\ .
\end{aligned}
$$
The justification of this formula relies on the easy identity
$$
\tag{$22$}
\frac 1{m^n} [z^n] f(z) \equiv [z^n] f\left(\frac zm\right)\ ,
$$
and on the fact that a probability is determined as the ratio between the number of favorable
cases (given by $(19)$) and the total number of cases ($m^n$).
Here, $e_b(z)$ is the truncated version of $\exp z$, the Taylor expansion of $\exp$ around $z=0$ stopping in degree $b$. The operator $[z^n]$ isolates from an analytic function $f$ the piece in $z^n$ in its Taylor series around zero.
In our problem, we have to change slightly the second line in $(21)$, we use four times the factor $\left(\exp\frac zm -e_b\left(\frac zm\right)\right)$ with $b=4$, $m=5$ (strictly more than $4$ balls in the first four bins) and one more factor $\exp\frac zm=\left(\exp\frac zm -e_{-1}\left(\frac zm\right)\right)$
(strictly more than $-1$ balls in the last bin), which is $(*)$ from awkward's answer:
$$
p_n:=\mathbb P[T\le n]
=n!\; [z^n]\ \left(\exp\frac z5 -e_4\left(\frac z5\right)\right)^4\exp\frac z5
\ ,
$$
where $T$ is the random variable whose mean $E=\Bbb E[T]$ is wanted. The coefficients $(p_n)$ fished in running degree $n$ from the above analytic function converge increasingly to one, $\nearrow 1$. To get the mean we build $q_n=1-p_n\searrow 0$, and add them. As in awkard's answer:
$$
\begin{aligned}
E&=\Bbb E[T]=\sum q_n
\\
&=
\sum_{n\ge 0}
n!\; [z^n]\ \exp z-\left(\exp\frac z5 -e_4\left(\frac z5\right)\right)^4\exp\frac z5
\\
&=
\int_0^\infty
\exp (-z)\left[
\exp z-\left(\exp\frac z5 -e_4\left(\frac z5\right)\right)^4\exp\frac z5
\right]\; dz\\
&=
\int_0^\infty
\left[
1-\left(1-\exp\left(-\frac z5\right)e_4\left(\frac z5\right)\right)^4
\right]\; dz
\\
&=
\int_0^\infty
\left[
1-\left(1-e^{-y}\left(1+y+\frac{y^2}2+\frac{y^3}6+\frac{y^4}{24}\right)\right)^4
\right]\; 5\;dy
\\
&=
5\int_0^\infty
4\cdot e^{-y}\left(1+y+\frac{y^2}2+\frac{y^3}6+\frac{y^4}{24}\right)\;dy
\\
&\qquad
-
5\int_0^\infty
6\cdot e^{-2y}\left(1+y+\frac{y^2}2+\frac{y^3}6+\frac{y^4}{24}\right)^2\;dy
\\
&\qquad
+
5\int_0^\infty
4\cdot e^{-3y}\left(1+y+\frac{y^2}2+\frac{y^3}6+\frac{y^4}{24}\right)^3\;dy
\\
&\qquad
-
5\int_0^\infty
1\cdot e^{-4y}\left(1+y+\frac{y^2}2+\frac{y^3}6+\frac{y^4}{24}\right)^4\;dy
\ .
\end{aligned}
$$
The last expression can be computed also manually,
by computing some powers of some polynomials, then replacing $y$-powers by factorials, using sage:
var('y');
e4 = 1 + y + y^2/2 + y^3/6 + y^4/24
f4 = exp(-y) * e4
E = 5*integrate( 4*f4 - 6*f4^2 + 4*f4^3 - f4^4, y, 0, oo )
E = QQ(E)
print(f'E = {E}\n = {E.factor()}\n ~ {E.n()}')
And we get:
E = 14127973228249375/380420285792256
= 2^-31 * 3^-11 * 5^4 * 11^3 * 16983288629
~ 37.1378019414152
(C) Simulation.
import numpy as np
r = np.random.default_rng(int(1234567890)) # randomizer
E = 0.0
N = 10*6 # trials
for trial in range(N):
a = r.integers(low=1, high=6, size=150) # random array with entries among 1,2,3,4,5
try:
T = 1 + max( [ int( np.argwhere( (a-2)(a-3)(a-4)(a-5) )[4] ),
int( np.argwhere( (a-1) (a-3)(a-4)(a-5) )[4] ),
int( np.argwhere( (a-1)(a-2) (a-4)(a-5) )[4] ),
int( np.argwhere( (a-1)(a-2)(a-3) *(a-5) )[4] ), ] )
E += T / N
except:
pass
print(E)
This time i've got:
37.13368599999
Here, a is an array of size $150$ with entries among $1,2,3,4,5$, and for instance
a - 3 is the array obtained from a by subtracting $3$ from each component.
The product (a-2)*(a-3)*(a-4)*(a-5) is also built componentwise, and then
np.argwhere fishes the positions with non-zero values.
So it fishes the positions of the 1.
The pythonically fourth, humanly fifth
position is then the index of the humanly fifth occurrence of 1 in the list a.
This index is also pythonic, so we need to add one to get its human version.
One can do better, but the above lazy code for the simulation is simpler to explain, and
should be easier to digest.
