Probability for finding a double-as-likely event Repeating an experiment with $n$ possible outcomes $t$ times independently, where all but one outcomes have probability $\frac{1}{n+1}$ and the other outcome has the double probability $\frac{2}{n+1}$, is there a good approximate formula for the probability that the outcome with the higher probability happens more often than any other one?
For me, $n$ is typically some hundreds, and $t$ is chosen depending on $n$ such that the probability that the most likely outcome occurs most often is between 10% and 99.999%.
In the moment I use a small program that calculates a crude approximation by assuming that the counts for how often each outcome shows up in $t$ trials are independent and approximate the counts using the Poisson distribution. How can I improve on this?
EDIT: I'd strongly appreciate comments/votes on the two (maybe soon more) answers given.
EDIT 2: As none of the two answers is convincing me, but as I don't want to let the 100 points bounty to vanish (and as nobody voted for/against one of the two answers), I'll just pick one of the answers. I'd still appreciate other answers.
 A: Partition the outcomes by the frequency of occurrences $x$ of the "double outcome", $0 \le x \le t$.  Conditional on this number, the distribution of the remaining $t-x$ outcomes is multinomial across $n-1$ equiprobable bins.  Let $p(t-x, n-1, x)$ be the chance that no bin out of $n-1$ equally likely ones receives more than $x$ outcomes.  The sought-for probability therefore equals
$$\sum_{x=0}^{t} \binom{t}{x}\left(\frac{2}{n+1}\right)^x \left(\frac{n-1}{n+1}\right)^{t-x} p(t-x,n-1,x).$$
In Exact Tail Probabilities and Percentiles of the Multinomial Maximum, Anirban DasGupta points out (after correcting typographical errors) that $p(n,K,x)K^n/n!$ equals the coefficient of $\lambda^n$ in the expansion of $\left(\sum_{j=0}^{x}\lambda^j/j!\right)^K$ (using his notation).  For the values of $t$ and $n$ involved here, this coefficient can be computed in at most a few seconds (making sure to discard all $O(\lambda^{n+1})$ terms while performing the successive convolutions needed to obtain the $K^{\text{th}}$ power).  (I checked the timing and corrected the typos by reproducing DasGupta's Table 4, which displays the complementary probabilities $1 - p(n,K,x)$, and extending it to values where $n$ and $K$ are both in the hundreds.)
Quoting a theorem of Kolchin et al., DasGupta provides an approximation for the computationally intensive case where $t$ is substantially larger than $n$.  Between the exact computation and the approximation, it looks like all possibilities are covered.
A: I agree with some comments, in that the Poisson approximation sounds nice here (not a 'crude' approximation). It should be asympotically exact, and it seems the most reasonable thing to do, as an exact analytic solutions seems difficult.
As an intermediate alternative, (if you really need it) I suggest I fisrt order correction to the Poisson approximation, in the following way (I've done something similar some time ago, and it worked).
As suggested by a comment, your model is (not approximately but exactly) Poisson if we condition on the sum. That is:
Let $X_t$ ($t$ is a parameter here) be a vector of $n$ independent Poisson variables, the first one with $\lambda = 2t/(n+1)$, the others with $\lambda = t/(n+1)$. Let $s=\sum x$, so $E(s)=t$. It is clear that $X_t$ is not equivalent to other model (because our model is restricted to $s=t$), but it is a good approximation. Further, the distribution of $X_t | s$ is equivalent to our model. Indeed, we can write
$ \displaystyle P(X_t) = \sum_s P(X_t | s) P(s)$ 
This can also be writen for the event in consideration (that $x_1 $ is the maximum).
We know to compute the LHS, and $P(s)$, but we are interested in the other term. Our first order Poisson approximation comes from assuming that $P(s)$ concentrates about the mean so that it can be assimilated to a delta, and then $ P(X_t) \approx  P(X_t | s=t) $ 
To refine the aprroximation, we can see the above as a convolution of two functions: our unknown $P(X_t | s)$, which we assume smooth around $s=t$, and a quasi delta function, say a gaussian with small variance. Now, we have our first order approximation (for continuous variables) :
$h(x) = g(x)  * N(x_0,\sigma^2)$ (convolution)
$h(x_0) \approx g(x_0)  + g(x_0)''\sigma^2/2$
$g(x_0) \approx h(x_0)  - h''(x_0)''\sigma^2/2$
Applying this to the previous equation can lead to a refined approximation to our desired probability.
A: Just a word of explanation: Part out of curiosity, part for lack of a better more theoretical method, i approached the problem in a completely empirical/inductive way. I'm aware that there is the risk of getting stuck in a dead end without gaining much insight, but i thought, i'll just present what i got so far anyway, in case it is useful to someone.
Starting by computing the exact probabilities for $n,t\in\{1,...,8\}$ we get

Due to the underlying multinomial distribution, multiplying the entries in the table by $(n+1)^t$ leaves us with a purely integer table:

Now we find that there is a polynomial in $n$ for every column which acts as the sequence function for that column:

Dividing the sequence functions by $(n+1)^t$ gives us sequence functions for the original probabilities for the first $t$'s. These rational polynomials can be simplified by decomposing them into partial fractions and substituting $x$ for $1/(n+1)$, leaving us with:

or as coefficient table

Starting with the $x^2$ column there are sequence functions for these coefficients again:

That's how far i got. There are definitely exploitable patterns here that allow sequence functions to occur, but i'm not sure if there is a nice closed form solution for these sequence functions.
