# Mean not working with Gamma convolution distribution

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.

• @JimB This might be right up your alley.
– Carl
Commented Jul 14, 2022 at 21:54
• mean = Assuming[a > 0 && b > 0 && α > 0 && β > 0 && b > β, Integrate[x PDF[Gdc, x], {x, 0, Infinity}]] // Expand gives the expected result. Commented Jul 15, 2022 at 0:05
• @BobHanlon Hummm. I get $\frac{(-1)^{-\alpha } a \cos (\pi \alpha )}{b}+\frac{(-1)^{-\alpha } \alpha \cos (\pi \alpha )}{\beta }$ in 13.0 when I try that. A bug perhaps?
– Carl
Commented Jul 15, 2022 at 1:58
• I posted my results with v13.1 on my Mac. I get these same results with v13.0 and v13.0.1 Commented Jul 15, 2022 at 2:23
• For v13.01 on Windows 10 I get $\frac{e^{-2 i \pi \alpha } (a \beta +\alpha b)}{b \beta }$ for the mean. The correct mean seems only to work for integer values of $\alpha$: FullSimplify[(E^(-2 I \[Pi] \[Alpha]) (b \[Alpha] + a \[Beta]))/(b \[Beta]), Assumptions -> \[Alpha] \[Element] Integers].
– JimB
Commented Jul 15, 2022 at 3:43

\$Version

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

Clear["Global*"]

dpa1 = DistributionParameterAssumptions[dist1] // Simplify

a > 0 && b > 0

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 + α/β *)

• (+1) Thanks for confirming my problem, next step is finding out why, and what else is bonkers.
– Carl
Commented Jul 15, 2022 at 2:26
• I also confirm your solution, i.e., mean = Assuming[dpa, Integrate[x PDF[Gdc, x], {x, 0, Infinity}]] // Expand (* a/b + α/β *) Still want to know the extent of the problem.
– Carl
Commented Jul 15, 2022 at 2:46
• Wow! What a tour de force. I accept your answer and will put in a bug` tag if you don't mind.
– Carl
Commented Jul 15, 2022 at 4:58
• The bug is still present in version 13.2 (November 18, 2022) Commented Dec 15, 2022 at 15:38