## 1 Introduction
This is known in statistics as the two-sample test of the hypothesis \(H_{\rm o}:\ E(X)=E(Y)\). Such a technique, however, is no simple task when the two averages are calculated from samples drawn from gamma populations because for the procedure to be applicable, the probability distribution of \(\overline{X}-\overline{Y}\) must be known under \(H_{\rm o}\) so that the important judgment can be made on whether \(\overline{X}-\overline{Y}\) has happened to take on a value that is very unlikely under \(H_{\rm o}\). If \(X=(X_{1},X_{2},\ldots,X_{n})\) is a random sample from the \(\Gamma(k=shape,\theta=scale)\) gamma distribution, the characteristic function of each \(X_{j}\) is  
\[\psi_{X_{j}}(t)=(1-i\theta t)^{-k}, \tag{1}\]  
and so  
\[\psi_{\sum X_{j}}(t)=(1-i\theta t)^{-kn}, \tag{2}\]  
\[\psi_{\overline{X}}(t)=(1-it\theta/n)^{-kn}, \tag{3}\]  
as is known from the theory of characteristic functions [2; 3; 4]. In other words, \(\overline{X}\) follows the distribution \(\Gamma(kn,\theta/n)\). Taking the analogous result for \(\overline{Y}\), \(\overline{X}-\overline{Y}\) is the difference of two independent gamma-distributed random variables if the two samples are independent. When the assumption of independence holds, it is straightforward to derive the characteristic function of the difference, but it is not clear anymore what distribution it belongs to. The problem of how the difference is distributed under \(H_{\rm o}\) could also be solved by deriving the density of \(-\overline{Y}\), applying the transformation theorem [5], since the general form of the density of \(\overline{Y}\) is known, and then one can use the convolution theorem [3; 6] to derive the density of \(\overline{X}+(-\overline{Y})\) and hence its distribution. But the result for this density in the form of an integral containing two nontrivial densities and their multiplication is too complex to evaluate in practice, even in theory. This means it is difficult to derive an exact test for the stated purpose since the distribution of \(\overline{X}-\overline{Y}\) under \(H_{\rm o}\) derived from its complex density is not very clear. On the other hand, it is certainly possible to analyze the distribution of \(\overline{X}-\overline{Y}\) under \(H_{\rm o}\) through simulations and, based on the simulation results, to formulate recommendations on how \(H_{\rm o}:E(X)=E(Y)\) could be tested. This is especially beneficial when dealing with smaller data samples since for very large data sets, the problem of evaluating the significance of the difference \(\overline{X}-\overline{Y}\) can be resolved with the central limit theorems and the associated higher-precision normality approximations. The parameters of the gamma distribution have been under the spotlight for quite some time now given the increasing practical importance of this stochastic model. Yet the interest has not always been focused on two-sample testing, and when it was, as well as when it was not, the scientific output either bore greater complexity for practitioners or did not always deliver precise results simply because the form of the gamma distribution was too involved to do so. To provide some excursion into the past in this regard, [PERSON] and [PERSON] [7] proposed an approximate test for the mean of a gamma distribution with both its parameters unknown, which is a different problem than the one analyzed in this text. [PERSON] and [PERSON] [8] then proposed an approximate test for testing the equality of the scale parameters of two gamma distributions with unknown and equal shape parameters. [PERSON] et al. [9] also proposed an approximate test for testing the equality of the means of two gamma distributions with unknown and unequal shape parameters. In the 1990s, [PERSON] et al. [10] suggested an asymptotic test, that is, only an approximate procedure for many small-sample situations for the parameters of two gamma distributions. [PERSON] et al.