Integral of a power and exponential of a sum of powers: generalizing Gradshteyn and Ryzhik's 3.478(4)

Question

Gradshteyn and Ryzhik in 3.478(4) give

$$\int_0^\infty x^{\nu-1} \exp\left( -\beta x^p -\gamma x^{-p} \right) dx = \frac{2}{p} \left(\frac{\gamma}{\beta}\right)^{\nu/2p} K_{\nu/p}\left( 2 \sqrt{\beta\gamma} \right) $$

for $\mathrm{Re}(\beta) > 0$ and $\mathrm{Re}(\gamma) > 0$.

My question is, is there a solution to the integral where we replace the $-p$ power with $-1$? I.e.,

$$\int_0^\infty x^{\nu-1} \exp\left( -\beta x^p -\gamma x^{-1}\right) dx = ?$$

(Or more generally any power, so replacing $\gamma x^{-p}$ in the original with $\gamma x^{-q}$?)

---

Some background—this integral arises when applying the distribution of the product of two random variables to

• a Gamma random variable multiplied by
• a Generalized Gamma random variable (specifically, a Gamma random variable to a power)

as described in this thread on the Stan forum. Specifically, when the second Gamma is raised to some integer power $n$, $p=1/n$ in my question while $q=1$.

Given that we can find closed-form expressions of the densities when $p=q=1$, i.e., when we have the product of two Gamma random variables, via G&R as well as right here on Math Stack Exchange, I'm somewhat hopeful we can find an analytical density for when $p\neq q$.

Answer 1

(Edit: I replaced the old discussion of specific integrals in this answer with a generalization. See the history for the original.)

Wikipedia references Sagias and Karagiannidis which give the solution to a slightly different integral:

$$ \int_0^\infty \frac{k}{l} x^n (x / l)^{k - 1} \exp(-t x - (x / l)^k) \, \mathrm{d}x = \frac{(p / t)^{k + n} \sqrt{k}}{l^k \left(\sqrt {2 \pi}\right)^{p + q - 2} } G^{q,p}_{p,q}\left( \left. \begin{matrix} \frac{1-k-n}{p}, \frac{2-k-n}{p}, \dots, \frac{p-k-n}{p} \\ \frac{0}{q}, \frac{1}{q}, \dots, \frac{q-1}{q} \end{matrix} \right| \frac{1}{q^q} \left( \frac{p}{l t}\right)^p \right) $$

This is actually the expectation $E[X^n e^{-tX}]$ where $X$ is a Weibull-distributed random variable with parameters $k>0$ and $l>0$, and crucially, for rational $k=p/q$, that is, $p$ and $q$ are positive integers. (I truly apologize for mixing up the notation between the question and this answer, which uses Wikipedia's notation.) $G$ is the Meijer G function and $n\geq0$ and $t>0$.

Note that, whereas the question asks for $e^{-t/x}$ inside the integral, the above expression is $e^{-tx}$. But we can show with Sympy that the integral of the desired expression is similar:

$$ \int_0^\infty \frac{k}{l} x^n (x / l)^{k - 1} \exp(-t / x - (x / l)^k) \, \mathrm{d}x = \frac{(t/p)^{k + n} \sqrt{k}}{l^k \left(\sqrt {2 \pi}\right)^{p + q - 2} } G^{p+q,0}_{0,p+q}\left( \left. \begin{matrix} - \\ \frac{1-k-n}{p}, \frac{2-k-n}{p}, \dots, \frac{p-1-k-n}{p}, \left(-\frac{1}{q}-\frac{n}{p}\right), \frac{0}{q}, \frac{1}{q}, \dots, \frac{q-1}{q} \end{matrix} \right| \frac{1}{q^q} \left( \frac{t}{l p}\right)^p \right) $$ Note here that the integrand here is almost exactly the same as the first, but with $e^{-t/X}$, so this result is the expectation $E[X^n e^{-t/X}]$ for Weibull $X$.

With such a small change in the integrand, it's not surprising that the results have much in common. The second reciprocates some leading constants and a part of the argument to the Meijer G function, and the array parameterizing the second Meijer G function is largely the union of the parameters of the first (with a strange exception: instead of $(p-k-n)/p$ which is in the first, the second has $(p-k-n)/p-1 = -1/q-n/p$).

This expression matches numerical integration for all values I've tested (mainly for low values of $q$, for example $k=1/5$ to $4/5$).

I'd love to be able to simplify this further, but for now this is acceptable, mpmath can evaluate the Meijer G function via a weighted combination of hypergeometric functions.

N.B. I don't quite understand many of the discussions around the convergence criteria for the integral represented by the Meijer G function, but I will note that many of the common expressions to simplify this to sums of hypergeometric functions fail because for example, with $n=0$, the difference between two of the Meijer G parameters, $(q-1)/q$ and $-1/q$ is an integer (1). I expect that the reason that some simplifications of this integral might fail is because they violate the convergence requirements, e.g., Przemo's solution?

Here's the Python code that I used in the above answer and that compares the Meijer G solution with numerical integration:

import mpmath as mp
import sympy as s
from sympy import S
def makeIntegrand2(n, k, l, t):
  return lambda x: xn * (x / l)(k - 1) * s.exp(-t * x - (x / l)*k)  k / l
def makeIntegrand3(n, k, l, t):
  return lambda x: xn * (x / l)(k - 1) * s.exp(-t / x - (x / l)*k)  k / l
k, l, x, t, n, v = s.symbols('k l x t n nu', real=True, positive=True)
f = s.exp(-t / x) * x(k - 1) * s.exp(-l * (x)k)
extraSub = {n: 2, l: 1, t: 10}
p, q = S(2), S(4)
for pp in range(1, q):
  p = S(pp)
  print(f'# p/q={p}/{q}')
  if False:
    # This is the Wikipedia/Sagias integral
    print('## E[X^n exp(-t X)]')
    f2 = s.simplify(xn * (x / l)(k - 1) * s.exp(-t * x - (x / l)k) * k / l)
    res2 = s.integrate(f2.subs({k: p / q}), (x, 0, s.oo)).simplify()
    s.pprint(res2)
    meijerList2 = ([(x - (p / q) - n) / p for x in range(1, p + 1)], [x / q for x in range(int(q))])
    prefix = (p / t)(p / q + n) * s.sqrt(p / q) / l(p / q) / s.sqrt(2 * s.pi)(p + q - 2)
    arg = (p / l / t)p / qq
    print(meijerList2) # parameters of the G function
    s.pprint(prefix) # this is before the G function
    s.pprint(arg) # this goes inside the G function
    # continue
print('## E[X^n exp(-t / X)]')
  f3 = s.simplify(xn * (x / l)(k - 1) * s.exp(-t / x - (x / l)*k)  k / l)
  res3 = s.integrate(f3.subs({k: p / q}), (x, 0, s.oo)).simplify()
  s.pprint(res3)
meijerList = [(x - (p / q) - n) / p for x in range(1, p)] + [-1 / q - n / p
                                                              ] + [x / q for x in range(int(q))]
  recon = (t / p)(p / q + n) * s.sqrt(p / q) / l(p / q) / s.sqrt(
      2 * s.pi)(p + q - 2) * s.meijerg(((), ()), (meijerList, ()), (t / l / p)p / q**q)
  s.pprint(recon)
qres3 = mp.quad(
      makeIntegrand3(n=extraSub[n], k=float(pp / q), l=extraSub[l], t=extraSub[t]), [0, mp.inf])
  print('comparison', [qres3, res3.subs(extraSub).evalf(), recon.subs(extraSub).evalf()])
```

Answer 2

Let us take $x >0 $ and $y>0$ and $\alpha >0$, $\gamma_1>0$ and $\gamma_2>0$. Then we define:

\begin{equation} {\mathfrak J}_{\alpha,\gamma_1,\gamma_2}(x,y):= \int\limits_0^\infty t^{\alpha-1} \exp(-y t^{\gamma_1}) \exp(-\frac{x}{t^{\gamma_2}}) dt \end{equation}

Then we have:

\begin{eqnarray} &&{\mathfrak J}_{\alpha,\gamma_1,\gamma_2}(x,y) = \tag{1} \\ &&\frac{y^{-\frac{\alpha}{\gamma_1}}}{\gamma_1} \cdot \frac{1}{2\pi \imath} \int\limits_{-\imath \infty+ \zeta}^{\imath \infty +\zeta} \Gamma\left( \frac{\alpha+\gamma_2 s}{ \gamma_1} \right) \cdot \Gamma\left(s\right) \cdot \left[ x y^{\frac{\gamma_2}{\gamma_1}} \right]^{-s} ds = \tag{2} \\ && \frac{y^{-\frac{\alpha}{\gamma_1}}}{\gamma_1} \cdot \sum\limits_{n=0}^{\infty} \frac{(-1)^n}{n!} \left[ \frac{\gamma_1}{\gamma_2} \cdot \Gamma[-\frac{n \gamma_1+\alpha}{\gamma_2}] \cdot (x y^{\frac{\gamma_2}{\gamma_1}})^{\frac{n \gamma_1+\alpha}{\gamma_2}} + \Gamma[\frac{\alpha-n \gamma_2}{\gamma_1}] \cdot (x y^{\frac{\gamma_2}{\gamma_1}})^n \right] \tag{3} \end{eqnarray}

In $(2)$ we expressed the integrand through the inverse Mellin transform, which boils down to using the approach from here the first answer to this question. Here $\zeta > 0 $.Then in $(3)$ we used Cauchy theorem to compute the integral in question. The integral is then equal to the sum of the residues of the integrand at poles of order one being located at firstly $(\alpha + \gamma_2 s_n)/\gamma_1 = - n $ and secondly at $s_n = - n $ where $n \in {\mathbb N} $.

In[632]:= {alpha, g1, g2} = 
 RandomReal[{1, 2}, 3]; M = 20;(*g1=1;g2=1;*)
{x, y} = RandomReal[{1, 2}, 2];
NIntegrate[t^(alpha - 1) Exp[-y t^g1] Exp[-x/t^g2], {t, 0, Infinity}]
ofst = -alpha/2;
y^(-alpha/g1) /
  g1 NIntegrate[ 
   Gamma[(alpha + g2 (I s - ofst))/g1] Gamma[
     I s - ofst] (x y^(g2/g1))^(-I s + ofst), {s, -Infinity, 
    Infinity}]/(2 Pi)
Take[Accumulate[
  y^(-alpha/g1) /
    g1 Table[(-1)^n/n! ( 
      g1/g2 Gamma[-(n g1 + alpha)/
          g2] (x y^(g2/g1))^((n g1 + alpha)/g2) + 
       Gamma[(alpha - n g2)/g1] (x y^(g2/g1))^n), {n, 0, M}]], -5]
Out[633]= 0.137716
Out[635]= 0.137716 + 0. I
Out[636]= {0.137716, 0.137716, 0.137716, 0.137716, 0.137716}