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.

  • 1
    $\begingroup$ With large n the independent Poisson approximation is probably fine. Did you try simulation studies of how well the formula is working? $\endgroup$ – Aniko May 6 '11 at 15:33
  • $\begingroup$ This question is closely related to the Generalized Birthday Problem. en.wikipedia.org/wiki/Birthday_problem $\endgroup$ – charles.y.zheng May 7 '11 at 20:40
  • $\begingroup$ @Aniko: I haven't run extensive simulations yet. But the examples I tried are roughly correct. $\endgroup$ – j.p. May 8 '11 at 16:51
  • 1
    $\begingroup$ 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. $\endgroup$ – charles.y.zheng May 8 '11 at 20:57
  • 1
    $\begingroup$ 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. $\endgroup$ – 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.

  • $\begingroup$ 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.) $\endgroup$ – j.p. May 16 '11 at 14:15
  • $\begingroup$ @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. $\endgroup$ – 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.

  • $\begingroup$ Could you please tell me how to find $\sigma^2$? $\endgroup$ – j.p. May 15 '11 at 15:10
  • $\begingroup$ $\sigma^2$ is the variance of $s$, which is the sum of $n$ independent poisson $\endgroup$ – leonbloy May 15 '11 at 15:39
  • $\begingroup$ @leonbloy: OK, in our case we have therefore $\sigma^2 = t$ (thanks!). And how do I get h"? $\endgroup$ – j.p. May 15 '11 at 15:56
  • $\begingroup$ 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$ $\endgroup$ – leonbloy May 15 '11 at 16:30
  • $\begingroup$ @leonbloy: I'm not really convinced of your answer (yet???), but before letting the bounty points vanish into nowhere, I'll accept your answer. $\endgroup$ – 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

Table of the first few probabilities for low n,t

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

Table of integerified probabilities

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

Sequence functions for different t's

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:

Sequence functions in x=1/(n+1)

or as coefficient table

sequence function coefficients

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

x^k coefficients sequence functions

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.

  • $\begingroup$ Thanks for the effort! I'm not sure what your results imply for $n$ and $t$ in the hundreds. $\endgroup$ – j.p. May 15 '11 at 15:14
  • $\begingroup$ 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. $\endgroup$ – Thies Heidecke May 15 '11 at 15:55
  • $\begingroup$ 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! $\endgroup$ – j.p. May 15 '11 at 16:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.