# How to solve following equation with multiple random variables

How can I program the following problem into Mathematica?

$$X,Y$$ are i.i.d. Exponentially distributed random variables. Such that $$X\sim\exp(1), \ Y\sim\exp(1)$$.

All other variables are positive real-valued.

$$\Phi =\mathbb{P}\left[P_v \geq A + \frac{B}{Y}\right]$$

## $$P_v= \left\{ \begin{array}{ll} a\left(\frac{b}{1+\exp\left(-\bar \mu\frac{P_s X}{r^\alpha}+\varphi\right)}-1\right), & \text{if}\ \frac{P_s X}{r^\alpha}\geq P_a,\\ 0, & \text{otherwise}. \end{array} \right.$$

My solution

$$\begin{multline} \Phi = \mathbb{P}\left[ a\left(\frac{b}{1+\exp\left(-\bar \mu\frac{P_s X}{r^\alpha}+\varphi\right)}-1\right) \geq A + \frac{B}{Y}\right]\mathbb{P}\left[X \geq \frac{P_ar^\alpha}{P_s}\right] + {\mathbb{P}\left[0 \geq A + \frac{B}{Y}\right] \mathbb{P}\left[X \geq \frac{P_ar^\alpha}{P_s}\right]} \end{multline}$$

Given that $$\mathbb{P}[0>A+\frac{B}{Y}] = 0$$

$$\Phi = \mathbb{P}\left[a\left(\frac{b}{1+\exp\left(-\bar \mu\frac{P_s X}{r^\alpha}+\varphi\right)}-1\right) \geq A + \frac{B}{Y}\right] \exp\left(-\frac{P_ar^\alpha}{P_s}\right)$$

Now consider only the first part of the expression on the right side of the equation, and letting $$c = \frac{\bar{\mu}P_s}{r^{\alpha}}$$.

$$\mathbb{P}\left[a\left(\frac{b}{1+\exp\left(-\bar \mu\frac{P_s X}{r^\alpha}+\varphi\right)}-1\right) \geq A + \frac{B}{Y}\right] = \mathbb{P}\left[\frac{ab}{1+\exp\left(-c X+\varphi\right)}-a \geq A + \frac{B}{Y}\right]$$

consider $$D = A + a$$, with mathematical manipulations and using the fact that $$\mathbb{E}[Y] = \frac{1}{\lambda} = 1$$ and $$\mathbb{E}[e^{-sX}] = \frac{1}{1+s}$$, we get:

$$= \mathbb{P} \left[Y \geq \frac{1}{ab}(DY + B). (1+e^\varphi e^{-c X}) \right] = \mathbb{P} \left[ Y\geq \frac{B (1 + e^\varphi e^{-cX})}{(ab - D)(1+e^\varphi e^{-cX})} \right]$$

$$= \mathbb{P} \left[ Y\geq \mathbb{E}_X\left[\frac{B (1 + e^\varphi e^{-cX})}{ab - D(1+e^\varphi e^{-cX})}\right] \right] \ \ \ (1)$$

My Code in Mathematica is as follows. But the Analysis does not meet Simulation.

Egx = Expectation[(B * (1 + E^\[CurlyPhi]*E^(-c X)))/(a b - D*(1 + E^\[CurlyPhi]*E^(-c X))), X \[Distributed] ExponentialDistribution[1]]

Probability[Y >= Egx, Y \[Distributed] ExponentialDistribution[1]]


My question is if my approach towards the solution correct?

• Have a look at Probability. Dec 23, 2020 at 8:59
• +1 for your suggestion. Please see the additions that I have performed. Can you please suggest first of all if my approach is correct?
– SJa
Dec 23, 2020 at 15:17
• Simultaneously posted at stats.stackexchange.com/questions/501958/….
– JimB
Dec 23, 2020 at 21:29

This is just an extended comment for now.

It seems that because $$P_v$$ and $$A+B/Y$$ are independent, then determining their respective pdf's and/or cdf's should occur first.

The pdf for $$A+B/Y$$ is easily found:

wDist =
TransformedDistribution[A + B/y, y \[Distributed] ExponentialDistribution[1];
PDF[wDist, w]


If maybe you're only interested in $$B>0$$, then there is much simplification. The point is that the signs of some of the parameters becomes important in getting a final answer.

Determining the distribution of $$P_v$$ is not so simple as the distribution is a mixture of a discrete distribution with a probability mass at 0 and the remainder being a continuous distribution with

$$\frac{a (b-1) e^{\mu P_a-\varphi}-a}{e^{\mu P_a-\varphi}+1}

My initial attempts at getting a direct symbolic result for the cdf of $$P_v$$ were not successful. So I plugged in specific values of the parameters for which Mathematica did provide a result. Looking for patterns I ended up with the following for the cdf of $$P_v$$:

Piecewise[
{{1, z >= a*(-1 + b)},
{1 - E^(-((Pa*r^α)/Ps)),
0 <= z <= (-a + a*(-1 + b)*E^(-φ + Pa*μ))/(1 + E^(-φ + Pa*μ))},
{1 - E^((r^α*(-φ + 2*ArcTanh[(-2 + b)/b - (2*z)/(a*b)]))/(Ps*μ)),
(-a + a*(-1 + b)*E^(-φ + Pa*μ))/(1 + E^(-φ + Pa*μ)) < z < a*(-1 + b)}}, 0]


From those distributions one might be able to get a symbolic result for $$\Phi$$.

Update: Still no luck with a completely symbolic solution (and I'm not convinced that a symbolic solution exists but I've been wrong before). Here is a numerical solution that matches with random sampling:

cdf[z_, a_, b_, μ_, Pa_, Ps_, r_, α_, φ_] :=
Piecewise[{{1, z >= a*(-1 + b)},
{1 - E^((r^α*(Pa*μ - φ))/(Ps*μ))*(-1 + (a*b)/(a + z))^(r^α/(Ps*μ)),
a*(-1 + b/(1 + E^(-(Pa*μ) + φ))) < z < a*(-1 + b)}},
0]

pdf[w_, A_, B_] := Piecewise[{{B/(E^(B/(-A + w))*(-A + w)^2), w > A}}, 0]

(* Numerical integration *)
prob[a_, b_, μ_, Pa_, Ps_, r_, α_, φ_, A_, B_] :=
Module[{p},
p = Exp[-((Pa r^α)/Ps)] - Exp[-((Pa r^α)/Ps)] If[A <= a (b/(1 + Exp[φ - Pa μ]) - 1),
NIntegrate[cdf[w, a, b, μ, Pa, Ps, r, α, φ] pdf[w, A, B], {w, a (b/(1 + Exp[φ - Pa μ]) - 1),
a (b - 1)}] + 1 - E^(B/(A - a (-1 + b))),
If[A <= a (b - 1), NIntegrate[cdf[w, a, b, μ, Pa, Ps, r, α, φ] pdf[w, A, B],
{w, A, a (b - 1)}] + 1 - E^(B/(A - a (-1 + b))), 1]];
p]

(* Random sampling *)
sample[a_, b_, μ_, Pa_, Ps_, r_, α_, φ_, A_, B_] := Module[{n, x, Pv, w, Pv2},
n = 1000000;
x = RandomVariate[ExponentialDistribution[1], n];
Pv = If[# >= r^α Pa/Ps, a (b/(1 + Exp[-(μ Ps/r^α) # + φ]) - 1), 0] & /@ x;
w = A + B/RandomVariate[ExponentialDistribution[1], n];
Length[Select[Pv - w, # >= 0 &]]/n // N
]


Now with an example:

prob[4, 6, 3, 1, 1, 2, 1, 3, 1, 3]
(* 0.105301 *)
SeedRandom[12345];
sample[4, 6, 3, 1, 1, 2, 1, 3, 1, 3]
(* 0.104986 *)


• First of all, I HIGHLY appreciate the work you put in for the problem that was not yours. I really do. I am trying the understand your solution and will see how the results match with my simulation. Even if it does not match, I appreciate you took time out.
– SJa
Dec 25, 2020 at 6:18
• One thing, how should I interpret the last case of the final output? I think I have seen this type of case only in Mathematica, but I don't know what does it means that Out is 0 when the case is True.
– SJa
Dec 25, 2020 at 6:19
• I take the "True" to mean a condition that isn't satisfied by any of the above conditions. So in this case the cdf is zero whenever $z\leq \frac{a (b-1) e^{\mu \text{Pa}-\varphi }-a}{e^{\mu \text{Pa}-\varphi }+1}$. (I'll clean up the notation sometime after Christmas. The $z$ should really be $P_v$.)
– JimB
Dec 25, 2020 at 6:38
• Are there any restrictions on the parameters $a$, $b$, $A$, $B$, $\varphi$, $P_a$, $P_s$, $\mu$, $r$, and $\alpha$? I'm guessing that $a \geq 0$, $B>0$, $\bar{\mu}P_s/r^\alpha>0$, and $P_a r^\alpha /P_s>0$.
– JimB
Dec 27, 2020 at 0:10
• Dear @JimB, all parameters are positive real numbers. Specifically, $a>0$, $b\geq 1$, $A,B,\varphi, P_a, P_s, \mu, r, \alpha$ are all greater than 0
– SJa
Dec 27, 2020 at 8:07