As evidenced by some of my previous answers, I like to write quick numerical simulations if they seem feasible. Bingo seems especially easy (Python code below).
I'm not sure if this is true, but I think the Bingo cards are essentially independent of each other. That is, if we can compute the probability distribution of a single player $N=1$ game length, we can use that to compute the joint probabilities for any number of players.
What I get seems to match with your playing experience, the mean game length for a single player was $42.4$ with a standard deviation of $9.6$. There is a slight skew in the PDF towards longer games. The full PDF is shown below:
from numpy import *
from collections import Counter
def new_board():
cols = arange(1,76).reshape(5,15)
return array([random.permutation(c)[:5] for c in cols])
def new_game():
for token in random.permutation(arange(1,76)):
yield token
def winning(B):
if (B.sum(axis=0)==5).any(): return True
if (B.sum(axis=1)==5).any(): return True
if trace(B)==5 or trace(B.T)==5: return True
return False
def game_length(board, game):
B = zeros((5,5),dtype=bool)
B[2,2] = True
for n,idx in enumerate(game):
if winning(B): return n
B[board==idx] = True
def simulation(trials):
C = Counter()
b = new_board()
for _ in xrange(trials):
C[game_length(b, new_game())] += 1
return C
repeats = 10**2
trials = 10**3
from multiprocessing import *
P = Pool()
sol = sum(P.map(simulation,[trials,]*repeats))
P.close()
P.join()
X = array(sorted(sol.keys()))
Y = array([float(sol[x]) for x in X])
Y/= repeats*trials
EX = array(list(sol.elements()))
print "Mean and stddev", EX.mean(), EX.std()
import pylab as plt
plt.fill_between(X, Y, lw=2, alpha=.8)
plt.plot([EX.mean(),EX.mean()], [0,1.2*max(Y)], 'r--',lw=2)
plt.ylim(ymax = 1.2*max(Y))
plt.xlabel("Expected game length")
plt.show()
Let $N^{(k)}_{a,b,c}$ be the number of possible unnumbered layouts of the last $k$ columns, given that there are $a$ numbers left to assign in the first row, $b$ in the second row, and $c$ in the third row. Given the rules (at least one number per column), you have the recursion
$$
\begin{eqnarray}
N^{(k)}_{a,b,c}&=&N^{(k-1)}_{a-1,b-1,c-1}+N^{(k-1)}_{a-1,b-1,c}+N^{(k-1)}_{a-1,b,c-1}+N^{(k-1)}_{a,b-1,c-1} + N^{(k-1)}_{a-1,b,c}+N^{(k-1)}_{a,b-1,c}+N^{(k-1)}_{a,b,c-1},
\end{eqnarray}
$$
with the boundary conditions that $N^{(0)}_{0,0,0}=1$ and $N^{(k)}_{a,b,c}=0$ unless all indices are between $0$ and $k$. You want to find $N^{(9)}_{5,5,5}$. The following Python function performs the calculation:
def N(a,b,c,k,cache={(0,0,0,0):1}):
if min(a,b,c)<0 or max(a,b,c)>k:
return 0
if not cache.has_key((a,b,c,k)):
val = N(a-1,b-1,c-1,k-1)
val += N(a-1,b-1,c,k-1) + N(a-1,b,c-1,k-1) + N(a,b-1,c-1,k-1)
val += N(a-1,b,c,k-1) + N(a,b-1,c,k-1) + N(a,b,c-1,k-1)
cache[(a,b,c,k)] = val
return cache[(a,b,c,k)]
N(5,5,5,9)
> 735210
The result is below the upper bound of ${{9}\choose{5}}^3=126^3=2000376$ obtained by allowing all-blank columns, as it should be.
For numbered layouts, the recursion is slightly different, because the number of ways to choose the numbers depends on the column. Here, let $m_k$ be the number of allowed numbers in the $k$-th column from the end: $m_1=11$, $m_9=9$, and $m_k=10$ otherwise. In a column with no blanks, there are ${m_k}\choose{3}$ ways to choose the numbers; with one blank, ${m_k}\choose{2}$; and with two blanks, $m_k$. The recursion is therefore
$$
\begin{eqnarray}
M^{(k)}_{a,b,c}&=&{{m_k}\choose{3}}M^{(k-1)}_{a-1,b-1,c-1}+ {{m_k}\choose{2}}\left(M^{(k-1)}_{a-1,b-1,c}+M^{(k-1)}_{a-1,b,c-1}+M^{(k-1)}_{a,b-1,c-1}\right) \\ && + m_k\left(M^{(k-1)}_{a-1,b,c}+M^{(k-1)}_{a,b-1,c}+M^{(k-1)}_{a,b,c-1}\right).
\end{eqnarray}
$$
Not unexpectedly, the result here is much larger:
$$
M^{(9)}_{5,5,5} = 3669688706217187500.
$$
As a sanity check on this result, one might consider the average branching factor for the $15$ numbers in an average unnumbered layout:
$$
\left(\frac{M^{(9)}_{5,5,5}}{N^{(9)}_{5,5,5}}\right)^{\frac{1}{15}} \approx 7.023.
$$
Since in a typical row (with $m_k=10$) the possible outcomes have average branching factors ${{10}\choose{3}}^{1/3}\approx 4.93$, ${{10}\choose{2}}^{1/2}\approx 6.71$, and $10$, this seems very reasonable.
Note that we can also enumerate the possible unnumbered layouts in an entirely different way, using inclusion-exclusion. The basic idea is to count all the ways of placing $5$ numbers in each of the three rows, then remove the cases with an all-blank column, then add back in the cases with two all-blank columns, and so on:
$$
\begin{eqnarray}
N^{(9)}_{5,5,5} &=& {{9}\choose{5}}^3 - {{9}\choose{1}}{{8}\choose{5}}^3 + {{9}\choose{2}}{{7}\choose{5}}^3 - {{9}\choose{3}}{{6}\choose{5}}^3 + {{9}\choose{4}} \\
&=& 735210,
\end{eqnarray}
$$
as above. This is an elegant solution to the unnumbered case, but I do not see a way to extend it to the numbered layouts.
Best Answer
The probability $P(A)$ of hitting the first row (all 5 numbers) is $$P(A)=\frac{\binom{55}{25}}{\binom{60}{30}}$$ Given $A$, the probability $P(B)$ of hitting the other 4 numbers in the first and last columns is $$P(B)=\frac{\binom{51}{21}}{\binom{55}{25}}$$ You want $P(A)\times(1-P(B))$ $$\frac{\binom{55}{25}}{\binom{60}{30}}\times\frac{\binom{55}{25}-\binom{51}{21}}{\binom{55}{25}}=\frac{\binom{55}{25}-\binom{51}{21}}{\binom{60}{30}}$$