Can this probability function speed be improved?

I have a process that as part of it requires a certain probability query, specifically, given a multinomial distribution with f equiprobable categories and r trials, some specific category has a given value v and that value is the (not necessarily unique) maximum across all categories.

This can obviously be done directly in Mathematica with

Probability[a[1] == v && v == Max[Array[a, f]], 
 Array[a, f] [Distributed] MultinomialDistribution[r, ConstantArray[1/f, f]]]

but it is quite sluggish, for example with {r,f,v}={25,12,6} it takes some… well, I don’t know, I aborted it after several minutes running.

I’ve come up with

f[r_, f_, v_] := 
  Block[{t}, 
   Coefficient[
    Expand[Sum[t^j/j!, {j, 0, v}]^(f - 1)] (f - 1)^(-(r - v))*
           (r - v)!*Binomial[r, v]/f^v (1 - 1/f)^(r - v), t^(r - v)]];

which is fairly snappy.

Can this query be done faster? I’d like to handle cases of r up to ~200, f up to ~25 and v generally ~r/3.

Edit/Update:

Using a technique from an old answer of mine, the new function

f2[r_, f_, v_, pf_ : Infinity] := Module[{pMax},
   pMax[q_, b_, max_, p_] := Module[{h}, h[0, n_, m_] = N[1, p];
     h[s_, n_, m_] := 
      h[s, n, m] = 
       Sum[(n*x + x - s) (h[s - x, n, m]/x!), {x, Min[s, m]}]/s;
     q! h[q, b, max]/b^q];
   pMax[r - v, f - 1, v, pf] Binomial[r, v]/f^v (1 - 1/f)^(r - v)];

is notably faster on larger cases and comparable on small ones. For example, with {r,f,v}={25,12,15} this takes ~0.36ms, vs ~2.5ms for my original function. By way of comparison, the direct calculation in Mathematica took ~2800 seconds.

If no better result arrives, I’ll self-answer with this function.

f[r_, f_, v_] := Block[{t}, Coefficient[Expand[Sum[t^j/j!, {j, 0, v}]^(f - 1)] (f - 1)^(-(r - v))*
  (r - v)!*Binomial[r, v]/f^v (1 - 1/f)^(r - v), t^(r - v)]];

f2[r_, f_, v_, pf_ : Infinity] := Module[{pMax}, 
   pMax[q_, b_, max_, p_] := Module[{h}, h[0, n_, m_] = N[1, p];
     h[s_, n_, m_] := h[s, n, m] = Sum[(n*x + x - s) (h[s - x, n, m]/x!), {x, Min[s, m]}]/s;
     q! h[q, b, max]/b^q];
   pMax[r - v, f - 1, v, pf] Binomial[r, v]/f^v (1 - 1/f)^(r - v)];

(* This function only useful when 2v > r *)
f3[r_, f_, v_] := If[2 v > r, Binomial[r, v] (f - 1)^(r - v)/f^r]

r = 200;
g = 25;  (* Really this is f but I get confused between the constant f and the function f *)
data = ConstantArray[{0, 0, 0, 0}, r];
Do[data[[v]] = {v, If[2 v > r, AbsoluteTiming[Binomial[r, v] (g - 1)^(r - v)/g^v][[1]], 10],
   AbsoluteTiming[f[r, g, v]][[1]],
   AbsoluteTiming[f2[r, g, v]][[1]]}, {v, 1, r - 1}]

ListLogPlot[{data[[All, {1, 3}]], data[[All, {1, 4}]], data[[All, {1, 2}]]},
 PlotRange -> {Automatic, {0, 9}}, Frame -> True,
 FrameLabel -> (Style[#, Bold, 18] &) /@ {"v", "Seconds"},
 PlotLegends -> {"f", "f2", "f3"}]