Integral involving a Gaussian and a fraction.

Question

This question is a generalization of An identity involving the incomplete Beta function. .

Let $x\ge 0$ and $\epsilon_\pm \in (1,\infty)$. We consider the following integral: \begin{equation} {\mathcal I}(x):= \int\limits_0^\infty \frac{ \eta^{\epsilon_+}}{B(\epsilon_-,\epsilon_+) \left( \eta + 1 \right)^{\epsilon_-+\epsilon_+}} \cdot \frac{\exp(-\frac{1}{2} x^2 \cdot \eta^2)}{\sqrt{2 \pi}} d\eta \end{equation}

Note that ${\mathcal I}(|x|)$ is the probability density function (pdf) of a random variable ${\mathcal R}:=(1-Z)/Z \cdot \Xi$ where $Z \sim Beta(\epsilon_+,\epsilon_-)$ and $\Xi \sim N(0,1)$ and $Z$ and $\Xi$ are independent.

Now I am going to present two different ways of computing the integral in question. Both methods are flawed for certain reasons, yet as it will turn out at the end, a sum of the results obtained by the two methods will give the correct result.

Method 1:

\begin{eqnarray} && \sqrt{2 \pi} B(\epsilon_+,\epsilon_-) {\mathcal I}(x)= \\ && \sum\limits_{p=0}^\infty \binom{-\epsilon_+-\epsilon_-}{p} \int\limits_0^\infty \eta^{-p-\epsilon_-} \exp(-\frac{x^2}{2} \eta^2) d\eta =\\ &&\frac{1}{2} \sum\limits_{p=0}^\infty \binom{-\epsilon_+-\epsilon_-}{p} (\sqrt{\frac{x^2}{2}})^{p+\epsilon_--1} \cdot \Gamma(-\frac{p+\epsilon_-}{2}+\frac{1}{2} ) =\\ &&\frac{1}{2} \sum\limits_{p=0}^\infty \binom{-\epsilon_+-\epsilon_-}{p} (\sqrt{\frac{x^2}{2}})^{p+\epsilon_--1} \cdot \frac{\pi}{\sin(\frac{1}{2} (1-\epsilon_-) \pi - p \pi/2) \Gamma(\frac{1}{2} (1+\epsilon_-+p) } =\\ &&x^{\epsilon_--1} \left(\right.\\ &&\left. \sqrt{2}^{-1-\epsilon_-} \Gamma(\frac{1}{2} -\frac{\epsilon_-}{2}) F_{2,2} \left[ \begin{array}{rr} \frac{1}{2}(\epsilon_++\epsilon_-) & \frac{1}{2}(\epsilon_++\epsilon_-+1) \\ \frac{1}{2} & \frac{1}{2}(1+\epsilon_-)\end{array} ; -\frac{x^2}{2} \right] - \right.\\ && \left. \sqrt{2}^{-2-\epsilon_-} \cdot x \cdot (\epsilon_++\epsilon_-) \Gamma(-\frac{\epsilon_-}{2}) F_{2,2} \left[ \begin{array}{rr} \frac{1}{2}(\epsilon_++\epsilon_-+1) & \frac{1}{2}(\epsilon_++\epsilon_-+2) \\ \frac{3}{2} & \frac{1}{2}(2+\epsilon_-)\end{array} ; -\frac{x^2}{2} \right] \right.\\ &&\left. \right) \end{eqnarray}

In the second line from the top we expanded the fraction in the integrand in a Taylor series and then we integrated term by term. Clearly this is wrong because neither the power series expansion converges on the whole positive half axis nor do the integrals in question exist. Yet if we forget about this and replace the divergent integral by the analytic continuation of the Gamma function we do get a convergent series (see third line from the top). Now in the forth line from the top we applied the reflection formula to the Gamma function and finally in the fifth line from the top we split the sum over $p$ into sums over even and odd values of $p$ and then we applied the duplication formula for the Gamma function to reduce the respective sums to hypergeometric functions.

Method 2: \begin{eqnarray} && \sqrt{2 \pi} B(\epsilon_+,\epsilon_-) {\mathcal I}(x)= \\ && \int\limits_0^1 v^{\epsilon_--2} (1-v)^{\epsilon_+} \exp(-\frac{x^2}{2} \frac{(1-v)^2}{v^2}) dv=\\ && \sum\limits_{p=0}^\infty \frac{(-\frac{x^2}{2})^p}{p!} \cdot \frac{\Gamma(\epsilon_--1-2 p) \Gamma(\epsilon_++1+2 p)}{\Gamma(\epsilon_++\epsilon_-)}\\ && B(\epsilon_+,\epsilon_-) \frac{\epsilon_+}{\epsilon_--1} F_{2,2} \left[ \begin{array}{rr} \frac{1+\epsilon_+}{2} & \frac{2+\epsilon_+}{2} \\ 1-\frac{\epsilon_-}{2} & \frac{3}{2} - \frac{\epsilon_-}{2}\end{array};-\frac{x^2}{2} \right] \end{eqnarray} In the first line from the top we substituted $v:=1/(1+\eta)$ and in the second line from the top we expanded the exponential in a power series and integrated term by term. Again this is not correct because the integrals in question do not exist if $p$ is large enough. Yet if we forget about that and replace the divergent integrals by the analytic continuation of the Gamma function then we obtain a convergent series.

Now if we add the two different results together we have:

\begin{eqnarray} &&{\mathcal I}(x):=\\ &&\frac{\epsilon_+}{\sqrt{2 \pi} (\epsilon_--1)} \cdot F_{2,2} \left[ \begin{array}{rr} \frac{1+\epsilon_+}{2} & \frac{2+\epsilon_+}{2} \\ 1-\frac{\epsilon_-}{2} & \frac{3}{2} - \frac{\epsilon_-}{2}\end{array};-\frac{x^2}{2} \right] +\\ && \frac{x^{\epsilon_--1}}{\sqrt{2 \pi} B(\epsilon_+,\epsilon_-)} \cdot \left(\right.\\ &&\left. \sqrt{2}^{-1-\epsilon_-} \Gamma(\frac{1}{2} -\frac{\epsilon_-}{2}) F_{2,2} \left[ \begin{array}{rr} \frac{1}{2}(\epsilon_++\epsilon_-) & \frac{1}{2}(\epsilon_++\epsilon_-+1) \\ \frac{1}{2} & \frac{1}{2}(1+\epsilon_-)\end{array} ; -\frac{x^2}{2} \right] - \right.\\ && \left. \sqrt{2}^{-2-\epsilon_-} \cdot x \cdot (\epsilon_++\epsilon_-) \Gamma(-\frac{\epsilon_-}{2}) F_{2,2} \left[ \begin{array}{rr} \frac{1}{2}(\epsilon_++\epsilon_-+1) & \frac{1}{2}(\epsilon_++\epsilon_-+2) \\ \frac{3}{2} & \frac{1}{2}(2+\epsilon_-)\end{array} ; -\frac{x^2}{2} \right] \right.\\ &&\left. \right) \end{eqnarray}

and it turns out that the result is correct as we check numerically below:

(*All together *)
{ep, em} = RandomReal[{1, 10}, 2, WorkingPrecision -> 50]; x = 
 RandomReal[{0, 5}, WorkingPrecision -> 50];
NIntegrate[(x^em eta^ep)/(x Beta[em, ep] (eta + x)^(em + ep))  
  Exp[-1/2 eta^2]/Sqrt[2 Pi], {eta, 0, Infinity}, 
 WorkingPrecision -> 20]
NIntegrate[
 eta^ep/( Beta[em, ep] (eta + 1)^(em + ep))  Exp[-1/2 x^2 eta^2]/Sqrt[
  2 Pi], {eta, 0, Infinity}, WorkingPrecision -> 20]

 ep /(Sqrt[2 Pi] (em - 1))
   HypergeometricPFQ[{1/2 + ep/2, 1 + ep/2}, {1 - em/2, 
    3/2 - em/2}, -(x^2/2)] + 
 x^em/(Sqrt[2 \[Pi]]
    x Beta[em, 
    ep]) (2^(1/2 (-1 - em))
      Gamma[1/2 - em/
       2] HypergeometricPFQ[{(ep + em)/2, (ep + em + 1)/2}, {1/2, 
       1/2 (1 + em)}, -x^2/2] -  
    2^(1/2 (-2 - em))
      x (ep + em) Gamma[-(em/
       2)] HypergeometricPFQ[{(ep + em + 2)/2, (ep + em + 1)/2}, {3/2,
        1/2 (2 + em)}, -x^2/2])

Update: Below I also plot the integral in question as a function of $x$ . Here we took $\epsilon_+=1.122$, $\epsilon_-=2.239$. In the top figure I plot the two different terms in the result , firstly the term being analytical in $x$ (Blue) and secondly the term being proportional toa fractional power of $x$ (Purple). In the bottom figure I plot the whole result (Green). We have:

SetOptions[ListPlot, LabelStyle -> {15, FontFamily -> "Arial"}, 
  BaseStyle -> {15, FontFamily -> "Bold"}];
{em, ep} = RandomReal[{1, 3}, 2];
f1[x_] := 
  ep/(Sqrt[2 Pi] (em - 1))
    HypergeometricPFQ[{1/2 + ep/2, 1 + ep/2}, {1 - em/2, 
     3/2 - em/2}, -(x^2/2)];
f2[x_] := x^(em - 1)/(Sqrt[2 \[Pi]] Beta[em, ep]) (
    2^(1/2 (-1 - em))
       Gamma[1/2 - em/
        2] HypergeometricPFQ[{(ep + em)/2, (ep + em + 1)/2}, {1/2, 
        1/2 (1 + em)}, -x^2/2] -  
     2^(1/2 (-2 - em))
       x (ep + em) Gamma[-(em/
        2)] HypergeometricPFQ[{(ep + em + 2)/2, (ep + em + 1)/2}, {3/
         2, 1/2 (2 + em)}, -x^2/2]);
mypdf[x_] := f1[x] + f2[x];
x = Array[# &, 100, {0, 5}];

pl1 = ListPlot[Transpose[{x, #}] & /@ {f1[x], f2[x]}, Joined :> True, 
   PlotMarkers -> Automatic, PlotRange -> All, ImageSize -> 700, 
   PlotLabel :> "ep,em=" <> ToString[{ep, em}]];
pl2 = ListPlot[Transpose[{x, #}] & /@ {mypdf[x]}, Joined :> True, 
   PlotMarkers -> Automatic, PlotRange -> All, PlotStyle -> Green, 
   ImageSize -> 700];
pl = GraphicsGrid[{{pl1}, {pl2}}];
Export["myGaussianIntegral.jpg", pl, "JPG"];
Import["myGaussianIntegral.jpg"]
x =.; NIntegrate[mypdf[x], {x, 0, Infinity}]

Now the obvious question is why is it that a sum of two "wrong" results gives something that turns out to be a correct result ?

Is there any other mathematically sound way of computing the integral in question ?