Solving an integral equation with inverse Laplace transform

Question

Let $\alpha,\beta,\mu>0$. I am looking for a solution, i.e. a function $g(x)$, that satisfies $$ \frac{\beta^{\alpha}}{\Gamma(\alpha)}\int_0^\infty g(x)x^{\alpha-1}e^{-\beta x}\,\mathrm dx=\left(\frac{\alpha}{\beta}-\mu\right)^{-1}, $$ where $\alpha/\beta>\mu$. Note that such a solution would yield an unbiased estimator for $\left(\frac{\alpha}{\beta}-\mu\right)^{-1}$, i.e. if $X\sim\operatorname{Gamma}(\alpha,\beta)$ then $\operatorname Eg(X)=\left(\frac{\alpha}{\beta}-\mu\right)^{-1}$. I tried solving this with an inverse Laplace transform by writing $$ \mathcal L\left\{x^{\alpha-1}g(x)\right\}(\beta)=\frac{\Gamma(\alpha)}{\beta^{\alpha}}\left(\frac{\alpha}{\beta}-\mu\right)^{-1}. $$ I recovered $g(x)$ by taking the inverse transform of both sides and then multiplying by $x^{1-\alpha}$. $$ \begin{aligned} g(x)% &=x^{1-\alpha}\mathcal L^{-1}\left\{\Gamma(\alpha)s^{-\alpha}\left(\frac{\alpha}{s}-\mu\right)^{-1}\right\}(x)\\ &=-\frac{x^{1-\alpha}}{\mu}\mathcal L^{-1}\left\{\Gamma(\alpha)s^{-\alpha}\left(1-\frac{\alpha/\mu}{s}\right)^{-1}\right\}(x). \end{aligned} $$ Using Bateman's Tables of Integral transforms, volume 1, $5.4.(9)$, this evaluates to $$ g(x)% =-\frac{1}{\mu}\Phi_2\left(1;\alpha;\frac{\alpha}{\mu}x\right), $$ where $$ \Phi_2(b_1,\dots,b_n;\gamma;z_1,\dots,z_n)=\sum_{m_1=0}^\infty \cdots\sum_{m_n=0}^\infty \frac{(b_1)_{m_1}\cdots (b_n)_{m_n}}{(\gamma)_{m_1+\cdots +m_n}m_1!\cdots m_n!}z_1^{m_1}\cdots z_n^{m_n} $$ is the hypergeometric function of $n$ variables. In this case we have a hypergeomatric function of a single variable; thus, $$ g(x)% =-\frac{1}{\mu}{_1}F_1\left(1;\alpha;\frac{\alpha}{\mu}x\right). $$

Unfortunately, this solution only yields sensible results if $\alpha/\beta<\mu$ (I have tried using some example parameters in MATLAB which demonstrates this). That said, what I am interested in is the case where $\alpha/\beta>\mu$. The formula in my table of integral transforms only has the restriction $\alpha>0$. Maybe there is an error? How can I get the solution to work for positive $\mu$?

We can check the solution which does seem to be correct. Using G&R formula $7.522.9$ we find $$ \begin{aligned} \operatorname Eg(X)% &=-\frac{\beta^{\alpha}}{\mu\Gamma(\alpha)}\int_0^\infty x^{\alpha-1}e^{-\beta x}{_1}F_1\left(1;\alpha;\frac{\alpha}{\mu}x\right).\,\mathrm dx\\ &=-\frac{1}{\mu}{_2}F_1\left({1,\alpha\atop\alpha};\frac{\alpha}{\beta\mu}\right)\\ &=-\frac{1}{\mu}{_1}F_0\left({1\atop -};\frac{\alpha}{\beta\mu}\right)\\ &=-\frac{1}{\mu}\left(1-\frac{\alpha}{\beta\mu}\right)^{-1}\\ &=\left(\frac{\alpha}{\beta}-\mu\right)^{-1}. \end{aligned} $$ So I am puzzled as to why this solution does not work for $\alpha/\beta>\mu$. One thing worth noticing is that when $\alpha/\beta<\mu$, the argument of the ${_1}F_0(1;-;\alpha/(\beta\mu))$ above is less than one and so the series defining it converges to $\left(1-\frac{\alpha}{\beta\mu}\right)^{-1}$ in the usual sense. For $\alpha/\beta\geq\mu$ the aruguement is greater than or equal to unity and the series defining the ${_1}F_0$ diverges; thus analytic continutation is used. Maybe this plays into the issue? Here is a test in MATLAB showing disagreement when $\alpha/\beta<\mu$:

alpha = sym(10);
beta = alpha/8;
mu = sym(5);

syms x g(x)
g(x) = -hypergeom(1,alpha,alpha*x/mu)/mu;
for i = 1:512
    X = gamrnd(double(alpha),double(1/beta));
    est(i) = vpa(g(X));
end

mean(est) = -11979.51
(alpha/beta-mu)^(-1) = 1/3

Answer 1

I have discovered that for my application the only useful result for $g(x)$ is one that does not contain $\beta$ and thus depends on $\alpha$, $x$, and $\mu$. This forces me to use Laplace transforms to solve the problem. I am posting this solution as a result of the comments from @reuns and @maxim which discuss using Mellin transforms. The equation numbers referenced to are from Bateman et. al. Tables of Integral Transforms.

To begin we have $$ \frac{\beta^\alpha}{\Gamma(\alpha)}\int_0^\infty x^{\alpha-1}e^{-\beta x}g(x)\,\mathrm dx=\left(\frac{\alpha}{\beta}-\mu\right)^{-1}. $$ In terms of the Mellin transform we write this integral equation as $$ \mathcal M\{e^{-\beta x}g(x)\}(\alpha)=\Gamma(\alpha)\beta^{-\alpha}\beta(\alpha-\beta\mu)^{-1}% \implies% g(x)=e^{\beta x}\mathcal M^{-1}\{h_1(\alpha)h_2(\alpha)\}(x), $$ where $h_1(\alpha)=\Gamma(\alpha)\beta^{-\alpha}$ and $h_2(\alpha)=\beta(\alpha-\beta\mu)^{-1}$. By $6.1.(14)$ we write the inverse transform as $$ g(x)=e^{\beta x}\int_0^\infty t^{-1}\mathcal M^{-1}\{h_1(\alpha)\}(x/t)\mathcal M^{-1}\{h_2(\alpha)\}(t)\,\mathrm dt. $$ By $6.3.(1)$ we find $\mathcal M^{-1}\{h_1(\alpha)\}(x)=e^{-\beta x}$. Furthermore, we want $\alpha/\beta>\mu$ so using $7.1.(3)$ we have $\mathcal M^{-1}\{h_2(\alpha)\}(x)=\beta x^{-\beta\mu}\mathbf 1_{(0,1)}(x)$. Substituting these results into the integral for $g(x)$ yields $$ g(x)=\beta e^{\beta x}\int_0^\infty t^{-1}e^{-\beta x/t}t^{-\beta\mu}\mathbf 1_{(0,1)}(t)\,\mathrm dt% =\beta e^{\beta x}\int_0^1 t^{-\beta\mu-1}e^{-\beta x/t}\,\mathrm dt. $$ Let $u=\beta x/t\implies t=\beta x/u$, $\mathrm dt=-\beta x/u^2\,\mathrm du$, then $$ g(x)=\beta (\beta x)^{-\beta\mu}e^{\beta x}\int_{\beta x}^\infty u^{\beta\mu-1}e^{-u}\,\mathrm du=\beta (\beta x)^{-\beta\mu}e^{\beta x}\Gamma(\beta\mu,\beta x). $$ Finally, using DLMF $8.5.3$ we express this result as $$ g(x) = \beta\, U(1,1+\beta\mu,\beta x), $$ where $U(a,b,z)$ is the confluent hypergeometric function of the second kind.

Here are some example parameters implemented in MATLAB:

mu = sym(1);
alpha = sym(10);
sigma = sym(sqrt(5));
beta = alpha/sigma^2;
syms x g(x)
g(x) = beta*kummerU(1,1+beta*mu,beta*x);
E = (alpha/beta-mu)^(-1);

for i = 1:512
    X = gamrnd(double(alpha),double(1/beta));
    est(i) = vpa(g(X));
end

vpa(E)
mean(est)

The results of the code show agreement with the derived result. Thanks to @reuns and @maxim for their thoughts.