mertcan said:
@Ray Vickson thanks for your return, I thought it was wrong due to fact that I saw different answer on a internet site while I was digging something related to my problem, there was no derivation of it on that site but just wrote down the answer differs from actual. But you know better of course that site actually not as good as here :) By the way what is your suggestion for the case of $$P(S^n_A<S^m_B<S^k_C)$$ What do we have to do here ??
I will write ##a,b,c## instead of ##\lambda_A, \lambda_B, \lambda_C##. Also, I will write ##A_n, B_m, C_k## instead of ##S^n_A, S^m_B, S^k_C##. Note that ##A_n## has distribution ##E_n(a)## = ##n##-Erlang with rate parameter ##a##. Its density is
$$f_{A_n}(t) = f_n(t,a) \equiv \frac{a^n t^{n-1}}{(n-1)!} e^{-at} \; 1_{\{ t > 0 \}} $$.
Similarly, ##B_m## has ##m##-Erlang density ##f_m(t,b)## and ##C_k## has density ##f_k(t,c)##. Thus
$$P(A_n < B_m < C_k) = \int_{t_1=0}^{\infty} d t_1 \,f_n(t_1,a) \int_{t_2=t_1}^{\infty} dt_2 \, f_m(t_2,b) \int_{t_3=t_2}^{\infty} dt_3 \, f_k(t_3,c)$$.
This can be evaluated by repeated application of the following result. If ##p_i(\mu) = \mu^i e^{-\mu}/i!## is a Poisson probability, then
$$\int_t^{\infty} f_L(x,r) \, dx = \sum_{i=0}^{L-1} p_i(rt)$$.
You can get this in two ways:
Method (1): Direct--using repeated integration-by-parts, starting with ##u = x^{L-1}, dv = e^{-rx} \, dx##, etc., so that the integral with ##x^{L-1}## in it is expressed in terms of an integral with ##x^{L-2}##, then integrate-by-parts again to get ##x^{L-3}##, etc. Eventually you will arrive at the previous formula.
Method (2): "Probabilistic" (more clever and a lot easier): The event that the ##L##th arrival is later than ##t## is the same as the event that the number of arrivals before time ##t## is ## \leq L-1##. That is, ##P(S_L > t) = P(N[0,t] \leq L-1)##, and for a Poisson process with rate ##r## we have ##P(N[0,t] \leq L-1) = \sum_{i=0}^{L-1} p_i(rt)##. That gives the same formula.
So
$$f_m(t_2,b) \int_{t_3=t_2}^{\infty} f_k(t_3,c) dt_3 = \sum_{i=0}^{k-1} f_m(t_2,b) (c t_2)^i e^{-c t_2}/ i! \\= \sum_{i=0}^{k-1} \frac{b^m c^i}{(m-1)! i!} t_2^{m+i-1} e^{-(b+c) t_2}$$
By the way, if you were to integrate this with respect to ##t_2## from ##t_2 = 0## to ##t_2 = +\infty## you would have the expression for ##P(B_m < C_k)##, and if you finish the algebra you will get the formula you previously cited for that probability.
Anyway, we now have
$$\int_{t_2=t_1}^{\infty}dt_2\, f_m(t_2,b) \int_{t_3=t_2}^{\infty} f_k(t_3,c) dt_3 = \sum_{j=0}^{k-1} \frac{b^m c^i}{i! (m-1)!}
\int_{t_1}^{\infty} t_2^{m-1+i} e^{-(b+c) t_2} \, dt_2$$
The integrand in the ##i##th term is ##[(m-1+i)!/(b+c)^{m+i}] f_{m+i}(t_2,b+c)##, so the ##t_2##-integral can be expressed as
$$\int_{t_1}^{\infty} t_2^{m-1+i} e^{-(b+c) t_2} \, dt_2 = \frac{(m-1+i)!}{(b+c)^{m+i} }\sum_{j=0}^{m+i-1} p_j((b+c)t_1).$$
So, finally when we multiply this by ##f_n(t_1,a)## and integrate from ##t_1=0## to ##t_1 = +\infty## we will have an answer that is a somewhat complicated double sum of the form ##\sum_{j=0}^{m+i-1} \sum_{i=0}^{k-1} ( \cdots ) ##; you can fill in the details. I don't know offhand if the double sum can be simplified (using some other known formulas) but you can work on that if you wish.
Note: when we integrate from ##t_1=0## to ##+\infty## we do not obtain a third sum, because when we integrate an Erlang density from 0 to ##\infty## we obtain 1; we would obtain another sum only when the lower limit of integration is ##> 0##.