Per this answer, if $X \sim Gamma(\alpha_1, \beta_1)$ and $Y \sim Gamma(\alpha_2, \beta_2)$, then $$P[X > Y] = H_{\alpha_2, \alpha_1} \left(\frac{\beta1}{\beta1+\beta2}\right)$$ where $H$ is the CDF of a beta distribution. My question is, how would I calculate $$P[X + k > Y]$$ where $k$ is a real number? I don't see a way to adapt this proof of the above. I recognize that this is equivalent to asking for the CDF of the difference distribution $Y - X$ evaluated at $k$, but I haven't had any luck finding a closed form solution for that. I would also be interested in a good approximation if no closed form solution exists.
$\begingroup$
$\endgroup$
5
-
2$\begingroup$ Does this answer your question? probability of gamma greater than exponential $\endgroup$– Xi'anCommented Nov 6, 2020 at 7:42
-
$\begingroup$ Those do not. They show how to calculate $P[X > Y]$ but don't generalize to $P[X > Y + k]$. $\endgroup$– JacobCommented Nov 6, 2020 at 17:09
-
$\begingroup$ Are we looking at different things? All three answers pretty explicitly give $P[X > Y]$. $\endgroup$– JacobCommented Nov 7, 2020 at 19:45
-
$\begingroup$ This answer contains the following$$P \left[ W = \frac{\beta_1Y}{\beta_1Y+\beta_2X}<\frac{\beta_1}{\beta_1+\beta_2} \right]=H_{\alpha_2,\alpha_1} \left( \frac{\beta_1}{\beta_1+\beta_2} \right),$$which obviously generalises to any $k$. $\endgroup$– Xi'anCommented Nov 8, 2020 at 9:38
-
1$\begingroup$ You're leaning hard on "obviously." I recommend attempting that derivation yourself. The algebra in the final step does not generalize to $W = \frac{\beta_1Y + \beta_1k}{\beta_1Y + \beta_2X + \beta_1k}$ or $W = \frac{\beta_1Y}{\beta_1Y + \beta_2X + \beta_2k}$ for $k \neq 0$. In fact, $W$ is not beta-distributed in that case. $\endgroup$– JacobCommented Nov 9, 2020 at 15:50
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
I found a decent approach, with the help of this paper, which as a lot of good info on the subject. You can use: $$P(X - Y \leq k) = \int_0^{\infty} CDF_X(x+k)*PDF_Y(x)dx$$
I don't think it has a closed form solution, but for my purposes numeric integration is enough, which I did with scipy.integrate.quad
in python:
from scipy.stats import gamma
from scipy.integrate import quad
import numpy as np
alpha1 = 24.99495
beta1 = 4.921394
alpha2 = 24.04853
beta2 = 4.435199
k = 10.0
def integrand(x):
return gamma.cdf(x+k, alpha1, 0, beta1)\
*gamma.pdf(x, alpha2, 0, beta2)
print(quad(integrand, 0, np.inf)[0]) ##P(X - Y <= k)
>>> 0.5746310306610378