Large Sample Comparison of Parameter Estimates in Gamma Raindrop Distributions

: Raindrop size distributions have been characterized through the gamma family. Over the years, quite a few estimates of these gamma parameters have been proposed. The natural question for the practitioner, then, is what estimation procedure should be used. We provide guidance in answering this question when a large sample size ( > 2000 drops) of accurately measured drops is available. Seven estimation procedures from the literature: ﬁve method of moments procedures, maximum likelihood, and a pseudo maximum likelihood procedure, were examined. We show that the two maximum likelihood procedures provide the best precision (lowest variance) in estimating the gamma parameters. Method of moments procedures involving higher-order moments, on the other hand, give rise to poor precision (high variance) in estimating these parameters. A technique called the delta method assisted in our comparison of these various estimation procedures.


Introduction
Knowledge of the raindrop size distribution (RSD) is essential in the retrieval of rainfall properties using radar remote sensing techniques and in understanding and describing the microphysical processes involved in the formation of precipitation. In particular, gamma RSDs have been used, among other investigations, to examine the time evolution of rainfall [1], characterize and quantify relations between such quantities as rainfall rate, radar reflectivity, and liquid water content [2][3][4][5], and help distinguish between stratiform and convective precipitation [6]. Throughout this article, we suppose that the RSD follows a gamma distribution. While the gamma distribution is not always appropriate, it often is (see e.g., [7] and [8]). Given a sample of raindrops, a variety of estimates have been proposed over the years for how this data may be used to estimate the gamma parameters. Among these, method of moments estimates are the most popular as the various moments may be understood as physical quantities [9]. For example, total drop surface area is related to the second moment, total drop volume or liquid water content to the third moment, and the reflectivity factor to the sixth moment. Another estimation procedure is maximum likelihood [10][11][12][13]. In this article, we compared seven different sets of gamma parameter estimates: five method of moments estimates [1][2][3][4][5][6]9,12], maximum likelihood estimates [10][11][12][13], and pseudo-maximum likelihood estimates [14], and did so only for the case of a "large" (>2000) sample of drop sizes, which were accurately measured across the full spectrum of size. All disdrometers are limited by truncation or quantization effects, and hence the ideal RSD, which is assumed in this article, will never be measured. However, a close approximation may be feasible by combining multiple disdrometers that enable sampling of both small and large drops. While not yet widely available, the technology to measure drop size quite accurately across nearly the full drop size 2 of 11 spectrum is now available and will continue to improve over time. Presently, meteorological particle spectrometers may be used to measure small drop size (100 microns to 1.5 mm) with a resolution of 50 microns [15]. Additionally, third generation 2D-video disdrometers may be used to measure larger drops (0.7 mm and larger) with a resolution of 170 microns [15,16]. The large sample estimate comparisons in this article, then, may be viewed as exact in the presence of perfect drop size information or as a near approximation using the best current technology for recording drop size.
As it turns out, all seven sets of gamma parameter estimates presented below are correct on average in the case of large sample sizes; that is, they are unbiased. How do we decide which of the seven estimation procedures should be used? To settle this issue, we could look at the variabilities, specifically the variances (or their square roots, the standard deviations), to help us choose between estimates. Small variability, of course, is desirable. A technique sometimes referred to as the delta method from the engineering and statistical literature [17][18][19] was used to establish the lack of bias and the variances in all but the case of maximum likelihood.

Raindrop Size Distribution and Parameter Estimates
Use of the gamma distribution was proposed by [4,20], among others, as it often gives an appropriate description of the natural variations of the observed RSDs. For the present study, we represent the raindrop sizes as a gamma distribution function where n(D) represents the number of raindrops per unit diameter interval and per unit volume of air and the gamma density is given by following the notation in [21]. Here, N T is the total drop number concentration; µ is the shape parameter; and λ is the rate parameter (in mm −1 if drop size is in mm). Additionally, Γ refers to the gamma function which may be thought of as a generalized factorial function since Γ(n) = (n − 1)! for n a positive integer. Assuming a gamma density for atmospheric samples is equivalent to assuming a gamma density for surface samples, albeit with modified shape and rate parameters, under standard models of terminal raindrop fall velocities as a function of size. Further details may be found in [22] (Section 8).
In the maximum likelihood (ML) and pseudo-maximum likelihood gamma parameter estimates to follow, we will need to refer to the digamma function as well as its derivative ψ (x) = d dx ψ(x), sometimes referred to as the trigamma function (see, e.g., [23]). The gamma parameter estimates we compared are listed in Table 1. As a number of these estimates are obtained by the method of moments (MM) some moment notation is appropriate. Given measured raindrop size values D 1 , D 2 , . . . , D n let denote the j th sample moment with n the total number of raindrops. In Table 1, the MMrst notation indicates the method of moments is implemented using moments of order r, s, t, which will have distinct values from 0-6. In what follows, we will use m 1 and D (the sample mean raindrop diameter) Atmosphere 2020, 11, 333 3 of 11 interchangeably. The weighting of moments using division by n above, by the way, is arbitrary for the method of moments estimates to be discussed. So, for example, volume-weighted moments may be used instead. Table 1. Selected gamma shape and rate parameter estimates from the literature. Corresponding references appear in the first column. The estimates in all but the last two rows make use of an auxiliary variable 'a' in the last column. The last two rows do not make use of an auxiliary variable a so this column is filled with an asterisk.

Approach
Gamma Parameter Estimates â µλ Various combinations of three of these moments have commonly been used by atmospheric scientists to estimate the gamma parameters µ, λ, N T . Among these, we find use of the zeroth, first, and second moments [1]; the second, third, and fourth moments [9,12]; the second, fourth, and sixth moments [2,3]; and the third, fourth, and sixth moments [4][5][6]. Table 1 lists the parameter estimates for these four different combinations of moments along with the use of the first, second, and third moments. For ease of presentation, the method of moments estimates for µ and λ in Table 1 are expressed in terms of values of an a appearing in the last column. The estimates in the last two rows of Table 1 do not make use of an auxiliary variable; to denote it these rows have an asterisk in the a column. We point out that each a value in Table 1 is at least one (note the (a − 1) denominators) because of Hölder's inequality [24] (Equation (3.2.10)) with a equal to one only in the pathological case where all the drop sizes are identical. Table 1 also includes two other types of gamma parameter estimates: (i) maximum likelihood (ML) estimates discussed, for example, in [11][12][13], and (ii) the pseudo maximum likelihood estimates by Ye and Chen in [14], who determined the maximum likelihood equations for a more encompassing generalized gamma distribution and then specialized these to the case of an "ordinary" gamma. The ML estimate for the shape parameter is obtained by numerically solving the equation (see e.g., [12]). The maximum likelihood estimate of the rate parameter is then given bŷ λ ML = (μ ML + 1)/D. Table 1 only lists estimates of the gamma parameters µ and λ because the large sample variance expressions forN T are especially complicated when using moment estimators, and partly because large sample maximum likelihood theory is developed only for estimates of density parameters (i.e., parameters appearing within just the density function f in Equation (1)).

Large Sample Behavior of the Estimates
For large sample sizes, each of the estimates in Table 1 is normal, unbiased, and has variance as given in Table 2. We used the delta method, described in detail in Appendix B, to establish these results for all but the maximum likelihood procedure, whose large sample results are given, for example, in [10]. It is not surprising that the delta method conclusions require large sample sizes as the delta method may be viewed as a multivariate version of the central limit theorem. For small samples, biases will occur and the variances will be larger than those stated in Table 2. To exemplify our large sample results more precisely, consider the method of moments estimate of, say, the shape parameter µ using moments 2, 3, and 4, from Table 1. Here, (1) The variance expressions for the Ye/Chen and ML estimates are all well-defined and positive for values of the shape parameter larger than −1. It is known, for example, that ψ (x) > 0 and xψ (x) − 1 > 0 for x > 0 from [25] (Equation (1.46)).
If E(μ MM234 ) denotes the expected value or mean of the estimated shape parameterμ MM234 and Var(μ MM234 ) denotes its variance, the claim is that E(μ MM234 ) µ, and where µ is the true value of the shape parameter in the RSD. The symbol indicates that the ratio of the left-and right-hand sides approaches one as the sample size increases. At the end of this section, more will be said about the sample size n needed for the large sample approximations in Table 2 to closely hold. For now, we mention that sample sizes as small as 2000 or more will suffice for several of the estimates. As stated in the Introduction, when choosing between several estimates, each of which are unbiased, we will generally prefer estimates with smaller variance (or smaller standard deviation, the square root of the variance). As the variance expressions in Table 2 are fairly complicated, it is best to compare the variabilities graphically.
We started by examining the estimates of the shape parameter µ. The estimate standard deviations which follow from Table 2 (by taking square roots) are all of the form 1/ √ n multiplied by a function of (just) the shape parameter µ. In Figure 1, we take the ordinates to be √ n multiplied by the standard deviation to compare the large sample variabilities of all of the estimates listed in Table 1. To clarify, Atmosphere 2020, 11, 333 5 of 11 the curve displayed in Figure 1 associated the MM123 estimate of the shape parameter, for example, has ordinates 2(µ + 2)(µ + 3)(µ + 6)/(µ + 1).
(1) The variance expressions for the Ye/Chen and ML estimates are all well-defined and positive for values of the shape parameter larger than −1. It is known, for example, that '( ) 0 x ψ > and '( ) 1 0 x x ψ − > for x > 0 from [25] (Equation (1.46)).
From Figure 1, we observe the moment estimates of the shape parameter incorporating the higherorder moments are least desirable as they have the highest variabilities, and the Ye/Chen and maximum likelihood estimates are most desirable as they have the lowest variabilities. The MM012 method has the smallest variability among the method of moments estimates and, in fact, compares favorably with the Ye/Chen and maximum likelihood methods.  From Figure 1, we observe the moment estimates of the shape parameter incorporating the higher-order moments are least desirable as they have the highest variabilities, and the Ye/Chen and maximum likelihood estimates are most desirable as they have the lowest variabilities. The MM012 method has the smallest variability among the method of moments estimates and, in fact, compares favorably with the Ye/Chen and maximum likelihood methods.
The Ye/Chen and maximum likelihood curves are practically identical and appear as a single curve in Figure 1. Consequently, only the single key entry of "ML/YC" is given in the Figure 1 legend. No attempt is made here to analytically quantify how close these two standard deviation curves are, but this is largely a consequence of being a near identity over much of the region x > −1.
Now we turn to the estimates of the rate parameter λ. It will again be easiest to compare estimates graphically. All estimate standard deviations that follow from Table 2 (by taking square roots) are of the form λ/ √ n multiplied by a function of (just) the shape parameter µ. In Figure 2, we take the ordinates to be √ n/λ multiplied by the standard deviation to compare the large sample variabilities of all of the estimates listed in Table 1. To clarify, the curve displayed in Figure 2 associated the MM123 estimate of the rate parameter, for example, has ordinates (2µ 2 + 23µ + 52)/(µ + 1)(µ + 2). μ μ μ μ + + + + As with the shape parameter, we see from Figure 2 that the moment estimates of the rate parameter incorporating the higher-order moments are least desirable as they have the largest variabilities, and the Ye/Chen and maximum likelihood estimates are the most desirable as they have the smallest variabilities. The MM012 method has the smallest variability among the method of moments estimates and compares favorably with the Ye/Chen and maximum likelihood methods. The large sample variabilities given in Table 2 are attained, in a limiting sense, as the sample size grows. A natural question, then, is how large the sample size must be for these variabilities to hold to good approximation. A careful analysis would involve looking at both of the estimates μ and λˆ across each of the seven different estimation procedures, with the behavior undoubtedly depending on the true parameter values of μ and λ. For our purposes here, we only very briefly discuss how empirical As with the shape parameter, we see from Figure 2 that the moment estimates of the rate parameter incorporating the higher-order moments are least desirable as they have the largest variabilities, and the Ye/Chen and maximum likelihood estimates are the most desirable as they have the smallest variabilities. The MM012 method has the smallest variability among the method of moments estimates and compares favorably with the Ye/Chen and maximum likelihood methods.
The large sample variabilities given in Table 2 are attained, in a limiting sense, as the sample size grows. A natural question, then, is how large the sample size must be for these variabilities to hold to good approximation. A careful analysis would involve looking at both of the estimatesμ andλ across each of the seven different estimation procedures, with the behavior undoubtedly depending on the true parameter values of µ and λ. For our purposes here, we only very briefly discuss how empirical values of the standard deviations through simulated observations from a gamma RSD match-up to values from Table 2 for large samples in two cases: µ = 2, λ = 5 and µ = 5, λ = 13 (essentially the values used in [22]).
For a sample size of 2000 and both pairs of the shape and rate parameter values just stated, we found the relative errors in the standard deviations for the MM012, MM123, ML, and Ye/Chenμ andλ estimates (but not the MM234, MM246, and MM346 estimates) to vary from about 1% to 3%. For the MM234 estimates with a sample size of 2000, we found relative errors for the standard deviations of the shape and rate parameters to each be about 1% in the case µ = 5, λ = 13, but were each about 10% in the case µ = 2, λ = 5. Sample sizes larger than 2000 are needed for the variabilities listed for MM246 and MM346 in Table 2 to be reasonably accurate.
Increasing the sample size to 10,000, the MM246 and MM346 estimates had relative errors for the standard deviations of the shape and rate parameters to each be about 5% in the case µ = 5, λ = 13 and to each be about 15% in the case µ = 2, λ = 5. Upon increasing the sample size to 100,000, these relative errors decreased to at most 0.8% in the case µ = 5, λ = 13 and to 4% to 5% in the case µ = 2, λ = 5 for both of these two method of moments procedures.
Atmosphere 2020, 11, 333 7 of 11 In practice, raindrop size samples are likely to be observed under stable conditions, corresponding to having fixed µ and λ parameter values, only over rather short time intervals, perhaps on the order of minutes. Multiple disdrometers would then be needed to collect sample sizes as large as, say 2000 drops. For more moderate sample sizes, the variabilities listed in Table 2 and displayed in Figures 1 and 2 should not, of course, be fully trusted. However, the ranking of estimates according to variability given in Table 2 and Figures 1 and 2 should persist to some degree as the sample size decreases and can be investigated through simulation. Such simulation for less than large sample sizes, however, is beyond the scope of this article.

Discussion and Conclusions
We used the delta method, a result known in the engineering and statistical literature, as a viable approach in deciding which of several gamma parameter estimates may be used when raindrop sizes can be accurately measured over the full spectrum of size. With the aid of the delta method, the estimates of the gamma shape and rate parameters listed in Table 1 were found, for large sample sizes, to be normal and unbiased with variances given in Table 2. Given several unbiased estimates, it is natural to prefer estimates with small variance.
If a method of moments procedure is to be used then to reduce estimate variance, we determined that gamma parameter estimates having the smallest order moments feasible should be preferred. We say 'feasible' as, for example, the lowest order moments might not always be reliable. Wind, for instance, can affect the reliability of the first moment. For some discussion related to this see, for example, [3,26,27].
We determined the maximum likelihood and pseudo maximum likelihood estimates of Table 1 have the smallest variances among the seven estimation procedures examined when estimating the gamma shape and rate parameters, but the MM012 method of moments procedure is very nearly as good in terms of variance.
When using instrumentation that does not accurately measure drop size over the full spectrum of drop size, then the variance results of Table 2 only hold approximately for large sample sizes. To compare gamma parameter estimates in this situation (or in the case of small sample size), samples of drop sizes from a gamma with known values of the shape and rate parameter may be randomly generated. Using these randomly generated samples, estimates may be compared in terms of their bias and variance. Such comparisons, for truncated and binned disdrometer data, were performed in [22] where we also suggested that the smallest order moments feasible be used when estimating gamma shape and rate parameters by method of moments.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix B Large Sample Calculations and the Delta Method
The large sample results for the maximum likelihood estimation given in Table 2 appear elsewhere in the literature (see, for example, [28] and [29](Equations (38) and (55))). The remaining large sample results in Table 2 all follow from the technique called the delta method [17][18][19]. The delta method details for the Ye and Chen estimates appear in [14]. As the remaining estimates in Table 1 are all method of moments estimates, and the delta method is perhaps not well-known, it seems appropriate to illustrate this technique for one of the method of moments procedures in Table 1. We did so for an MM234 estimate. We chose this particular combination of moments in large part because the calculations are a bit easier and the process is more transparent. With this example, it will be clear, in principle, how to compute the large sample properties for the other method of moments estimates in Table 1.
Before doing so, we mention what the delta method amounts to in the one-dimensional case. Readers may recognize this under the phrase "propagation of errors." Letting θ = E(D) denote the expected value or mean drop size, we have where g denotes the derivative of g, g is a function of a single variable. We conclude the normality of g(D) as a consequence of that of D (technically, use the fact that the right-hand side above is an affine transformation of D) and that g(D) has mean with the large sample normality of g(D) following from that of D by the central limit theorem.
The general delta method deals with a function of several arguments, each of which is approximately normal. The large sample mean and variance expressions of the delta method are generalizations of the example of propagation of errors above. The large sample normality depends on the fact that each of the sample moments m j are approximately normal, which is a consequence of the central limit theorem.
For the MM234 method of moments estimates, the general delta method concludes for large samples that the gamma parameter estimates given by (we drop the MM on the 234 subscript) are approximately normal with corresponding mean and with covariance matrix forμ 234 ,λ 234 given by the expression 1 n A A t where A is the Jacobian matrix and is the covariance matrix of the moments used The details are best accomplished using a computer algebra system that includes basic matrix computation (we used Maple 18 [30]), but we provide a few details showing how the rest of the computation is accomplished.
To carry out the above calculations, we require the population moments of raindrop size for our gamma distribution setting. Using the gamma density in Equation (2), these are given by for positive integer values of k. To determine, for example, the large sample mean ofμ 234 (the top entry in Equation (4)) we evaluate, from Table 1, ). Straightforward calculation shows that a evaluates to (µ + 4)/(µ + 3) giving (4 − 3a)/(a − 1) = µ. That is, E(μ 234 ) µ.