# How to get a more compact form of this probability calculation?

Inspired by the probability calculation here, I am trying to solve a little general one:

$$\mathbb{P} (\sum_{i=1}^{m-1} A_i + \sum_{i=1}^{m} S_i < L < \sum_{i=1}^{m} A_i + \sum_{i=1}^{m+1} S_i)$$

where,

$$A_i \sim \textrm{exp}(\lambda), \; S_i \sim \textrm{exp}(\mu), \; L \sim \textrm{exp}(\lambda), \textrm{ and } \lambda \neq \mu \textrm{ are two integers}.$$

All $A_i, S_i, \textrm{ and } L$ are mutually independent. $m$ is an integer parameter.

In this calculation, I fix $m = 4$. Because $A_i \sim \textrm{exp} (\lambda)$, we have $\sum_{i=1}^{i=3} A_i \sim \textrm{Erlang}(3, \lambda)$.

Probability[
a3 + s4 < middle < a4 + s5,
{a3 \[Distributed] ErlangDistribution[3, λ],
s4 \[Distributed] ErlangDistribution[4, μ],
a4 \[Distributed] ErlangDistribution[4, λ],
s5 \[Distributed] ErlangDistribution[5, μ],
middle \[Distributed] ExponentialDistribution[λ]}]


After a few minutes, I got a long result from which I could not read anything useful:

$\frac{\mu ^4 \left(186624 \lambda ^{16}+3825792 \lambda ^{15} \mu +35971776 \lambda ^{14} \mu ^2+205589664 \lambda ^{13} \mu ^3+798109200 \lambda ^{12} \mu ^4+2217713192 \lambda ^{11} \mu ^5+4537481548 \lambda ^{10} \mu ^6+6954729890 \lambda ^9 \mu ^7+8071898695 \lambda ^8 \mu ^8+7133161040 \lambda ^7 \mu ^9+4799247376 \lambda ^6 \mu ^{10}+2441453824 \lambda ^5 \mu ^{11}+923808416 \lambda ^4 \mu ^{12}+252053248 \lambda ^3 \mu ^{13}+46855424 \lambda ^2 \mu ^{14}+5307904 \lambda \mu ^{15}+276224 \mu ^{16}\right)}{11664 (\lambda +\mu )^5 (2 \lambda +\mu )^7 (\lambda +2 \mu )^8}$

In addition, the command FullSimplify[Probability[]] does not help.

Edit (2014-11-17):

By paper-and-pencil, I have obtained a closed form $\frac{\mu + 2 \lambda}{\mu + \lambda} (\frac{1}{2} \frac{\mu}{\mu + \lambda})^{m}$. Notice that the sum $\sum_{m=1}^{\infty} \frac{\mu + 2 \lambda}{\mu + \lambda} (\frac{1}{2} \frac{\mu}{\mu + \lambda})^{m} = \frac{\mu}{\mu + \lambda}.$ The remaining $\frac{\lambda}{\mu + \lambda}$ for special case $m = 0$ has been omitted here.

Therefore, the question is

How to get a more compact form of the above probability? Specifically, is it consistent with my manual calculation (which may be wrong)?

• I think your program is wrong. In your program, a3 and a4, s4 and s5 are independent, but actually they're not, a3 is part of a4, s4 is part of s5! Nov 17, 2014 at 12:23
• @xzczd Thanks. Even "error" is helpful to me. How should I program the non-independent? Is it right to first define each $A_i$ and $S_i$ separately and then sum them to get the Erlang distributions? Nov 17, 2014 at 12:37
• Network is a bit slow today. See my answer. Nov 17, 2014 at 13:39

The approach that kguler suggest in your precedent question is completely suitable for the current one:

FindSequenceFunction[Table[With[{
d1 = HypoexponentialDistribution[Flatten[Table[{λ, μ}, {i - 1}]]],
d2 = HypoexponentialDistribution[{λ, μ}],
d3 = ExponentialDistribution[μ],
d = ExponentialDistribution[λ]},
Probability[l1 + l3 < l < l1 + l2 + l3,
{l1 \[Distributed] d1, l2 \[Distributed] d2,
l3 \[Distributed] d3, l \[Distributed] d}]], {i, 2, 5}]][m - 1] // Simplify //
AbsoluteTiming


$\left\{108.573200,\frac{2^{-m} (2 \lambda +\mu ) \left(\frac{\mu }{\lambda +\mu }\right)^m}{\lambda +\mu }\right\}$

Of course ErlangDistribution can also be used, and it's faster than the above approach actually:

FindSequenceFunction[
Table[Probability[
a + s < middle < a + aadd + s + sadd,
{a \[Distributed] ErlangDistribution[i, λ],
s \[Distributed] ErlangDistribution[i + 1, μ],
middle \[Distributed] ExponentialDistribution[λ]}],
{i, 4}]][m - 1] // Simplify // AbsoluteTiming


$\left\{37.932000,\frac{2^{-m} (2 \lambda +\mu ) \left(\frac{\mu }{\lambda +\mu }\right)^m}{\lambda +\mu }\right\}$

You can use ParallelTable instead of Table to speed up even more.

BTW, there exists a more straight-forward but much slower approach:

With[{i = 2}, Probability[
Total@Array[a, i - 1] + Total@Array[s, i] <l< Total@Array[a, i] + Total@Array[s, i + 1],
Join[{l \[Distributed] ExponentialDistribution[λ]},
a[#] \[Distributed] ExponentialDistribution[λ]&/@Range[i],
s[#] \[Distributed] ExponentialDistribution[μ]&/@Range[i + 1]]]] // AbsoluteTiming


$\left\{96.188400,\frac{\mu ^2 (2 \lambda +\mu )}{4 (\lambda +\mu )^3}\right\}$

• A long way for me to go to learn Mathematica. Thanks a lot. (I have an even worse network these days. I should have accepted your answer yesterday.) Nov 18, 2014 at 2:13
• @hengxin It's OK to wait 24 hours or even longer before accepting so your question may attract better answers :) . I added a ErlangDistribution approach, have a look. Nov 20, 2014 at 4:24