# 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! – xzczd Nov 17 '14 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? – hengxin Nov 17 '14 at 12:37
• Network is a bit slow today. See my answer. – xzczd Nov 17 '14 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.) – hengxin Nov 18 '14 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. – xzczd Nov 20 '14 at 4:24