# How to calculate the PDF of product of two random variables from generalized gamma distributions?

Namely, we want to find the explicit formula of PDF of double Generalized Gamma distribution.

MellinConvolve[
x,y]

TransformedDistribution[
x*y
, {
}
] // PDF[#, t] &

pdfx = (E^(-((x^k1 φ1)/Ω1))  k1 x^(-1 +  k1 φ1) (φ1/Ω1)^φ1) /Gamma[φ1]

pdfy = (E^(-((y^k2  φ2)/Ω2)) k2  y^(-1 + k2  φ2)  (φ2/Ω2)^φ2)/Gamma[φ2] /. {y -> z/x};

Assuming[
{k1 > 0, φ1 > 0, Ω1 > 0},
Integrate[
pdfx*pdfy*1/x
, {x, 0, Infinity}
, PrincipalValue -> True
]
]


None of them can produce explicit result.

The expected result is (I have VERIFIED some test cases under some parameter settings, it's true):

$$f_T(t)=\frac{k_2 \alpha \beta^{\varphi_1-1 / 2} \alpha^{\varphi_2-1 / 2}(2 \pi)^{1-(\alpha+\beta) / 2}}{t \Gamma\left(\varphi_1\right) \Gamma\left(\varphi_2\right)} \times G_{q, 0}^{0, q}\left[\left.\frac{\alpha^\alpha \beta^\beta \Omega_1^\beta \Omega_2^\alpha}{t^{k_2 \alpha} \varphi_1^\beta \varphi_2^\alpha} \right\rvert\, \begin{array}{c} \Delta(\beta,1-\varphi_1),\Delta(\alpha,1-\varphi_2) & - \\ - & - \end{array}\right]$$

Where $$q=\alpha+\beta$$, $$\alpha,\beta \in \mathbb{N}^+$$ and $$\alpha/\beta=k_1/k_2$$; $$\Delta(m, n) \triangleq \frac{n}{m}, \frac{n+1}{m}, \cdots, \frac{n+m-1}{m}$$

On[Assert];
Assert[ResourceFunction["RationalQ"][Rationalize[k1/k2]]]
MyDeltaSign[m_Integer /; (m > 0), n_] :=
Table[(n + j)/m, {j, 0, m - 1}]
{\[Alpha], \[Beta]} = Rationalize[k1/k2] // NumeratorDenominator
k2*\[Alpha]*\[Beta]^(\[CurlyPhi]1 - 1/2)*\[Alpha]^(\[CurlyPhi]2 -
1/2)*(2*Pi)^(1 - (\[Alpha] + \[Beta])/2)/(t*Gamma[\[CurlyPhi]1]*
Gamma[\[CurlyPhi]2])*
MeijerG[{MyDeltaSign[\[Beta], 1 - \[CurlyPhi]1]~Join~
MyDeltaSign[\[Alpha],
1 - \[CurlyPhi]2], {}}, {{}, {}}, (\[Alpha]^\[Alpha]*\[Beta]^\
\[Beta]*\[CapitalOmega]1^\[Beta]*\[CapitalOmega]2^\[Alpha])/(t^(k2*\
\[Alpha])*\[CurlyPhi]1^\[Beta]*\[CurlyPhi]2^\[Alpha])]


Rewrite the above integral expression in FoxH: (firsly replace /. {x :> \[Tau]^(1/k1)}, then using formula (07.34.21.0012.01) from MeijerG.pdf)

k2 (\[CurlyPhi]1/\[CapitalOmega]1)^(1/k1) (\[CurlyPhi]2/\[CapitalOmega]2)^(1/k2) Gamma[\[CurlyPhi]1] Gamma[\[CurlyPhi]2]*FoxH[
(*Upper List*){
(*Upper Front list*){{1-(-(1/k1)+\[CurlyPhi]1), -k2/k1}},
(*Upper Rear List*){}
},
(*Lower List*){
(*Lower Front List*){{-(1/k2)+\[CurlyPhi]2, 1}},
(*Lower Rear List*){}
},
(t^k2 \[CurlyPhi]2 (\[CurlyPhi]1/\[CapitalOmega]1)^(k2/k1))/\[CapitalOmega]2]


• Making your code readable helps people understand what you are trying to do. Please elaborate on what you mean by "None of them work.". What is varI? Commented Jun 4 at 9:17
• The reference you cite states the steps taken (in words rather than explicit mathematical steps) to obtain the pdf. One of the steps consists of a change in variables. Maybe the experimental function IntegrateChangeVarianbles might help.
– JimB
Commented Jun 5 at 3:06

Even for concrete values of parameters the result is in terms of  MeijerG:

\[CurlyPhi]1 = 1; \[CurlyPhi]2 = 2; k1 = 1; k2 = 2; \[CapitalOmega]1
= 3; \[CapitalOmega]2 = 4; pr = TransformedDistribution[x*y,
PDF[pr, t]


Piecewise[ {{(t^3*MeijerG[{{}, {}}, {{-3/2, -1, 0}, {}}, t^2/72])/ (2592*Sqrt[Pi]), t > 0}}, 0]

If we help Mathematica a little, we will achieve the goal and obtain the solution to this integral.

I use Mathematica version 13.3, because version 14.0 have some Bugs in FunctionExpand to calculate a numeric values FoxH function see here and here

You want calculate this integral:

$$\int_0^{\infty } \frac{e^{-\frac{x^{\text{k1}} \text{\varphi 1}}{\text{\Omega 1}}-\frac{\left(\frac{t}{x}\right)^{\text{k2}} \text{\varphi 2}}{\text{\Omega 2}}} \text{k1} \text{k2} x^{-2+\text{k1} \text{\varphi 1}} \left(\frac{t}{x}\right)^{-1+\text{k2} \text{\varphi 2}} \left(\frac{\text{\varphi 1}}{\text{\Omega 1}}\right)^{\text{\varphi 1}} \left(\frac{\text{\varphi 2}}{\text{\Omega 2}}\right)^{\text{\varphi 2}} \Gamma (\text{\varphi 1})}{\Gamma (\text{\varphi 2})} \, dx$$

Using Mellin Transform and Inverse Mellin Transform then:

\$Version
(*"13.3.0 for Microsoft Windows (64-bit) (June 3, 2023)"*)

F1 = MellinTransform[(E^(-A*(x^k1 \[CurlyPhi]1)/\[CapitalOmega]1 - (
t^k2 x^-k2 \[CurlyPhi]2)/\[CapitalOmega]2)
k1 k2 t^(-1 + k2 \[CurlyPhi]2)
x^(-1 + k1 \[CurlyPhi]1 -
k2 \[CurlyPhi]2) \[CurlyPhi]1^\[CurlyPhi]1 \[CurlyPhi]2^\
\[CurlyPhi]2 \[CapitalOmega]1^-\[CurlyPhi]1 \[CapitalOmega]2^-\
\[CurlyPhi]2 Gamma[\[CurlyPhi]1])/Gamma[\[CurlyPhi]2], A, s] //
PowerExpand(*Where: A=1*)

(*(E^(-((t^k2 x^-k2 \[CurlyPhi]2)/\[CapitalOmega]2)) k1 k2 t^(-1 +
k2 \[CurlyPhi]2) x^(-1 - k1 s + k1 \[CurlyPhi]1 -
k2 \[CurlyPhi]2) \[CurlyPhi]1^(-s + \[CurlyPhi]1) \[CurlyPhi]2^\
\[CurlyPhi]2 \[CapitalOmega]1^(
s - \[CurlyPhi]1) \[CapitalOmega]2^-\[CurlyPhi]2 Gamma[s] Gamma[\[CurlyPhi]1])/Gamma[\[CurlyPhi]2]*)

F2 = Integrate[F1, {x, 0, \[Infinity]}, Assumptions -> {\[CapitalOmega]2 > 0, \[CapitalOmega]1 > 0, k1 > 0,
k2 > 0, \[CurlyPhi]1 > 0, \[CurlyPhi]2 > 0, t > 0}]

(*ConditionalExpression[(k1 t^(-1 +
k2 \[CurlyPhi]2) (\[CurlyPhi]1/\[CapitalOmega]1)^(-s + \
\[CurlyPhi]1) (\[CurlyPhi]2/\[CapitalOmega]2)^\[CurlyPhi]2 ((
t^k2 \[CurlyPhi]2)/\[CapitalOmega]2)^(-((
k1 s - k1 \[CurlyPhi]1 + k2 \[CurlyPhi]2)/k2))
Gamma[s] Gamma[\[CurlyPhi]1] Gamma[(k1 (s - \[CurlyPhi]1))/
k2 + \[CurlyPhi]2])/Gamma[\[CurlyPhi]2],
k1 \[CurlyPhi]1 < k2 \[CurlyPhi]2 + k1 Re[s]]*)

InverseMellinTransform[F2[[1]], s, A] /. A -> 1(*Can't compute.Probably in newer version like 14.1 can computed*)


Using another method to compute Inverse Mellin Transform:

 Sum[Assuming[m \[Element] Integers && m >= 0,
Residue[F2[[1]]*A^(-s) /. A -> 1, {s, -m}]], {m, 0,
Infinity}](*Can't compute.Probably in newer version like 14.1 I hope can computed*)


$$\sum _{m=0}^{\infty } \frac{(-1)^m \text{k1} t^{-1+\text{k2} \text{\varphi 2}} \left(\frac{\text{\varphi 1}}{\text{\Omega 1}}\right)^{m+\text{\varphi 1}} \left(\frac{\text{\varphi 2}}{\text{\Omega 2}}\right)^{\text{\varphi 2}} \left(\frac{t^{\text{k2}} \text{\varphi 2}}{\text{\Omega 2}}\right)^{\frac{\text{k1} m}{\text{k2}}+\frac{\text{k1} \text{\varphi 1}}{\text{k2}}-\text{\varphi 2}} \Gamma (\text{\varphi 1}) \Gamma \left(-\frac{\text{k1} (m+\text{\varphi 1})}{\text{k2}}+\text{\varphi 2}\right)}{m! \Gamma (\text{\varphi 2})}=\frac{\text{k1} t^{-1+\text{k2} \text{\varphi 2}} \left(\frac{\text{\varphi 1}}{\text{\Omega 1}}\right)^{\text{\varphi 1}} \left(\frac{\text{\varphi 2}}{\text{\Omega 2}}\right)^{\text{\varphi 2}} \left(\frac{t^{\text{k2}} \text{\varphi 2}}{\text{\Omega 2}}\right)^{\frac{\text{k1} \text{\varphi 1}}{\text{k2}}-\text{\varphi 2}} \sum _{m=0}^{\infty } \frac{\left(-\frac{\text{\varphi 1} \left(\frac{t^{\text{k2}} \text{\varphi 2}}{\text{\Omega 2}}\right)^{\text{k1}/\text{k2}}}{\text{\Omega 1}}\right)^m \Gamma \left(-\frac{\text{k1} (m+\text{\varphi 1})}{\text{k2}}+\text{\varphi 2}\right)}{m!}}{\Gamma (\text{\varphi 2})}$$

Using:

 Sum[ Gamma[a1 + A1 m]*z^m/m!, {k, 0, Infinity}] == FoxH[{{{1 - a1, A1}}, {}}, {{{0, 1}}, {}}, -z]
(*And again, it's impossible to compute, it's a shame.
I have the impression that Mathematica can't do anything.*)


Then we have the Solution expressed by FoxH function:

 F3=(k1*t^(-1 + k2*\[CurlyPhi]2)*(\[CurlyPhi]1/\[CapitalOmega]1)^\[CurlyPhi]1*(\[CurlyPhi]2/\[CapitalOmega]2)^\[CurlyPhi]2*((t^k2*\[CurlyPhi]2)/\[CapitalOmega]2)^((k1*\[CurlyPhi]1)/k2 - \[CurlyPhi]2)*FoxH[{{{1 + (k1*\[CurlyPhi]1)/k2 - \[CurlyPhi]2, -(k1/k2)}}, {}}, {{{0, 1}}, {}}, (\[CurlyPhi]1*((t^k2*\[CurlyPhi]2)/\[CapitalOmega]2)^(k1/k2))/\[CapitalOmega]1]*Gamma[\[CurlyPhi]1])/Gamma[\[CurlyPhi]2]

\[CurlyPhi]1 = 1; \[CurlyPhi]2 = 2; k1 = 1; k2 = 2;
\[CapitalOmega]1 = 3; \[CapitalOmega]2 = 4; t = 1;

F3 // FunctionExpand // N
(*0.162355*)

pr = TransformedDistribution[x*y, {x \[Distributed]

Mathematica 13.3 or 14.0 can't convert FoxH to MeijerG function!