Let $X \sim \text{Gamma}(a,b)$ with pdf $f(x)$:
and $Y \sim \text{Gamma}(\alpha,\beta)$ be independent with pdf $g(y)$:
Then, the pdf of the product $Z = X Y$ can be obtained as $h(z)$:
where I am using the TransformProduct function from mathStatica/Mathematica to do the nitty-gritties, and BesselK[n,z] denotes the modified Bessel function of the second kind. This is much simpler than requiring MeijerG functions. I should note that I am one of the authors of the software function used.
Quick Monte Carlo check
- against theoretical solution derived above when $a =2$, $b = 3$, $\alpha = 4$, and $\beta = 1.1$
Looks fine :)
If you know the mean is $\mu$ and the standard deviation is $\sigma$, then the shape parameter of a Gamma distribution is $\dfrac{\mu^2}{\sigma^2}$ and the scale parameter is $\dfrac{\sigma^2}{\mu}$, making the corresponding rate parameter $\dfrac{\mu}{\sigma^2}$
As an illustration of what is possible, suppose you knew that the mean is $40$ and you had an interval of $[30,50]$ representing about $2$ standard deviations either side of the mean
Then the standard deviation is about $\frac{50-40}{2}=5$, and the variance is therefore about $5^2=25$
For a Gamma distribution with shape parameter $k$ and scale parameter $\theta$, the mean would be $k\theta$ and the variance $k\theta^2$, suggesting with these numbers that $\theta \approx \frac{25}{40} = 0.625$ (equivalent to a rate of $1.6$) and $k \approx \frac{40^2}{25}=64$
As a check, we can look at the corresponding interval for these parameters in R
> pgamma(50,shape=64,scale=0.625) - pgamma(30,shape=64,scale=0.625)
[1] 0.9553145
> c(qgamma(0.025,shape=64,scale=0.625),qgamma(0.975,shape=64,scale=0.625))
[1] 30.80487 50.37773
which shows this approach is not exact, but is not that far away.
$k=59.3749$ and $\theta=0.66312$ would get you closer to the confidence interval with $2.5\%$ each side but at the cost (due to the asymmetry of the Gamma distribution) of a corresponding mean of $39.372$ rather than $40$
Best Answer
I think the mistake is in your code. It looks like you tried to square with
*2instead of^2. Maybe you were thinking of Python, which would use**2.