2
$\begingroup$

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 enter image description here

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.

$\endgroup$
6
  • $\begingroup$ @JimB This might be right up your alley. $\endgroup$
    – Carl
    Jul 14, 2022 at 21:54
  • $\begingroup$ mean = Assuming[a > 0 && b > 0 && α > 0 && β > 0 && b > β, Integrate[x PDF[Gdc, x], {x, 0, Infinity}]] // Expand gives the expected result. $\endgroup$
    – Bob Hanlon
    Jul 15, 2022 at 0:05
  • $\begingroup$ @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? $\endgroup$
    – Carl
    Jul 15, 2022 at 1:58
  • $\begingroup$ I posted my results with v13.1 on my Mac. I get these same results with v13.0 and v13.0.1 $\endgroup$
    – Bob Hanlon
    Jul 15, 2022 at 2:23
  • 1
    $\begingroup$ 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]. $\endgroup$
    – JimB
    Jul 15, 2022 at 3:43

1 Answer 1

4
$\begingroup$
$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]]

enter image description here

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]]

enter image description here

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]

enter image description here

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 + α/β *)
$\endgroup$
4
  • $\begingroup$ (+1) Thanks for confirming my problem, next step is finding out why, and what else is bonkers. $\endgroup$
    – Carl
    Jul 15, 2022 at 2:26
  • $\begingroup$ 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. $\endgroup$
    – Carl
    Jul 15, 2022 at 2:46
  • $\begingroup$ Wow! What a tour de force. I accept your answer and will put in a bug tag if you don't mind. $\endgroup$
    – Carl
    Jul 15, 2022 at 4:58
  • $\begingroup$ The bug is still present in version 13.2 (November 18, 2022) $\endgroup$
    – user58955
    Dec 15, 2022 at 15:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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