Mean not working with Gamma convolution distribution

Question

Using v 13.01 or 13.1: bug in integrate, see accepted answer below.

First off, the distribution Simplify[PDF[GammaDistribution[a, 1/b], x], Assumptions -> b > 0] $$\begin{array}{cc} \Bigg\{ & \begin{array}{cc} \dfrac{b^a x^{a-1} e^{-b x}}{\Gamma (a)} & x>0 \\ 0 & \text{True} \\ \end{array}\;\;, \\ \end{array}$$ which is an alternative form of the gamma distribution, when convolved with the same function but with parameters $\alpha$ and $\beta$ yields $$\text{Gdc}=\begin{array}{cc} \Bigg\{ & \begin{array}{cc} \dfrac{b^a \beta ^{\alpha } e^{b (-x)} x^{a+\alpha -1} \, _1F_1(\alpha ;a+\alpha ;x (b-\beta ))}{\Gamma (a+\alpha )} & x>0 \\ 0 & \text{True} \\ \end{array} \\ \end{array}\;\;,$$ where $_1F_1$ is Hypergeometric1F1. For documentation and application of the Gdc density function see https://ejnmmiphys.springeropen.com/articles/10.1186/s40658-016-0166-z

In case someone wants to try this I have replaced $\alpha$ and $\beta$ that do not render well on this site with A and B here (x^(-1+a+A) b^a E^(-x b) B^A Hypergeometric1F1[A,a+A,x (b-B)])/Gamma[a+A].

Next, please recall that the mean value of a convolution of two density functions is the sum of the mean values of each function (provided they exist), https://stats.stackexchange.com/q/342023/99274, or for this parameterization as each mean takes the form $\alpha\,/\beta$ (see https://en.wikipedia.org/wiki/Gamma_distribution), $$\bar{x}=\dfrac{a}{b}+\dfrac{\alpha}{\beta}\;\;.$$ Couldn't be simpler, right? However, if we define $$\text{Gdc}=\text{ProbabilityDistribution}\left[ \begin{array}{cc} \Bigg\{ & \begin{array}{cc} \frac{x^{a+\alpha -1} b^a e^{-x b} \beta ^{\alpha } \, _1F_1(\alpha ;a+\alpha ;x (b-\beta ))}{\Gamma (a+\alpha )} & x>0 \\ 0 & \text{True} \\ \end{array} \\ \end{array} ,\{x,0,\infty \},\text{Assumptions}\to \{a>0,b>0,\alpha >0,\beta >0,b>\beta ,x\in \mathbb{R},a\in \mathbb{R},b\in \mathbb{R},\alpha \in \mathbb{R},\beta \in \mathbb{R}\}\right]\;\;,$$ where $b>\beta$ without loss of generality, and then ask for Mean[Gdc] without assigning values to the parameters we get of all things $$\frac{(-1)^{-\alpha } \cos (\pi \alpha ) (a \beta +\alpha b)}{b \beta }\;\;,$$ which is obviously incorrect as the mean value of definite positive values, as required to define a density, can be neither complex valued nor negative.

However, if we put in values for the parameters of $\left(a=1;b=1;\alpha =\frac{1}{2};\beta =\frac{1}{10};\right)$, which, when plotted Plot[PDF[Gdc,x],{x,0,10},PlotRange->Full] looks like

we get the correct values of 6 for both Mean[Gdc] and $\dfrac{a}{b}+\dfrac{\alpha }{\beta }$, but 0 for $$\dfrac{(-1)^{-\alpha } \cos (\pi \alpha ) (a \beta +\alpha b) }{b \beta }$$

So, obviously, something is wrong, but what? This looks like a bug of some sort. Any insight into what this means would be helpful, in particular, I would like to be able to attempt using FindDistributionParameters and am uncertain what I would be getting from that.

Note As pointed out by @JimB in his comment below, when $\alpha$ is exactly a positive integer, then and only then does

$$\frac{(-1)^{-\alpha } \cos (\pi \alpha ) (a \beta +\alpha b)}{b \beta }=\dfrac{a}{b}+\dfrac{\alpha}{\beta}\iff\alpha\in\mathbb{Z}^+\;,$$ which is vanishingly improbable for a real valued function.

Answer 1

$Version

(* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *)

Clear["Global`*"]

dist1 = GammaDistribution[a, 1/b];

dpa1 = DistributionParameterAssumptions[dist1] // Simplify

a > 0 && b > 0

dist2 = GammaDistribution[α, 1/β];

dpa2 = DistributionParameterAssumptions[dist2] // Simplify;

dpa = dpa1 && dpa2 && b > β;

pdf = Assuming[dpa,
  Convolve[PDF[dist2, y], PDF[dist1, y], y, x]]

Gdc = ProbabilityDistribution[
   pdf, {x, 0, Infinity},
   Assumptions -> dpa];

As you indicated, Mean fails

Mean[Gdc]

(* (E^(-2 I π α) (b α + a β))/(b β) *)

However, direct calculation provides the expected result.

mean = Assuming[dpa,
   Integrate[x PDF[Gdc, x], {x, 0, Infinity}]] // Expand

(* a/b + α/β *)

EDIT: Note that the integral of pdf does not return unity

Assuming[dpa, Integrate[pdf, {x, 0, Infinity}]]

(* E^(-2 I π α) *)

I believe that this is a bug in Integrate. Express the pdf using WhittakerM

altRep = Entity["MathematicalFunction", "Hypergeometric1F1"][
   "AlternativeRepresentations"][[5]]

altRep[a, b, z] // Activate // FullSimplify

(* True *)

pdf is then

pdf2 = Assuming[dpa,
  (pdf // FunctionExpand) /.
    Hypergeometric1F1[a_, b_, z_] :>
     E^(z/2) z^(-b/2) WhittakerM[(b - 2 a)/2, (b - 1)/2, z] // 
  Simplify]

The total probability is now correctly unity.

Assuming[dpa, Integrate[pdf2, {x, 0, Infinity}]]

(* 1 *)

Gdc2 = ProbabilityDistribution[pdf2, {x, 0, Infinity}, 
  Assumptions -> dpa];

With this representation both Mean and direct calculation provide the expected result.

Mean[Gdc2] // Expand

(* a/b + α/β *)

mean = Assuming[dpa, Integrate[x PDF[Gdc2, x], 
  {x, 0, Infinity}]] // Expand

(* a/b + α/β *)