Gamma distribution

Error probability or chance to beat

Given two distributions \(X_A \sim \mathcal{G}(\alpha_A, \beta_A)\) and \(X_B \sim \mathcal{G}(\alpha_B, \beta_B)\) such that \((\alpha_A, \beta_A, \alpha_B, \beta_B) \in \mathbb{R}_+^4\), \(P[X_B > X_A]\) is given by

\[P[X_B > X_A] = 1 - \frac{\beta_A^{\alpha_A}\beta_B^{\alpha_B}}{(\beta_A + \beta_B)^{\alpha_A+\alpha_B}}\frac{_2F_1\left(1, \alpha_A + \alpha_B; \alpha_B + 1; \frac{\beta_B}{\beta_B + \beta_A}\right)}{\alpha_B B(\alpha_A, \alpha_B)} = I_{\frac{\beta_A}{\beta_A + \beta_B}}(\alpha_A, \alpha_B),\]

where \(_2F_1(a,b;c;z)\) is the Gauss hypergeometric function and \(I_x(a,b)\) is the regularized incomplete beta function.

Expected loss function

The expected loss function can easily be calculated from the definition yielding

\[\mathrm{EL}(X_B) = \frac{\alpha_A}{\beta_A} I_{\frac{\beta_B}{\beta_A + \beta_A}}(\alpha_B, \alpha_A + 1) - \frac{\alpha_B}{\beta_B} I_{\frac{\beta_B}{\beta_A + \beta_A}}(\alpha_B + 1, \alpha_A).\]

A similar expression is obtained for \(\mathrm{EL}(X_A)\),

\[\mathrm{EL}(X_A) = \frac{\alpha_B}{\beta_B} I_{\frac{\beta_A}{\beta_A + \beta_A}}(\alpha_A, \alpha_B + 1) - \frac{\alpha_A}{\beta_A} I_{\frac{\beta_A}{\beta_A + \beta_A}}(\alpha_A + 1, \alpha_B)\]

Credible intervals

Credible intervals are employed to account for uncertainty in the expected loss and relative expected loss measures. Let us considered the relative expected loss if variant B is chosen, which follows the distribution \((X_A - X_B)/X_B = X_A / X_B - 1\). This requires the distribution of the ratio of two random gamma variables, \(U = X_A / X_B\). The probability density function is given by

\[f(u) = \left(\frac{\beta_B}{\beta_A}\right)^{\alpha_B} \frac{u^{\alpha_B - 1} (1 + \frac{\beta_B}{\beta_A}u)^{-\alpha_A - \alpha_B}}{B(\alpha_A, \alpha_B)}.\]

Note that this is the probability density function of the generalized beta prime distribution. The cumulative distribution function is given by

\[F(u) = I_{\frac{u}{u + \frac{\beta_B}{\beta_A}}}(\alpha_A, \alpha_B).\]

Note

Credible intervals are computed by solving \(F(u) = p\), \(p \in [0, 1]\). A reasonable starting point is the normal approximation of the gamma distribution.

The expected value and variance of the distribution \(Z = (X_A - X_B)/X_B = X_A / X_B - 1\) can be computed using

\[\mathrm{E}\left[\frac{X_A}{X_B} \right] = \frac{\alpha_A}{(\alpha_B - 1)}\frac{\beta_B}{\beta_A}.\]
\[\mathrm{Var} \left[\frac{X_A}{X_B} \right] = \frac{\alpha_A (\alpha_A + \alpha_B - 1)}{(\alpha_B - 2)(\alpha_B - 1)^2} \left(\frac{\beta_B}{\beta_A}\right)^2.\]

Proofs

Error probability

Integrating the joint distribution over all values of \(X_B > X_A\) we obtain the integral

\[\begin{split}P[X_B > X_A] &= \int_0^{\infty} \int_{x_A}^{\infty} \frac{\beta_A^{\alpha_A}}{\Gamma(\alpha_A)} x_A^{\alpha_A - 1} e^{-\beta_A x_A} \frac{\beta_B^{\alpha_B}}{\Gamma(\alpha_B)} x_B^{\alpha_B - 1} e^{-\beta_B x_A} \mathop{dx_B}\mathop{dx_A}\\ &= 1 - \int_0^{\infty}\frac{\beta_A^{\alpha_A}}{\Gamma(\alpha_A)} x_A^{\alpha_A - 1} e^{-\beta_A x_A} P(\alpha_B, \beta_B x_A)\mathop{dx_A},\end{split}\]

where \(P(a,z)\) is the regularized lower incomplete gamma function defined by

\[P(a, z) = \frac{\gamma(a, z)}{\Gamma(a)} = 1 - Q(a,z),\]

and \(Q(a,z)\) is the regularized upper incomplete gamma function with series expansion

\[Q(a,z) = \frac{\Gamma(a,z)}{\Gamma(a)} = 1 - z^a e^{-z} \sum_{k=0}^{\infty}\frac{z^k}{\Gamma(a+k+1)}.\]

Hence, the integral is rewritten in the form

\[P[X_B > X_A] = \int_0^{\infty}\frac{\beta_A^{\alpha_A}}{\Gamma(\alpha_A)} x^{\alpha_A - 1} e^{-\beta_A x} Q(\alpha_B, \beta_B x) \mathop{dx}.\]

Interchange of integration and summation leads to a representation in terms of the Gauss hypergeometric function \(_2F_1(a,b;c;z)\), also expressible in terms of the incomplete beta function

\[\begin{split}P[X_B > X_A] &= 1 - \frac{\beta_A^{\alpha_A}\beta_B^{\alpha_B}}{\Gamma(\alpha_A)}\int_0^{\infty} x^{\alpha_A + \alpha_B - 1} e^{-(\beta_A + \beta_B) x} \sum_{k=0}^{\infty}\frac{(\beta_B x)^k}{\Gamma(\alpha_B + k + 1)} \mathop{dx}\\ &= 1 - \frac{\beta_A^{\alpha_A}\beta_B^{\alpha_B}}{\Gamma(\alpha_A)} \sum_{k=0}^{\infty} \frac{\beta_B^k}{\Gamma(\alpha_B + k + 1)}\int_0^{\infty}x^{\alpha_A + \alpha_B + k - 1} e^{-(\beta_A + \beta_B) x}\mathop{dx}\\ & =1 - \frac{\beta_A^{\alpha_A}\beta_B^{\alpha_B}}{(\beta_A + \beta_B)^{\alpha_A+\alpha_B}}\sum_{k=0}^{\infty}\frac{\Gamma(\alpha_A + \alpha_B + k)}{\Gamma(\alpha_B + k + 1) \Gamma(\alpha_A)} \left(\frac{\beta_B}{\beta_A + \beta_B}\right)^k\\ &=1 - \frac{\beta_A^{\alpha_A}\beta_B^{\alpha_B}}{(\beta_A + \beta_B)^{\alpha_A+\alpha_B}}\frac{_2F_1\left(1, \alpha_A + \alpha_B; \alpha_B + 1; \frac{\beta_B}{\beta_B + \beta_A}\right)}{\alpha_B B(\alpha_A, \alpha_B)}= I_{\frac{\beta_A}{\beta_A + \beta_B}}(\alpha_A, \alpha_B).\end{split}\]