Ytukyg

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}-muright)^{-1},
$$
where $alpha/beta>mu$. Note that such a solution would yield an unbiased estimator for $left(frac{alpha}{beta}-muright)^{-1}$, i.e. if $Xsimoperatorname{Gamma}(alpha,beta)$ then $operatorname Eg(X)=left(frac{alpha}{beta}-muright)^{-1}$. I tried solving this with an inverse Laplace transform by writing
$$
mathcal Lleft{x^{alpha-1}g(x)right}(beta)=frac{Gamma(alpha)}{beta^{alpha}}left(frac{alpha}{beta}-muright)^{-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}-muright)^{-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_2left(1;alpha;frac{alpha}{mu}xright),
$$
where
$$
Phi_2(b_1,dots,b_n;gamma;z_1,dots,z_n)=sum_{m_1=0}^infty cdotssum_{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_1left(1;alpha;frac{alpha}{mu}xright).
$$

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}}{muGamma(alpha)}int_0^infty x^{alpha-1}e^{-beta x}{_1}F_1left(1;alpha;frac{alpha}{mu}xright).,mathrm dx\
&=-frac{1}{mu}{_2}F_1left({1,alphaatopalpha};frac{alpha}{betamu}right)\
&=-frac{1}{mu}{_1}F_0left({1atop -};frac{alpha}{betamu}right)\
&=-frac{1}{mu}left(1-frac{alpha}{betamu}right)^{-1}\
&=left(frac{alpha}{beta}-muright)^{-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/(betamu))$ above is less than one and so the series defining it converges to $left(1-frac{alpha}{betamu}right)^{-1}$ in the usual sense. For $alpha/betageqmu$ 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

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}-muright)^{-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-betamu)^{-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-betamu)^{-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^{-betamu}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^{-betamu}mathbf 1_{(0,1)}(t),mathrm dt%
=beta e^{beta x}int_0^1 t^{-betamu-1}e^{-beta x/t},mathrm dt.
$$
Let $u=beta x/timplies t=beta x/u$, $mathrm dt=-beta x/u^2,mathrm du$, then
$$
g(x)=beta (beta x)^{-betamu}e^{beta x}int_{beta x}^infty u^{betamu-1}e^{-u},mathrm du=beta (beta x)^{-betamu}e^{beta x}Gamma(betamu,beta x).
$$
Finally, using DLMF $8.5.3$ we express this result as
$$
g(x) = beta, U(1,1+betamu,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.