5
$\begingroup$

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
$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – reuns Dec 22 '18 at 3:25
  • 1
    $\begingroup$ I think some conditions might be missing here. For instance, $g(x) = (\alpha/\beta - \mu)^{-1}$ is a solution. $\endgroup$ – Maxim Dec 22 '18 at 21:30
  • $\begingroup$ @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? $\endgroup$ – Aaron Hendrickson Dec 22 '18 at 22:01
  • $\begingroup$ @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. $\endgroup$ – Aaron Hendrickson Dec 22 '18 at 22:04
  • 1
    $\begingroup$ 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)$. $\endgroup$ – Maxim Dec 24 '18 at 14:27
-1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.