0
$\begingroup$

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]]

Mathematica solution

My question is if my approach towards the solution correct?

Simulation and Analysis

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

1 Answer 1

3
$\begingroup$

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]

pdf of A+B/y

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}<P_v<a (b-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]

CDF for Pv

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 *)

I'll add more explanation later.

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – SJa
    Commented Dec 25, 2020 at 6:18
  • $\begingroup$ 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. $\endgroup$
    – SJa
    Commented Dec 25, 2020 at 6:19
  • $\begingroup$ 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$.) $\endgroup$
    – JimB
    Commented Dec 25, 2020 at 6:38
  • $\begingroup$ 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$. $\endgroup$
    – JimB
    Commented Dec 27, 2020 at 0:10
  • $\begingroup$ 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 $\endgroup$
    – SJa
    Commented Dec 27, 2020 at 8:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.