I believe that it is quite instructive to examine this problem using
generating functions and combinatorial classes as presented in
Analytic Combinatorics by Flajolet and Sedgewick.
Using the notation from the accepted answer we have for partitions of
cell size at least $\gamma$ and the maximum element exactly $\mu$
the combinatorial class
$$\def\textsc#1{\dosc#1\csod}
\def\dosc#1#2\csod{{\rm #1{\small #2}}}
\textsc{SET}(\textsc{SET}_{=\gamma}(\mathcal{Z}))
\times \textsc{SET}(\textsc{SET}_{=\gamma+1}(\mathcal{Z}))
\times \textsc{SET}(\textsc{SET}_{=\gamma+2}(\mathcal{Z}))
\\ \times \textsc{SET}(\textsc{SET}_{=\gamma+3}(\mathcal{Z}))
\times \cdots
\times \textsc{SET}(\textsc{SET}_{=\mu-1}(\mathcal{Z}))
\times \textsc{SET}_{\ge 1}(\textsc{SET}_{=\mu}(\mathcal{Z})).$$
This immediately produces the generating function which is
$$G(z) = \left(\exp\left(\frac{z^\mu}{\mu!}\right)-1\right)
\prod_{q=\gamma}^{\mu-1} \exp\left(\frac{z^q}{q!}\right).$$
This is (here we recognize the class
$\textsc{SET}(
\textsc{SET}_{\gamma\le\cdot\le\mu-1}(\mathcal{Z}))$)
$$G(z) = \left(\exp\left(\frac{z^\mu}{\mu!}\right)-1\right)
\exp\left(\sum_{q=\gamma}^{\mu-1} \frac{z^q}{q!}\right).$$
We can extract an closed formula from this, getting
$$n! [z^n] G(z)
= n! \sum_{k=1}^{\lfloor n/\mu \rfloor}
\frac{1}{(\mu!)^k k!}
[z^{n-k\mu}]
\exp\left(\sum_{q=\gamma}^{\mu-1} \frac{z^q}{q!}\right).$$
Now let $\lambda\vdash_\gamma^{\mu-1} (n-k\mu)$ be a partition
of $n-k\mu$ with parts at least $\gamma$ and at most $\mu-1$
and let the multiplicities of $\lambda$ be denoted by
$1^{f_1} 2^{f_2} 3^{f_3} \ldots$ to obtain
$$n! \sum_{k=1}^{\lfloor n/\mu \rfloor}
\frac{1}{(\mu!)^k k!}
\sum_{\lambda\vdash_\gamma^{\mu-1} (n-k\mu)}
\prod_{q=\gamma}^{\mu-1} \frac{1}{(q!)^{f_q} (f_q)!}.$$
Observe that this formula produces the same values as the recurrence
from the accepted answer, as can be verified using the following Maple
code:
with(combinat);
B :=
proc(g, m, n)
option remember;
if n < m then return 0 fi;
if n = m then return 1 fi;
add(binomial(n-1,k-1)*B(g, m, n-k), k=g..m-1)
+ add(binomial(n-1, m-1)*B(g, j, n-m), j=g..m);
end;
G :=
proc(g, m, n)
option remember;
local gf;
gf := (exp(z^m/m!)-1)
*exp(add(z^q/q!, q=g..m-1));
n!*coeftayl(gf, z=0, n);
end;
Q :=
proc(g, m, n)
option remember;
local p, k, res, interm, f;
res := 0;
for k to floor(n/m) do
interm := 0;
p := firstpart(n-m*k);
while type(p, list) do
if min(p) >= g and max(p) <= m-1 then
f := convert(p, 'multiset');
interm := interm +
1/mul((ff[1]!)^ff[2]*ff[2]!, ff in f);
fi;
p := nextpart(p);
od;
res := res + interm/(m!)^k/k!;
od;
n!*res;
end;