After reading your last post I tried checking my formula with simulations myself (Scilab code below), and it happens to match pretty well... So maybe we're not dealing with exactly the same problem. Here is the problem that I answered:

Let $\displaystyle Y_1,\ldots,Y_M$ be independent exponential variables of parameter $\displaystyle \lambda$, and $\displaystyle Z$ be an exponential variable of parameter $\displaystyle \lambda N$ independent of the previous ones. We denote by $\displaystyle (I_1,\ldots,I_M)$ the ordered sample $\displaystyle (Y_1,\ldots,Y_M)$, i.e. $\displaystyle I_1<I_2<\cdots<I_M$ and $\displaystyle \{I_1,\ldots,I_M\}=\{Y_1,\ldots,Y_M\}$. Then, for $\displaystyle k=1,\ldots,M$, $\displaystyle P(Z>I_k)=\frac{M(M-1)\cdots(M-k+1)}{(M+N)(M+N-1)\cdots(M+N-k+1)}$.

Code:

M=20;
N=8;
lambda=10;
// Simulation
S=zeros(M,1); // S(k) will contain the empirical P(Z>I_k)
nb=50000; // nb of simulations (probability is an average over nb simulations)
for i=1:nb,
Z=-log(rand())/(N*lambda);
Y=-log(rand(M,1))/lambda;
Y=gsort(Y,'g','i'); // ordering of the sample
I=find(Z>Y); // vector I contains the indices i such that Z>Y(i)
for j=I,
S(j)=S(j)+1;
end;
end;
S=S/nb;
clf;plot(1:M,S,'bx'); // Simulated plot in blue
// Theoretical computation
T=zeros(M,1); // T(k) will contain theoretical P(Z>I_k)
u=M; v=N+M;
T(1)=u/v;
for i=2:M,
u=u-1;
v=v-1;
T(i)=T(i-1)*u/v;
end;
plot(1:M,T,'rx'); // Theoretical plot in red