# Solving an integral equation with inverse Laplace transform

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

• Let $h_b(x) = g(x) e^{-b x}$ and $H_b(a) = \int_0^\infty x^{a-1} h_b(x)dx$ Then $H_b(a) = b^{-a} \Gamma(a) \frac{b}{a-bu}$ implies $h_b$ is equal to the Mellin convolution $h_b(x) = \int_0^\infty u_b(y) v_b(x/y)dy$ where $\int_0^\infty u_b(x)x^{a-1}dx = b^{-a} \Gamma(a)$ and $\int_0^\infty v_b(x)x^{a-1}dx = \frac{b}{a-bu}$ so $u_b(x) = e^{-bx}$ or $u_b(x)= e^{-bx} - \sum_{k \le K} \frac{(-b)^k}{k!} x^k$ and $v_b(x) = -b x^{-bu} 1_{x < 1}$ or $v_b(x) = b x^{-bu} 1_{x > 1}$, the "or" depending on the strip of convergence of those Mellin transforms. – reuns Dec 22 '18 at 3:25
• I think some conditions might be missing here. For instance, $g(x) = (\alpha/\beta - \mu)^{-1}$ is a solution. – Maxim Dec 22 '18 at 21:30
• @Maxim you are right that that is a solution. That said, the goal is to come up with a function of the random variable that is an unbiased estimator of the r.h.s. I am not sure what other conditions could be specified other than the ones stated, I.e. a/b>u, a,b,u>0. Was there some conditions you had in mind that were omitted? – Aaron Hendrickson Dec 22 '18 at 22:01
• @Maxim I would also point out that the solution g(x) in my question works for u<0 but what I really need is one for when u is positive. I'm not sure how to adjust the approach to make that happen. – Aaron Hendrickson Dec 22 '18 at 22:04
• At the very least, $g(x)$ cannot be positive for all positive $x$. If it is and the integral converges for $\beta < \alpha/\mu$, we'll have convergence for larger values of $\beta$ as well. But the equation says the integral diverges for $\beta = \alpha/\mu$. The inverse Laplace transform takes $\beta$ to the right of all singularities of $F(\beta)$. – Maxim Dec 24 '18 at 14:27

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.