## Abstract

This paper proposes an estimation method of joint size and terminal velocity distribution on the basis of sampling data of precipitation particles containing multiple types. Assuming that the velocity follows the normal distribution and the size follows the gamma distribution, the method searches a locally maximum logarithmic likelihood within a realistic parameter range using the expectation–maximization algorithm. Several test populations were prepared with a realistic number of elements, and then the method was evaluated by retrieving the populations from their sample. The results showed that the original parameters were successfully estimated in most cases of the test population containing some of liquids, graupels, and rimed and unrimed aggregates. The original number of elements was also estimated with an adjustment of the number of elements in a manner such that each of their minority fractions exceeded a threshold. Applied to the two-dimensional disdrometer observation data, the method was helpful to discard frequently observed erroneous data with unrealistically large fall velocity.

## 1. Introduction

The precipitation particle size distribution (PSD) was already well known for the liquid phase (Marshall and Palmer 1948) and for the solid phase (Gunn and Marshall 1958). The number concentration rate of precipitation particles diameter *N*_{d} (m^{−3} mm^{−1}) follows a parametric gamma distribution on the diameter *D* (mm) as

where *N*_{0} (m^{−3} mm^{−1}), *μ*, and *λ* (mm^{−1}) denote the intercept parameter, the distribution shape parameter, and the slope parameter, respectively. A long-lasting problem is how observed samplings of precipitation particle size are fit to a parametric distribution. Any observation instruments available at present tended to miss small-size particles for their insufficient sensitivity and resolution, unless disdrometers like Meteorological Particle Sensor (Baumgardner et al. 2002) are combined to cover small sizes; this is hereinafter referred to as truncation problem (Handwerker and Straub 2011). Moreover large-size particles were rarely detected, which often brings us a kind of sampling problems. The moment method has been widely used for the fitting to size distribution (e.g., Vivekanandan et al. 2004; Smith and Kliche 2005; Yuter et al. 2006; Brandes et al. 2007; Kliche et al. 2008; Cao and Zhang 2009; Smith et al. 2009; Handwerker and Straub 2011), but a shortage of sampling length biases the estimate of parameters (Smith and Kliche 2005). The maximum likelihood and L-moment method can provide an unbiased estimate, but the truncation problem destabilizes the statistics (Kliche et al. 2008; Cao and Zhang 2009). Smith et al. (2009) showed that the above methods were difficult to fit limited sampling data to a gamma distribution of a small shape parameter. Recently, Yano et al. (2018) relieved this problem in the fitting to a gamma distribution solely on the basis of bulk quantities of PSD by applying the maximum entropy principle (Yano et al. 2016). Another problem comes up when fitting the mixture of multiple particle types in a single sampling (Brandes et al. 2007). For example, sampling data possibly contain graupels and aggregates in a time interval (Yuter et al. 2006) that should be long enough to avoid the sampling problem (Smith and Kliche 2005).

On the other hand, the particle terminal velocity *V* (m s^{−1}) is assumed to be proportional to the power of the particle’s diameter as

where *a*_{υ} and *b*_{υ} depend on particle phase, shape, and degree of riming (Atlas et al. 1973; Locatelli and Hobbs 1974). It is remarked that another functional form may be possible for liquid particles (Atlas et al. 1973) but we chose the simpler form as Eq. (2) to avoid a superfluous parameter throughout the paper. Many publications have been devoted to find an optimal fitting curve for particle diameter and velocity, on the basis of observed sampling data with optical or video disdrometers such as Particle Size Velocity (PARSIVEL; Löffler-Mang and Joss 2000; Battaglia et al. 2010; Angulo-Martínez et al. 2018), Laser Precipitation Monitor (LPM; Angulo-Martínez et al. 2018), and two-dimensional video disdrometer (2DVD; Kruger and Krajewski 2002). Recently, Molthan and Colle (2012) compared on-site observation with model outputs in every separated cell on the diameter–velocity quarter plane, instead of the direct comparison of velocity–diameter relationship or probability density function (PDF) on the plane. Bernauer et al. (2016) classified snow particles observed by 2DVD into three ranges of riming degree and four crystal types. The above fitting to a power law like Eq. (2) is, however, very difficult in a case of mixed particle types, because the least squares method conventionally used assumed a single type of precipitation particles. Enough of a short time interval, like 1 min (Bernauer et al. 2016), could partly, not perfectly, avoid the mixing, around the period when two independent distributions are clearly transited, by not sampling data from both the distributions, but would cause the sampling problem. A fitting to mass flux, a production of size–velocity and size–mass relations, partially resolved this problem (Ishizaka et al. 2013), but it is still difficult in a case of a clear bimodal distribution due to coexistence of liquid and solid phases (Brandes et al. 2007; Yuter et al. 2006) or coexistence of graupels and aggregates. Another problem in the fitting to size–velocity relation is a sampling error in small-size particles with an unrealistically high speed. This was often found in the output of single-sensor type disdrometers such as PARSIVEL and LPM (Angulo-Martínez et al. 2018), when particles passed through the edge of the sensitivity frame (Minda et al. 2016). This was also found in a double line-sensor type disdrometer of 2DVD for a failure of matching (Huang et al. 2010; Bernauer et al. 2015). One inevitably excluded the part of sampling data out of the velocity range between 0.5 and 5 m s^{−1} (Huang et al. 2010) and excluded data of small size less than 1 mm (Yuter et al. 2006; Bernauer et al. 2015, 2016), although the method filtering drops outside the ±40% of the theoretical fall velocity may be also applicable if all particles are guaranteed to be liquid (Kruger and Krajewski 2002). This treatment deteriorated the truncation and sampling problems in the fitting to a size distribution, however.

This study aims to relieve the above problems of the mixing and the erroneous data often observed by disdrometers by fitting diameter–velocity data including multiple precipitation particle types and phases to mixed joint PDF on the diameter–velocity quarter plane; We conveniently utilized the gamma and normal distributions and the velocity–diameter relationship to evaluate the PDF in the comparison with a conventional method and well-known velocity–diameter relationships previously proposed. The fitting problem can be then reduced to a numerical search of an optimal set of parameters of the mixed joint PDF from the diameter–velocity data. We here develop the method with an aid of the expectation–maximization (EM) algorithm (Dempster et al. 1977). A performance of the new method is evaluated by a PDF retrieval from test data randomly sampled from a given population PDF. We next emphasize the strength of the new method to discern the coexistence of precipitation particle types and the sampling error in a subset of 2DVD data. This paper is organized as follows. Section 2 describes the observed 2DVD data, and section 3 elucidates the proposed method including a nutshell of the EM algorithm. We design experiments to check performance of the new method in section 4 and detect a bundle of erroneous data by applying the method to sample datasets obtained by the 2DVD in section 5. Section 6 discusses the limitation of the new method and section 7 concludes this paper.

## 2. Observation data

We used the precipitation diameter–velocity data obtained by a 2DVD installed at a site of the Institute of Low Temperature Science, Hokkaido University, Sapporo, Japan, located at 43.083°N, 141.339°E (Nagumo and Fujiyoshi 2015; Campbell et al. 2018). The 2DVD (Kruger and Krajewski 2002), widely used in the scientific community (Huang et al. 2010; Bernauer et al. 2015, 2016), measures precipitation particle diameter and velocity by matching two images captured by double line sensors installed at an upper and a lower position of the instrument. Image processing algorithm originally implemented in the 2DVD sensor matched the two particle images. As a preprocessing, we discarded erroneous sampling particles with their velocity larger than 15 m s^{−1}, much larger than the terminal velocity of rain droplets (Atlas et al. 1973). We used surface air temperatures and precipitation records by Sapporo Meteorological Observatory of the Japan Meteorological Agency (JMA) located at 43.0605°N, 141.329°E, 2.6 km away from the 2DVD site (JMA 2019a,b).

## 3. Method

### a. Mixed joint PDF

The PDF of particle diameter was here assumed to be the gamma distribution, a normalized version of Eq. (1):

where Γ is the gamma function. Here *P*(*D*|*μ*, *λ*) indicates conditional probability density of *D* under parameters *μ* and *λ*. We assume that the particle’s terminal velocity follows the normal distribution with its mean being Eq. (2):

where *σ*^{2} is variance of velocity. This assumption was supported by theoretical experiments (Schmitt et al. 2019) and observations (Sasyo and Matsuo 1980). In multiplying Eqs. (3) and (4), a joint PDF of particle diameter and velocity is

Linearly combining independent joint PDFs, a mixed joint PDF with *K* elements was written as

where *ω*_{k} is a mixing fraction of the *k*th PDF element and should satisfy $\u2211k=1K\omega k=1$. Now we introduced the vector ** θ** whose components are conditional parameters as

Many possible PDFs can be constructed by sweeping a set of parameters ** θ** within a realistic range. Given a sampling dataset, we may search an optimal parameter set

**to estimate the PDF to be fit to the data. The distribution of particle number concentration rate within diameter–velocity quarter plane,**

*θ**N*

_{vd}(m

^{−3}s m

^{−1}mm

^{−1}), more suitable unit for the traditional definition in meteorology, was easily obtained by multiplying the number concentration of total precipitation particles

*N*

_{tot}(m

^{−3}):

Note that the PDF [Eq. (6)] and the number concentration rate [Eq. (7)] explained diameter–velocity data by a disdrometer directly observing particle number concentration (e.g., Muramoto and Shiina 1989).

The joint PDF following diameter–velocity data by disdrometers observing precipitation particles passing a finite area such as 2DVD, PARSIVEL, and LPM is described as

From this equation, the shape parameter *μ* is shifted by *b*_{υ} with a diameter–velocity relationship form of Eq. (2) (Ignaccolo and De Michele 2014; Adirosi et al. 2016) and *P*_{S}(*V*, *D*) extends to the tendency of higher velocity due to a factor of $V/a\upsilon Db\upsilon $. The mean fall velocity for the joint PDF of Eq. (8), *V*_{S}, is also different from Eq. (2):

The mixed joint PDF based on Eq. (8) is

The distribution of particle number concentration rate is

where *A*, Δ*t*, and *L* are, respectively, an observing area (m^{2}), the sampling time (s), and the number of total particles observed by a disdrometer. See appendix A for derivations of Eqs. (8) and (11).

Equations (5) and (8) should be properly assumed as a joint PDF form that diameter–velocity data by a disdrometer should follow. Because we could not formulate the EM algorithm for parameters *a*_{υ} and *b*_{υ} with this joint PDF form of Eq. (8) (cf. section 3b), however, the PDF form of Eq. (5) was applied to the observed sampling data by 2DVD. The problem on using Eq. (5) for 2DVD data will be discussed in section 6.

### b. EM algorithm

The EM algorithm (Dempster et al. 1977) was used to search an optimal parameter set of Eq. (6). First, we posed a kernel PDF from the *l*th particle for the *k*th PDF element as $P\u2061(Vl,Dl|a\upsilon k,b\upsilon k,\sigma k2,\u2009\mu k,\u2009\lambda k)$. The logarithmic (log) likelihood was then constructed by summing the logarithm of this kernel PDFs for all of *L* particles in the sample:

The EM algorithm finds a parameter set of ** θ** that satisfies

to maximize the log likelihood [Eq. (12)] under the constraint of $\u2211k=1K\omega k=1$. Equation (13) is reduced to the iteration relations below:

where

is the digamma function, the responsibility fraction *γ*_{lk} ranging over 0–1 as

represents how much the *l*th particle explains the *k*th PDF element, and $D\xafk$ and $ln\u2061(D)\xafk$, respectively, represent the mean diameter [$=\u2061(\u2211l=1L\gamma lkDl)/\u2061(\u2211l=1L\gamma lk)$] and the mean log diameter {$=[\u2211l=1L\gamma lk\u2061ln\u2061(Dl)]/\u2061(\u2211l=1L\gamma lk)$} corresponding to the *k*th PDF element. We numerically solved Eq. (15) by the secant method (Householder 1970), and solved Eq. (17) by iteration (Kliche et al. 2008):

Note that Eqs. (17) and (18) with *K* = 1 are the completely same form as the equations of the conventional maximum-likelihood method in estimating particle size distribution (Kliche et al. 2008) so that the EM algorithm with *K* = 1 estimates *μ* and *λ* as comparably as the maximum-likelihood method.

Next, the EM algorithm estimates an optimal parameter set including responsibility fraction [Eq. (20)] in the iteration procedure. Initially, we give the number of PDF elements *K* and parameters for the *k*th PDF element as

where

with overbar and prime denoting the average and the variance for all the particles. See appendix B for details of this initial parameter determination. We call providing the responsibility fraction for the parameter set with Eq. (20) the E step and updating the parameter set with a sequential calculation order of *b*_{υ} with Eq. (15) substituting Eq. (14) for *a*_{υ}, *a*_{υ} with Eq. (14), *σ*^{2} with Eq. (16), *μ* with Eq. (17), *λ* with Eq. (18), and *ω* with Eq. (19) the M step. The E step and M step are iterated if the log likelihood increases by more than 10^{−2} relative to the value before the steps. Otherwise the iteration stops. No more iterations than 5000 avoids computation for divergent cases. It is noted that a convergence solution by this iteration procedure is solely a locally most likelihood nearest from the initial and is not always the most likelihood *globally*, and therefore the EM algorithm is not referred to as the maximum-likelihood method.

### c. Estimation of the number of PDF elements

A remaining parameter that should be estimated is the number of PDF elements, *K*, even though the EM algorithm itself can be performed with given *K*. We now introduce the Akaike information criterion (AIC; Akaike 1974)

and Bayesian information criterion (BIC; Schwarz 1978)

An appropriate number of PDF elements can be searched by maximizing AIC or BIC by sweeping *K* from 1 to 4. We also applied minority rejection method to find the appropriate number of elements by an adjustment such that the smallest mixing fraction in the elements is larger than a threshold following 0.1 K. See appendix C for the rationality of this choice.

## 4. Performance test

### a. Design

Performance of the proposed method is examined in light of (i) computational stability of the EM algorithm, (ii) the accuracy on the parameter forming the PDF elements, (iii) the accuracy on the estimate of the optimal PDF-element number, and (iv) comparison with the least squares method conventionally used. All points can be checked by experiments, similarly designed, to estimate a given PDF from random sampling data from the PDF, just like a Monte Carlo simulation (Smith and Kliche 2005; Kliche et al. 2008; Cao and Zhang 2009; Smith et al. 2009; Handwerker and Straub 2011).

For points i–iii, we prepared nine population distributions containing PDF elements including liquid precipitation particle L with (*a*_{υL}, *b*_{υL}) = (3.78, 0.67) (Atlas et al. 1973), lump graupel G with (*a*_{υG}, *b*_{υG}) = (1.3, 0.66) (Locatelli and Hobbs 1974), densely rimed aggregates R with (*a*_{υR}, *b*_{υR}) = (1.1, 0.15) (Ishizaka 1995), or side-plane unrimed aggregates U (*a*_{υU}, *b*_{υU}) = (0.82, 0.12) (Locatelli and Hobbs 1974). The mixed joint PDFs prepared here are mixture of liquid particle, graupel, and rimed aggregates (LGR), mixture of graupel and rimed aggregates (GR), and mixture of rimed and unrimed aggregates (RU), with the same mixing fraction for all elements (i.e., *ω* = 1/3 for LGR and *ω* = 1/2 for GR and RU). Note that these combinations of particle types were subjectively chosen to represent some of possible observed examples. For point iv, a single PDF element with liquids (L), graupels (G), rimed aggregates (R), or unrimed aggregates (U) was prepared. We set the shape parameter as *μ* = 1.5 for all of the populations. Rather than the parameters (*a*_{υ}, *b*_{υ}, *ω*, *μ*), the PDF elements depend on data variance *σ*^{2} and slope *λ*. We set three different combinations of these two parameters to check a variety of mixed joint PDFs as Table 1. Three kinds of PDF element combined with three parameter sets of *σ*^{2}, *μ*, and *λ* made nine mixed joint PDFs (Fig. 1), and each experiment is nominated as the prefix of LGR, GR, or RU combined with the suffix of 0, s, or v. For example, the experiment LGRv means a given PDF of elements LGR with a parameter set of v. Next, we randomly sampled diameter–velocity data from the population distribution by von Neumann’s method within ranges of 0 ≤ *D* ≤ 30 mm and 0 ≤ *V* ≤ 10 m s^{−1}. The sensitivity to the sample size was checked by taking 100, 1000, 2000, 3000, 6000, and 10 000 particles. Last, the EM algorithm was applied to the sampled data and the obtained PDF was compared with the original. Each estimation was repeated 1000 times with different sampling data. To keep the same element order between estimate and population, the estimated elements are matched to the elements in the population by selecting the combination of (*j*, *k*) with the nearest neighbor measured by normalized Euclidian distance of

Here *j* represents the *j*th PDF elements in the estimation (*j =* 1, 2, 3 in LGR and *j* = 1, 2 in GR and RU), and *k* represents PDF element L, G, R, or U.

Point i is to check computational stability measured by the convergence rate of EM algorithm, in which we counted how many trials were numerically converged in 1000 trials with a different sample dataset. This check should be required because the EM algorithm, searching a local optimal solution nearest from an initial value, not a *global* optimal solution, is not always converged. We regarded that the trial was successful when the computation was converged for every number of PDF elements given from 1 to 4. Point ii was to check the accuracy of the parameters forming the PDF elements. This was evaluated as the root-mean-square error (RMSE) of parameters, *a*_{υ}, *b*_{υ}, *σ*^{2}, *μ*, *λ*, and *ω*, of the estimated PDF elements from the reference of the population PDF given (Fig. 1). The relative error of log likelihood was also used as the measure. Point iii was to check the accuracy of the estimate of the number of PDF elements *K* from 1 to 4. When the number of PDF elements in the estimation was matched with the given number, we regarded it as hit. The hitting rate is defined as the hit number divided by the total trials that are converged. Point iv was to compare our method that could generally estimate a mixed joint PDF with the conventional method that estimated only a single velocity–diameter relationship. We then used the single-element PDF as the population with three different combination patterns of variance and slope parameters.

### b. Results

The computational stability was first checked by the convergence rate of estimation. Whatever we gave the PDF element number, the EM algorithm was converged by mostly 100% of 1000 trials if the sample size was greater than 1000 (Fig. 2a). In a small sample size of 100, the EM algorithm sometimes did not find a convergence solution with greater PDF elements given (Fig. 2b). Therefore, it is not recommended to give the number of PDF elements greater than 2 when the sample size is less than 100. Hereinafter we show the results with the sample data size at 2000.

Second, the accuracy of the estimation of PDF elements in comparison with the given population was checked by RMSE of each parameter (Fig. 3), in the estimation with the truth number of PDF elements *K* given. Greater the sample size was, more accurate the estimate was in most cases; the RMSE of any parameters plummeted with *L* exceeding 2000 (not shown). Parameters *a*_{υ} and *b*_{υ}, directly determined a diameter–velocity relationship, were successfully estimated, especially in the retrieval of graupel and rimed aggregates with small variance (GR0 and GRs; Figs. 3a,b). The population containing rimed and unrimed aggregates with larger variance (RUv) was estimated with relatively larger RMSE. The variance parameter, *σ*^{2}, of PDF elements was also well retrieved with its RMSE almost less than a quarter of variance parameter in the PDF elements in the given population except for liquids for LGRv (Fig. 3c). The parameter *μ* was estimated with its RMSE mostly less than 0.5 (Fig. 3d). The slope parameter *λ* was well retrieved in all the cases with its RMSE less than 10% of the slope parameter in the PDF elements in the given population (Fig. 3e). This estimate from smaller sample size than 2000 was much worse in RU0, RUs, and RUv (not shown). The RMSE of mixing fraction *ω* ranged over 0.01–0.05 except for the estimation for RUv population (Fig. 3f).

Figure 4 shows the PDF of the relative error of log likelihood based on 1000 trials, in which positive and negative values mean overfitting and mismatching, respectively. None of the trials fell into a wrong optimum except for RU0 and RUs. In the cases of RU0 and RUs, the distribution of relative error shows bimodality separated at zero (Figs. 4c,f). If the relative error of log likelihood below −0.2% was regarded as a wrong optimum, about 10% of the trials failed to estimate the true PDF. Because the error majority ranged from 0% to 0.5%, the PDF estimate was slightly biased to the overfitting side. It is worthwhile noting that this width of range was probably attributed to the uncertainty due to a sampling randomness indwelling in the sampled data, which was small enough here (Fig. 4). Therefore, the overfitting was the primary cause for the RMSE in estimated parameters (Fig. 3), but the mismatching was the secondary cause only in the cases of RU0 and RUs.

Third, we evaluated the number of PDF element by the hitting rate (Fig. 5). The result showed that the minority rejection method based on smallest mixing fraction provided the best performance in all the population patterns given (Fig. 5c). The hitting rate was asymptotic to 100% for greater sample size. The sample size of 1000 or more made the hitting rate of 100% for the populations of LGR0, LGRs, LGRv, GR0, GRs, and GRv based on the method maximizing BIC (Fig. 5b). The method maximizing AIC provided the worst estimates (Fig. 5a), very sensitive to the sample size; the hitting rate is too low for the population with the large variance parameter even the sample size is large.

Fourth, our method for the estimate from the population containing the single PDF element was compared with the conventional method to find an optimal velocity–diameter relationship by the least squares method: the comparison was done for two parameters *a*_{υ} and *b*_{υ} (Fig. 6). The populations regardless of parameters *σ*^{2} and *λ* were retrieved with the comparable RMSE between the methods. The difference between methods was insensitive to the sample size (not shown). Hence, our method can be also applicable to a single-element population in addition to a multiple-elements population. Comparing the result of single PDF element cases (Fig. 6) and multiple PDF elements cases (Fig. 3), the estimations for the joint PDF population of GR0 and GRs had as a good performance as two individual estimations for the single PDF population G0 or R0 (Figs. 3 and 6). In contrast, the estimation for the RUs and RUv populations had a worse performance than two estimations for Rs, Rv, Us, or Uv population.

## 5. Results

In addition to the performance test above, we here demonstrated how the proposed method fit to diameter–velocity data observed by the 2DVD at Sapporo (section 2), with two cases not to be considered in the performance test. The first case included erroneous and normal data during 1400–1430 Japan standard time (JST) 22 January 2017, when the JMA reported the snow particles fell in the air temperature of −4.9°C (JMA 2019a). The second case included liquid and wet snow particles during 0100–0130 JST 6 December 2016, when the JMA reported that the rain turned to sleet and the air temperature changed from 3° to 2.2°C (JMA 2019b).

The observed data by the 2DVD during 1400–1430 JST 22 January 2017, included erroneous data with unrealistically large velocity probably due to the matching problem (Fig. 7a). With two PDF elements given, our method estimated one element (red in Fig. 7a) with its center at *V* = 2.61*D*^{0.54} and the variance of 6.43 m^{2} s^{−2}, which was so unrealistic. Moreover, this erroneous PDF element underestimated probability around 1.5 m s^{−1} (red in Fig. 7c) due to the assumption that the velocity followed the normal distribution in this study. In contrast, the erroneous PDF well represented diameter distribution (red in Fig. 7b), because the matching error affected the velocity only. The estimation of the other PDF element (blue in Fig. 7a) was successfully fit to *V* = 0.83*D*^{0.21}, near to the optimal curve for graupel-like snow of hexagonal type, *V* = 0.86*D*^{0.25}, or aggregates of densely rimed radiating assemblages of dendrites, *V* = 0.79*D*^{0.27}, in Locatelli and Hobbs (1974) (Fig. 7a). This PDF element well represented the observed distributions of both the diameter and velocity (blue in Figs. 7b,c). However, because the estimated velocity–diameter relationship roughly diagnosed the characteristics of precipitation particles (Bernauer et al. 2016), the estimated velocity–diameter relationship *V* = 0.83*D*^{0.21} did not always guarantee the precipitation particles containing graupel-like snow and aggregates of densely rimed radiating assemblages. While the fitting result was obtained with *K* = 2 subjectively given, the minority rejection method estimated the number of PDF elements at *K* = 1 because the minimum mixing fraction with *K* = 2, *ω* = 0.18 (Fig. 7a), was less than the threshold of 0.2 determined in advance (section 3c).

Interesting is a different case for mixing rain droplets and wet snow particles observed with the 2DVD at Sapporo during 0100–0130 JST 6 December 2016 (Fig. 8). Here the minority rejection method estimated two PDF elements. If we subjectively fixed the number of elements at 3, we estimated the first element with its center at *V* = 2.33*D*^{0.03} and the variance of *σ*^{2} = 0.13 (green in Fig. 8), the second element with its center at *V* = 2.9*D*^{0.53} and the variance of *σ*^{2} = 0.12 (blue in Fig. 8), and the third element with its center at *V* = 2.83*D*^{0.62} and the variance of *σ*^{2} = 3.46 (red in Fig. 8). The first element resembles a distribution by wet snow particles reported by Yuter et al. (2006), and the second element was similar to a diameter–velocity curve of rain droplets (*V* = 3.78*D*^{0.67} in Atlas et al. 1973) with small variance and good fitness of velocity distribution (blue in Fig. 8c). The third element probably contained erroneous data because of its large variance and inappropriate fitting to velocity distribution (red in Fig. 8c), however. The minority rejection method estimated the number of elements at *K* = 2, probably because the number of erroneous data was smaller as shown in Fig. 7.

## 6. Discussion

We showed that the proposed method successfully reproduced the original PDF from the sampling data of its size larger than 2000 under a limitation of possible PDF element number considered (Fig. 3). However, as in the application to observation data including error (section 5), we must take a treatment for erroneous data as the preprocessing because an arithmetic underflow often occurred in representing their probability. For example, the probability of particles at (*D*, *V*) = (2, 15) under a condition of (*a*_{υ}, *b*_{υ}, *σ*^{2}, *μ*, *λ*) = (0.82, 0.12, 0.09, 1.5, 0.3) in Eq. (5) approximates 10^{−481} that causes a double-precision underflow. This might cause the zero division error in the calculation of the responsibility fraction [Eq. (20)] with uniform mixing fractions among the PDF elements. This underflow problem can be relieved by setting the initial variance parameter larger than *σ*^{2} = 0.16 and by discarding particles with their velocity faster than 15 m s^{−1}. If one successfully avoided this problem by the preprocessing like above, the EM algorithm would make light of erroneous extreme data and search an optimal set of PDF elements that follow data majority.

Moreover, we still suffer from the truncation effect attributed to not only the small size particles but also slow velocity ones. To relieve the problem, we now consider the truncated version of mixed joint PDF. The truncated version of particle size distribution is

where *D*_{min} and *D*_{max} are, respectively, the minimum and maximum diameters that a disdrometer is capable to measure, Π is the rectangle function (1 between *D*_{min} and *D*_{max} and 0 elsewhere), and *Q* is the cumulative density function (Mallet and Barthes 2009; Johnson et al. 2014). Similar to Eq. (27), the truncated version of terminal velocity distribution is

where *V*_{min} and *V*_{max} are, respectively, the minimum and maximum velocity that a disdrometer is capable to measure. Therefore, the truncated version of mixed joint PDF is derived by linearly combing the joint PDF of *P*_{T}(*V*|*D*, *a*_{υ}, *b*_{υ}, *σ*^{2}) × *P*_{T}(*D*|*μ*, *λ*). Modifying the EM algorithm using this truncated version of mixed joint PDF, the proposed method might relive the truncation problem, but this is beyond the scope of this paper.

We next discuss the time interval for data sampling in the observation. The time interval that the total number of precipitation particles passing through an area of *A* (m^{2}) reaches *L* is computed in appendix D as

where *R* (mg m^{−2} s^{−1}) is precipitation intensity, *a*_{m} and *b*_{m} are parameters of the mass–diameter relationship: $M=amDbm$ (mg). Note that (*a*_{m}, *b*_{m}) for liquids, graupels and densely rimed and unrimed aggregates are (0.52, 3), (0.078, 2.8), (0.094, 1.9), and (0.04, 1.4), respectively (Atlas et al. 1973; Locatelli and Hobbs 1974; Ishizaka 1995). Figure 9 shows the time interval to obtain 2000 particles per a unit area of ~0.01 m^{2} composed by the typical 2DVD (Kruger and Krajewski 2002). Since the sample size should be greater than 2000 for computational stability of the algorithm (Fig. 2) and for accuracy of the parameter estimation (Figs. 3 and 5), in 1 mm h^{−1} precipitation intensity, the time interval is roughly 5 min for LGRs, GRs, and RUs populations, 30 min for LGR0, LGRv, GR0, GRv, RU0, and RUv. Moreover, with an intake of PARSIVEL and LPM of almost one-half of 2DVD, the time interval is doubled (Angulo-Martínez et al. 2018).

We similarly compared our method with the conventional maximum-likelihood method by Kliche et al. (2008) in light of the shape *μ* and slope *λ* parameters from the single population of liquid precipitation particles following the gamma distribution. They reported that the RMSE of shape parameter was 0.13 with respect to the population of *μ* = 2 for the sample size of 1000. This result was consistent with our result of the estimated PDF element corresponding to the liquid particles for the LGR0, LGRs, and LGRv populations (Fig. 3d). As for the slope parameter, Kliche et al. (2008) computed a median value of estimated *λ* normalized with its population being ~1.0 with a sample size of 1000. In this study, the normalized *λ* of the estimated PDF element corresponding to the liquid particles was ~1.02 for LGR0, ~1.03 for LGRs, and ~1.05 for LGRv (not shown). Hence, the performance of our proposed method with *K* = 3 was slightly worse than the maximum-likelihood method, completely same as the proposed method with *K* = 1 (section 3b). Moreover, Kliche et al. (2008) computed a median value of normalized *λ* estimated with the moment method being ~1.05 with a sample size of 1000, which is mostly the same performance as that of the proposed method with *K* = 3. Therefore, the proposed method can estimate the particle size distribution discriminating PDF elements as accurately as the moment method with a single PDF.

Last, we discuss the problem that the diameter–velocity data by disdrometers such as 2DVD, PARSIVEL, and LPM originally follows the different form of joint PDF [Eq. (8)] from the joint PDF assumed in this study [Eq. (5)]: a shifted shape parameter as much as *b*_{υ} [Eq. (8)] and an offset of $\sigma 2/\u2061(a\upsilon Db\upsilon )$ in the mean terminal velocity [Eq. (9)]. The difference in terminal velocity is larger for the smaller particle diameter (Fig. 10), so that parameter *a*_{υ} and *b*_{υ} should be, respectively, overestimated and underestimated under an assumption of the joint PDF form by Eq. (5) comparing with the estimation assuming the joint PDF form by Eq. (8). However, the difference in the mean velocity is mostly less than 0.1 m s^{−1} with a larger diameter than 1 mm in the parameter range used in this study (Fig. 10). Therefore, the estimation of *a*_{υ} and *b*_{υ} might not be largely affected regardless of the joint PDF forms expressed with Eqs. (5) and (8) if the sampling data consists of large particles. If so, the shape parameter estimated by the proposed method can be simply converted to the shape parameter estimated with an assumption of the joint PDF form of Eq. (8) by an offset of *b*_{υ}.

## 7. Conclusions

A new method with the EM algorithm was developed to estimate an optimal mixed joint PDF from precipitation particles’ diameter–velocity data. The method was verified by the test to retrieve the given PDF from a randomly sampling data based on a population that contains some of PDF elements: liquid particles, graupels, densely rimed aggregates, and unrimed aggregates (Fig. 1). The performance test suggested that the sample size should exceed 2000 for the computational stability (Fig. 2) and for accuracy of parameter estimation (Fig. 3). In the estimate of the number of PDF elements, minority rejection method made a better performance than maximizing AIC or BIC (Fig. 5). Our estimate made a comparable performance to the conventional least squares method (Fig. 6) and the maximum-likelihood method (section 3b) in the estimate of a single-element population. Applying this method to the observation data with 2DVD, our method was able to detect erroneous data (Fig. 7) and to distinguish liquid-phase from solid-phase precipitation particles (Fig. 8).

## Acknowledgments

The authors thank Dr. Naoto Nakano for giving us insightful comments on a mathematical limitation on the EM algorithm. The 2DVD data used in this study were provided by Prof. Yasushi Fujiyoshi. Observation data of surface air temperature and precipitation were obtained online from the JMA. Author Katsuyama is fully supported by JSPS KAKENHI Grant 18J12196. Author Inatsu is partly supported by JSPS KAKENHI Grants (18H03819, and 19H00963); by the Social Implementation Program for Climate Change Adaptation Technology funded by the Ministry of Education, Culture, Sports, Science, and Technology; by the Environment Research and Technology Development Fund 2-1905 of the Environmental Restoration and Conservation Agency of Japan; and by Research Field of Hokkaido Weather Forecast and Technology Development (Endowed by Hokkaido Weather Technology Center Co., Ltd.). A software source code of the developed fitting algorithm will be available online (https://humet.sci.hokudai.ac.jp/~meteo/product-e.html).

### APPENDIX A

#### Derivation of Eqs. (8) and (11)

Since the disdrometers such as 2DVD, PARSIVEL, and LPM observes precipitation particles passing a finite area, the observed diameter–velocity data should follow the joint PDF as follows (Ignaccolo and De Michele 2014):

Here, *P*(*V*, *D*) is the joint PDF expressed with Eq. (5). The total number of precipitation particles is represented as

### APPENDIX B

#### Initial Parameters for the EM Algorithm

The initial parameters for the EM algorithm should be reasonably determined for an accurate parameter estimation. First, we determine two velocity–diameter curves approximately passing maximum and minimum bounds of velocity on the diameter–velocity quarter plane: one is passing a point of $\u2061(D\xaf,\u2009max1\u2264l\u2264LVl)$ on the diameter–velocity quarter plane with *b*_{υ} = 1, and another is constant at $V=V\xaf/4$. These conditions are satisfied with two parameters sets of $\u2061(a\upsilon ,b\upsilon )=(a^\upsilon ,1)$ and $\u2061(a\upsilon ,b\upsilon )=(\u0103\upsilon ,0)$, where $a^\upsilon $ and $\u0103\upsilon $ are given by Eq. (23). Then, *K* sets of (*a*_{υ}, *b*_{υ}) are equally spaced out between $a^\upsilon $ and $\u0103\upsilon $ and between 1 and 0 as per Eq. (22). Parameters *σ*^{2}, *μ*, and *λ* are simply determined by the mean and variance of diameter–velocity data and are uniformly given to all the PDF elements. The mixing fraction *ω* is also uniformly given.

### APPENDIX C

#### Threshold for the Minority Rejection Method

The threshold value for minority rejection method is given by *cK*, where *c* is an arbitrary positive constant and *K* is the number of PDF elements. The threshold value by this formulation effectively reduces a possible range of *K* because the *K* estimation with the minority rejection method never exceeds an expectation value: *cK* ≤ 1/*K*. Hence

This limitation is probably helpful to avoid an overfitting problem.

Figure C1 shows a hitting rate of the number of PDF elements averaging all cases of sample size (*L* = 100, 1000, 2000, 3000, 6000, 10 000) with the minority rejection method as a function of *c*. The estimation were achieved by mostly 100% for GR0, GRs, and GRv populations within *c* ranging over 0.1–0.2. The number of PDF elements in RU0, RUs, and RUv populations were also successfully estimated with a relatively narrower range of *c*. On the other hand, for LGR0, LGRs, and LGRv populations, the PDF element number were estimated only within a narrow range of *c* around 0.05. The results suggested that the hitting rate becomes more sensitive to *c* in a case for estimating a population constructed from more number of PDF elements. Eventually, an averaged hitting rate over all population tests was approximately 90% at *c* = 0.1 (gray thick line in Fig. C1), which limits a range of *K* to 1–3 [Eq. (C1)]. We hence used *c* = 0.1 for the minority rejection method in the text.

### APPENDIX D

#### Derivation of Eq. (29)

The relationship between mass and diameter of precipitation particles is empirically described by the power law as

where *M* (mg) is the mass of a particle and *a*_{m} and *b*_{m} are empirical coefficients (Locatelli and Hobbs 1974). Multiplying *M* [Eq. (D1)] and *V* [Eq. (2)], a mass flux *F* (mg m s^{−1}) is obtained:

Assuming that a PSD follows the gamma form represented by Eq. (1), precipitation intensity *R* (mg m^{−2} s^{−1}) is obtained by integrating *N*_{d}*F* over *D*:

Similarly, particle number flux $L\u02dc$ (m^{−2} s^{−1}) is obtained by integrating *N*_{d}*V* over *D*:

For *K* elements, a linear combination of independently estimated particle number flux $L\u02dck$ following Eq. (D5) brings

The total particle number passing through an area *A* (m^{2}) within a time-interval Δ*t*, $L=A\u222b0\Delta tL\u02dc\u2009d\tau $, provides Eq. (29).

## REFERENCES

*Adv. Water Resour.*

*IEEE Trans. Autom. Control*

^{2}and Thies LPM optical disdrometers

*Hydrol. Earth Syst. Sci.*

*Rev. Geophys.*

*J. Atmos. Oceanic Technol.*

*Atmos. Meas. Tech.*

*Atmos. Res.*

*J. Appl. Meteor. Climatol.*

*Mon. Wea. Rev.*

*J. Appl. Meteor. Climatol.*

*J. Roy. Stat. Soc.*

*J. Meteor.*

*J. Atmos. Oceanic Technol.*

*The Numerical Treatment of a Single Nonlinear Equation.*McGraw-Hill, 216 pp

*J. Atmos. Oceanic Technol.*

*J. Appl. Meteor. Climatol.*

*Seppyo*

*J. Meteor. Soc. Japan*

*Quart. J. Roy. Meteor. Soc.*

*J. Appl. Meteor. Climatol.*

*J. Atmos. Oceanic Technol.*

*J. Geophys. Res.*

*J. Atmos. Oceanic Technol.*

*J. Atmos. Oceanic Technol.*

*J. Meteor.*

*IEEJ Trans. Electr. Electron. Eng.*

*Mon. Wea. Rev.*

*J. Inst. Electron. Inf. Commun. Eng.*

*Mon. Wea. Rev.*

*Pap. Meteor. Geophys.*

*J. Appl. Meteor. Climatol.*

*Ann. Stat.*

*J. Appl. Meteor.*

*J. Appl. Meteor. Climatol.*

*J. Appl. Meteor.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Appl. Meteor. Climatol.*

## Footnotes

^{a}

Current affiliation: Tohkamachi Experimental Station, Forestry and Forest Products Research Institute, Tohkamachi, Japan.

11th Conf. on Cloud Physics, Ogden, UT, Amer. Meteor. Soc., 8.6, https://ams.confex.com/ams/11AR11CP/techprogram/session_12983.htm