1

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.

Jacob
  • 127
  • 2
  • 2
    Does this answer your question? [probability of gamma greater than exponential](https://stats.stackexchange.com/questions/264861/probability-of-gamma-greater-than-exponential) – Xi'an Nov 06 '20 at 07:42
  • Those do not. They show how to calculate $P[X > Y]$ but don't generalize to $P[X > Y + k]$. – Jacob Nov 06 '20 at 17:09
  • Are we looking at different things? All three answers pretty explicitly give $P[X > Y]$. – Jacob Nov 07 '20 at 19:45
  • This [answer](https://stats.stackexchange.com/a/264882/7224) contains the following$$P \left[ W = \frac{\beta_1Y}{\beta_1Y+\beta_2X} – Xi'an Nov 08 '20 at 09:38
  • 1
    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. – Jacob Nov 09 '20 at 15:50

1 Answers1

0

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
Jacob
  • 127
  • 2