# 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 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.

• With large n the independent Poisson approximation is probably fine. Did you try simulation studies of how well the formula is working? – Aniko May 6 '11 at 15:33
• This question is closely related to the Generalized Birthday Problem. en.wikipedia.org/wiki/Birthday_problem – charles.y.zheng May 7 '11 at 20:40
• @Aniko: I haven't run extensive simulations yet. But the examples I tried are roughly correct. – j.p. May 8 '11 at 16:51
• The central difficulty in your problem (and in the birthday problem) is the difficulty of determining the distribution of the maximum (supremum norm) of a multinomial random variable, which involves summing over partitions. – charles.y.zheng May 8 '11 at 20:57
• A hard bound on the Poisson model is as follows. Let $Z_1 \sim \mathrm{Poi}(2t/(n+1))$ and $Z_i \sim \mathrm{Poi}(t/(n+1))$ for $2 \leq i \leq n$. All the $Z_i$ are mutually independent. Then $\mathbb{P}(Z_1 > \max_{i\geq 2} Z_i) \geq 1 - (n-1) \exp(-c t / (n+1))$ where $c = (\sqrt{2}-1)^2$. As you can see, it only works well for $t \geq 6 n \log n$ or so. – cardinal May 15 '11 at 15:22

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.

• Thanks for the answer! Looks very good, but I have to check the details. What do you mean with "I ... corrected the typos by reproducing DasGupta's Table 4, ..."? (By the way, if you had answered 2-3 hours earlier, you'd save me some headaches about what to do with my bounty.) – j.p. May 16 '11 at 14:15
• @pul His inequalities are in the wrong direction: what he claims are $p(n,K,x)$ are really $1 - p(n,K,x-1)$. Sorry about the bounty problem: I knew how to answer this one when it first appeared but needed to check the results first and had no time to do anything about it until the weekend. – whuber May 16 '11 at 15:41

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.

• Could you please tell me how to find $\sigma^2$? – j.p. May 15 '11 at 15:10
• $\sigma^2$ is the variance of $s$, which is the sum of $n$ independent poisson – leonbloy May 15 '11 at 15:39
• @leonbloy: OK, in our case we have therefore $\sigma^2 = t$ (thanks!). And how do I get h"? – j.p. May 15 '11 at 15:56
• I'd approximate the second derivative by the second difference : $A_{t+1}-2A_t+A{t-1}$, the probability of your 'success' event evaluated at different values of $t$ – leonbloy May 15 '11 at 16:30
• @leonbloy: I'm not really convinced of your answer (yet???), but before letting the bounty points vanish into nowhere, I'll accept your answer. – j.p. May 15 '11 at 17:00

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.

• Thanks for the effort! I'm not sure what your results imply for $n$ and $t$ in the hundreds. – j.p. May 15 '11 at 15:14
• I actually tried to approximate the probability for bigger $n,t$ by just taking the $x^2$ part of the series into account, but that only works for small probabilities, for the probabilities that you're interested in, the approximation is way off. – Thies Heidecke May 15 '11 at 15:55
• I'm not convinced of neither your answer (as small $n$ don't seem to help for the $n$'s I need) nor the other (as I don't understand (yet??) its correctness/helpfulness). As I have more hope to get something out of the other answer and as I don't want to let the 100 points vanish, I'll probably accept the other answer. Sorry for not picking yours! – j.p. May 15 '11 at 16:55