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 be independent exponential variables of parameter , and be an exponential variable of parameter independent of the previous ones. We denote by the ordered sample , i.e. and . Then, for , .

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