Research Article  Open Access
YoungJin Kang, Yoojeong Noh, "Development of Hartigan’s Dip Statistic with Bimodality Coefficient to Assess Multimodality of Distributions", Mathematical Problems in Engineering, vol. 2019, Article ID 4819475, 17 pages, 2019. https://doi.org/10.1155/2019/4819475
Development of Hartigan’s Dip Statistic with Bimodality Coefficient to Assess Multimodality of Distributions
Abstract
In general, although some random variables such as wind speed, temperature, and load are known to have multimodal distributions, input or output random variables are considered to follow unimodal distributions without assessing the unimodality or multimodality of distributions from samples. In uncertainty analysis, estimating unimodal distribution as multimodal distribution or vice versa can lead to erroneous analysis results. Thus, whether a distribution is unimodal or multimodal must be assessed before the estimation of distributions. In this paper, the bimodality coefficient (BC) and Hartigan’s dip statistic (HDS), which are representative methods for assessing multimodality, are introduced and compared. Then, a combined HDS with BC method is proposed. The proposed method has the advantages of both BC and HDS by using the skewness and kurtosis of samples as well as the dip statistic through a link function between the BC values in BC and significance level in HDS. To verify the performance of the proposed method, statistical simulation tests were conducted to evaluate the multimodality for various unimodal, bimodal, and trimodal models. The implementation of the proposed method to real engineering data is shown through case studies. The results demonstrate that the proposed method is more accurate, robust, and reliable than the BC and original HDS alone.
1. Introduction
In the field of engineering, most random variables are treated as having unimodal distributions with one mode so that the estimated probability density function (PDF) or cumulative distribution function (CDF) with unimodality from samples is applied to probabilistic statistical analysis or design methods such as reliability analysis, reliabilitybased design optimization, and robust design [1, 2]. However, it is observed that some random variables such as wind speed, temperature, amplitude of saving tools, shear strength of seafloor sediment, and vehicle weight have multimodal distributions [3–7]. Furthermore, output random variables generated for responses of systems occasionally have multimodal distributions owing to the highly nonlinear performance function or multimodal input variable [8, 9]. To determine the absolute appropriateness of candidate unimodal distributions for given data, goodnessoffit (GOF) tests such as Kolmogorov–Smirnov, Anderson–Darling, etc., can be used. Although GOF tests can identify whether given data are fit to candidate unimodal distributions, they cannot be used to assess the unimodality or multimodality of estimated distributions from data [2]. If users know in advance whether a random variable has a unimodal or multimodal distribution, they can select the appropriate parametric or nonparametric modeling method [1, 2, 10]. However, in many cases, the modality of random variable distributions is unknown. Of course, if a sufficient amount of data is available, it might be possible to assess the modality of data using histograms. The histogram can represent the distribution of numerical data so that unimodality or multimodality can be easily identified. Still, the histogram can be seen as a unimodal or multimodal distribution according to the number of bins, bin size, and starting points. Thus, a statistical measure is needed to numerically quantify the modality of data.
Methods of assessing the unimodality or multimodality of distributions have been developed in the field of statistics, which have been widely applied to various fields such as medical and brain sciences as well as statistical field [11–14]. Wolfe introduced a method using the likelihood ratio with a normal mixture model to assess unimodality [15]. Hartigan and Hartigan developed Hartigan’s dip statistic (HDS) method by introducing the dip statistic to assess unimodality [16]. SAS Institute Inc. proposed the bimodality coefficient (BC) to assess multimodality by using information of the third and fourth statistical moments of data [17]. The Akaike information criterion between onecomponent and twocomponent distribution models (AIC_{diff}) was introduced to determine bimodality by comparing the AIC values of a onecomponent Gaussian mixture model with that of a twocomponent Gaussian mixture model [18]. Freeman and Dale compared three methods—AIC_{diff}, BC, and HDS—for the assessment of dual cognitive process in behavioral science [11]. Although HDS is a more accurate, reliable, and appropriate method for assessing a unimodality or multimodality, it slowly converges to the modality of the true model as sample size increases, and its results are sensitive to the significance level. On the other hand, BC is less reliable than HDS but quickly converges to the true model for multimodal distribution. AIC_{diff} is quite sensitive to the normality assumption and bimodality; thus, it detects nearly all distributions as a bimodal distribution [11, 19].
Freeman and Dale revealed that AIC_{diff} was worse than BC and HDS because it is highly sensitive to bimodality unlike BC and HDS and that HDS was a more suitable method than BC [11]. Kang and Noh compared BC and HDS through simulation tests by considering the randomness of data, sample size, and shape of the distribution function. They also observed that HDS was more accurate and reliable than BC [19]. However, HDS also has some limitations in estimating modality according to significance level; hence, an advanced HDS method must be developed.
The modality assessment method have been recently applied into various research fields; the BC has been employed in the psychiatry science, public opinion, chemistry, and physics [11, 13, 20–22], and HDS has been implemented in the brain science, biological science, microbiology, and ecology [23–26]. Even though there are many applications with multimodality in the engineering fields, the methods of evaluating the multimodality have been rarely applied, and graphical methods such as histograms have been simply employed [3–9, 27–32].
This study proposed a new modality assessment method (HDS with BC), which was first applied in the field of engineering, by combining the existing modality estimation methods HDS and BC. The HDS with BC defines the relationship between BC values and the significance level as a function in order to select the significance level appropriate for estimating the modality of the distribution, which results in HDS determining the accurate modality at the optimized significance level. HDS with BC has the advantages of both BC and HDS by doublechecking the modality of data, preventing errors caused by mistakenly selecting data with unimodality as a multimodal distribution or vice versa, which are commonly made with BC and HDS.
It was verified that the proposed method is more accurate, reliable, and quickly converges to the true unimodality or multimodality through the assessment of multimodality of unimodal, bimodal, and trimodal distributions in simulations and case studies of real measurements and engineering data. Accordingly, the proposed HDS with BC method can improve the accuracy of quantification and propagation such as calculations of variation in dynamic stiffness of rubber mounts [4], reliabilitybased design optimization of helicopters [8], reliability assessment [32], and fatigue reliability assessment for steel decks [7] through correct identification of data modality. It can be highly likely to be used in the engineering field.
2. Method to Assess Multimodality
2.1. Bimodality Coefficient
The BC is an empirical method which assumes that a bimodal distribution will have high skewness, low kurtosis, or both [17]. BC is a very simple and easy method to assess unimodality or multimodality since BC values are easily calculated from only sample size, skewness, and excess kurtosis. The BC value is calculated as [17]where is the sample size, is the skewness, and is the excess kurtosis estimated from given data. Since both skewness and kurtosis are sensitive to sample bias, the corrected skewness and kurtosis for sample bias are used as [33, 34]
If (BC value for uniform distribution), the data are considered to follow a unimodal distribution. Otherwise, they follow a bimodal or multimodal distribution [11, 17]. BC has a normalized value on , and it is a deterministic value; thus, it is easy to utilize and intuitively evaluate the unimodality or multimodality of distribution.
2.2. Hartigan’s Dip Statistic
HDS is based on the hypothesis test that the given data have more than one mode in their distribution using the dip statistic. The “dip” refers to the maximum difference between the empirical distribution function and unimodal distribution function that minimizes that maximum difference where a uniform distribution is used as a reference unimodal distribution. The dip for unimodal samples converges asymptotically to zero; otherwise, it has a positive value [16]. The mathematical definition and computational algorithm of dip are described in more detail in the original papers [16, 35]. In this study, the dip algorithm for testing multimodality is briefly explained. HDS requires two dip statistic values, i.e., the dip of given samples () and r^{th} uniform samples () randomly generated from a uniform distribution, where the number of uniform samples (bootstrap samples) should be equal to the number of given data. In other words, is the dip between eCDF and estimated uniform distribution from given samples and is the dip from arbitrary uniform samples for bootstrap samples. In this paper, the arbitrary samples are generated from . Both dips are iteratively compared using equation (3) until r becomes the total number of bootstrap samples. Finally, the pvalue can be calculated aswhere is the total number of bootstrap samples and is the indicator function for each uniform sample. If , then is closer to a unimodal distribution than a uniform distribution; is one, and if , then is closer to a multimodal distribution.
The null hypothesis of HDS is that a distribution is unimodal. Thus, if the value is greater than or equal to a specific significance level (), then HDS assesses a sample as having a unimodal distribution; otherwise, it has a multimodal distribution. Notice that HDS requires more computation time than BC in calculating the value because the bootstrapping process is repeatedly performed to calculate the value.
3. Hartigan’s Dip Statistic with Bimodality Coefficient
Both BC and HDS have characteristics that contradict each other. BC can be calculated by a simple mathematical formula using only sample size and statistical moments; thus, it returns a deterministic result. In addition, it can quickly converge to the multimodality of the true model, especially with strong multimodality. However, it is too sensitive to the third and fourth statistical moments; thus, the results could inconsistently converge to the true modality as the sample size increases. Although the true model is a unimodal distribution, BC frequently assesses it as multimodal, i.e., it makes a Type I error when the null hypothesis is that a distribution is a unimodal distribution. This is a critical disadvantage in the assessment of unimodality for unimodal distributions. On the other hand, HDS tests the multimodality by a dip statistic; thus, it returns the probability result (value) based on a specific significance level, and it does not depend on the statistical moments of samples unlike BC. Therefore, it tends to converge to the true modality as the number of samples increase. It is also reliable and stable and thus widely used for the assessment of unimodality or multimodality [11, 12, 14, 26].
Even if HDS is more useful than BC, HDS is empirical, subjective, or very difficult to use because users are required to select the appropriate significance level based on their experience or amount of given data [36, 37]. The higher the significance level, the higher the probability of identifying data with unimodality as one with multimodality. On the other hand, the lower the significance level, the higher the probability of identifying data with multimodality as one with unimodality. Accordingly, the appropriate significance level should be used to avoid Type I or Type II errors, which decrease the assessment accuracy according to data modality. Furthermore, the convergence rate of HDS is relatively slower than that of BC for a multimodal model because of repeated bootstrapping and selecting significance level, and small sample sizes, low significance levels, or both are occasionally assessed as unimodal distribution even when the true model is a multimodal distribution, i.e., Type II error occurs. As shown in Table 1, BC and HDS have different advantages and disadvantages in estimating the modality of data. BC is based on statistical moments, whereas HDS is based on the shape of the distribution function using the dip statistic. Both moments and dip statistic are important information to assess the multimodality. Therefore, it is necessary to develop a method of correctly identifying data modality by combining them.
 
H_{0}: X is a unimodal distribution. 
In this study, to compensate for the differences between the two methods, BC is correlated with HDS to determine the optimal significance level in HDS; thus, the proposed method is called HDS with BC (or HDSw/BC). The proposed HDSw/BC method first deterministically calculates BC values. Then, if the BC values are low, i.e., the data are more likely to have unimodality, it uses a low significance level to increase the probability that HDS judges the data as a unimodal distribution. On the other hand, higher BC values indicate the higher probability of data having multimodality; thus, HDS uses a high significance level to increase the probability that HDS judges the data as a multimodal distribution. Although HDSw/BC uses the skewness and kurtosis, it is not sensitive to the skewness and kurtosis of samples since these are not directly used to assess multimodality and are only used to select the optimal significance level in HDS. Consequently, the proposed method can achieve more reliable, stable, useful, and accurate performance than BC and HDS by doublechecking the unimodality and multimodality of data, thereby adaptively reducing either Type I or II errors based on information on moments and distributions.
For this, it is necessary to define the relationship between the significance level and BC values to determine the optimal significance level for HDSw/BC. In HDS, if the significance level increases, the probability of assessing multimodality increases, i.e., Type II error decreases but Type I error increases; otherwise, the probability of assessing unimodality increases, i.e., Type I error decreases but Type II error increases (see Appendix A). If the BC value is large, it means that multimodality is clearly observed in data; hence, Type II error needs to be reduced using a high significance level. Otherwise, unimodality is observed in data; hence, Type I error needs to be reduced using a low significance level. The relationship between the significance level and BC values can be defined such that the optimal significance level is selected according to BC values. The significance level is subjectively determined by users in most cases and can vary from 5% to 32%. A significance level of 5% is often used as the default value [36–38] and was used in the original HDS paper [11, 16, 35, 39]. The significance level of 32% is the corresponding value to the 1sigma of 689599.7 rule [40]. Four types of functions, i.e., linear, quadratic, exponential, and irrational functions, were used to define the relationship between the BC values and significance level as shown in Table 2, where is the lower significance level with , is the upper significance level with , and BC is the BC values. Figure 1 shows the graphs of the four types of functions according to BC values.

Consequently, HDSw/BC can automatically select the optimal significance level using the information on sample size, skewness, and kurtosis by BC and then assess multimodality using a dip statistic with the optimal significance level. Even though the HDSw/BC uses BC in an empirical manner, it finally and reasonably assesses the modality through the hypothesis test using the dip statistic. It was proved by the original paper for HDS where the dip asymptotically converges to a positive value when samples are multimodal distribution as sample size increases; otherwise, it converges to zero [16, 35]. It assesses a multimodal distribution if the value is less than the significance level, which varies by BC value in equation (1) according to the four types of functions; otherwise, it assesses a unimodal distribution. Since the estimation of modality changes depending on function type, it needs to be tested for the four function types. The simulation tests and their results for different function types are explained in Appendix B. Assuming that the population distribution is unimodal, bimodal, and trimodal, HDS was tested at the optimum significance level obtained from the four types of functions. Quadratic and exponential functions showed similar performance owing to their similar function shapes, while the results of the linear function lies between that of the quadratic and irrational functions. Thus, only the quadratic and irrational functions were used in this study.
4. Simulation to Assess Multimodality
To compare the performance of BC, HDS, and HDSw/BC, statistical simulation tests were conducted to assess the multimodality of samples when the true models are unimodal, bimodal, and trimodal distributions. The samples were randomly generated from the true model n = 5, 7, 10, 20, 30, 50, 100, 200, and 300 with 1,000 repetitions; thus 1,000 sample sets were generated for each sample size (n). Subsequently, the assessment of multimodality was conducted for each sample set as n increased, and the simulation results were compared using different types of multimodality measures.
4.1. Unimodal Distribution
Modality simulation tests were conducted for two unimodal distributions: Birnbaum–Saunders (BS) distribution and lognormal (LN) distributions, where is a skewed distribution with extremely heavy tails and is an extremely skewed distribution with heavy tails. Figure 2 shows PDFs of two true unimodal distributions, and Figure 3 shows the percentages of identifying unimodality (accuracy) and multimodality (error) according to the number of samples (n) with 1,000 repetitions. In Figure 3, the light gray bar and dark gray bar indicate the accuracy and error, respectively. The first, second, third, and fourth bars indicate the results using BC and HDS with , HDSw/BC with the quadratic model (HDSw/BC(Quad)), and HDSw/BC with irrational model (HDSw/BC(Irr)), respectively.
(a)
(b)
When the BS distribution is the true model, HDS and HDSw/BC(Quad) correctly identified the unimodality of the true model for most sample sizes except ; thus, it has the highest accuracy among all measures. On the other hand, HDSw/BC(Irr) has the lowest accuracy for , while BC has the lowest accuracy for . The error of HDSw/BC(Irr) consistently decreased as n increased, while that of BC increased for and then decreased for . This is because sample skewness and kurtosis are often incorrectly estimated with large variations since BC is sensitive to the sample size and the distribution of sampled data. BC values increased for high skewness or low kurtosis of samples or both; especially, skewness significantly affects the increase of BC values (see Appendix C). HDSw/BC(Irr) uses a higher significance level than HDSw/BC(Quad), which often leads to Type I error because it identifies a unimodal distribution as a multimodal one. However, as the number of samples increased, the accuracy of BC increased or decreased without a consistent tendency owing to inaccurate estimation of skewness and kurtosis, while HDSw/BC(Irr) consistently decreased Type I error. HDSw/BC(Quad) has a slightly larger significance level than HDS; hence, its accuracy is slightly lower than that of HDS with . However, their performances are similar, are more accurate than HDSw/BC(Irr), and are more reliable than BC regardless of sample size.
When the true model is a lognormal distribution with LN(3,0,8), similar to BS(50,0.4), HDS and HDSw/BC(Quad) show similar performance and their errors are the lowest for most sample sizes except in Figure 3. However, the error of BC increased as n increased and thus that of BC is the highest for . For LN distribution with high skewness, the disadvantage of BC is more clearly observed because of the large estimation error of sample skewness, especially for smallsized samples. For skewed unimodal distributions, the HDS performance that does not rely on skewness and kurtosis becomes more powerful than BC. Since HDSw/BC method indirectly uses the sample skewness and excess kurtosis to determine the optimal significance level, the results of the proposed method are rarely affected by sample skewness and kurtosis, unlike the BC method. In particular, the performance of HDSw/BC(Quad) is almost similar to the original HDS since the estimated optimal significance levels are close to the significance level of the original HDS. Accordingly, HDS and HDSw/BC(Quad) methods are the most accurate, robust, and reliable.
4.2. Multimodal Distribution
In this section, tests are implemented for bimodal and trimodal distributions. Figure 4 illustrates two bimodal distributions, a mixture of two normal distributions with 0.5 N(150,10)&0.5 N(200,10), and a mixture of two lognormal distributions with 0.5LN(3,0.4)&LN(4,0.2).
Figure 5 depicts the results of the multimodality assessment. In Figure 5(a), when a normal mixture model is a true model, the accuracy of identifying the multimodality increases as increases. Then, finally all the methods correctly identify the multimodality of the true model for . HDSw/BC(Irr) has the highest accuracy for all the sample sizes, followed by HDSw/BC(Quad); BC and HDS have the lowest one for . Especially, the accuracy using HDSw/BC(Irr) remarkably increases as increases, whereas the accuracy using HDS gradually increases and BC does not tend to identify the modality.
(a)
(b)
For the two lognormal mixture models depicted in Figure 5(b), although the overall tendency is similar to that of the normal mixture model, the convergent rate to the true modality is lower than that of the normal mixture model because of the weak bimodality of the lognormal mixture model. The true skewness of the two lognormal mixture models is 0.26; however, the original HDS does not use the skewness, whereas the proposed methods use it to determine the optimal significance level for HDS. In addition, although the number of sample sizes increases, BC cannot converge to the true modality and its convergence rate is extremely low because of the low skewness and high kurtosis.
Subsequently, the multimodality test is conducted for the trimodal distributions, mixture of normal distributions with 0.4 N(100,20), 0.3 N(200,20), and 0.3 N(300,20), and mixture of Weibull distributions with 0.4 W(100,4), 0.3 W(200,10), and 0.3 W(300,12). Figure 6 depicts the PDFs of two trimodal distributions, and Figure 7 depicts the accuracy and errors obtained using BC, HDS, HDSw/BC(Quad), and HDSw/BC(Irr) over 1,000 repetitions. In Figure 7, a tendency for the two trimodal distributions is similarly observed in the one for the bimodal distributions; HDSw/BC(Irr) converges to the true multimodality most rapidly, followed by HDSw/BC(Quad); however, their convergence rates are quite different from each other. For the normal mixture model, the convergent rate of BC is the lowest for and that of HDS is the lowest for . Because the BC value of the true model is 0.630, which is high for a multimodal distribution, BC quickly converges to the true modality as n increases. However, for the trimodal distribution with Weibull distributions, BC does not converge to the true modality because its true BC value is 0.585, which is similar to that of a uniform distribution. The multimodality of the Weibull mixture model is not as severe as that of the normal mixture model; hence, the modality assessment results of the normal mixture model are better than those of the Weibull mixture model.
(a)
(b)
4.3. Summary of Simulation Tests
In summary, the unimodal and multimodal distributions are tested when the null hypothesis distribution is a unimodal distribution. Table 3 presents the accuracy (%), which is calculated as the ratio of the number of correct modality identifications to the total repletion number for unimodal distributions, where the values in bold indicate the highest accuracy for each sample size and the values in italics indicate the lowest accuracy. In Table 3, the accuracy of BC is the lowest for unimodal models, except when and for BS and LN models, respectively. BC is extremely inaccurate for the LN distribution because the skewness and kurtosis of the population are quite high; thus, the sample skewness and kurtosis are incorrectly estimated. In addition, as mentioned before, BC is inconsistent for both the unimodal models as n increases. HDS is the most accurate, followed by HDSw/BC(Quad) and HDSw/BC(Irr). Although the conventional HDS is more accurate than both HDSw/BC methods, the accuracy of HDSw/BC(Quad) is almost equal to that of HDS and the accuracies of both the methods are higher than 90% for .

Table 4 presents the accuracy for the multimodal distributions. In Table 4, the accuracy of assessing multimodality using HDSw/BC(Irr) is always the highest; however, that using BC or HDS is the lowest according to the number of samples (n). The results reveal the limitations of the BC and conventional HDS methods in assessing the multimodality for the multimodal distributions, whereas HDSw/BC(Quad) and HDSw/BC(Irr) are considerably accurate, convergent, and consistent methods in comparison with the BC and HDS methods.

Therefore, using BC is extremely risky in assessing the multimodality for highly skewed or heavytailed unimodal distributions and multimodal distributions, whereas HDS and HDSw/BC methods are more reliable and useful. Although HDSw/BC(Quad) and HDSw/BC(Irr) yield high accuracies of multimodality assessment for multimodal distributions, the accuracies of both the methods are lesser than that of HDS. Figures 8 and 9 depict the difference in accuracies between the original HDS and HDSw/BC methods. If HDSw/BC identifies the true modality with a higher accuracy than HDE, then it is positive; otherwise, it is negative. In Figure 8, the accuracies of HDSw/BC(Quad) for the BS and LN distributions are lower than those of the original HDS by approximately 1–3% for ; however, the accuracies of both the methods are almost the same for . On the other hand, the accuracy of HDSw/BC(Quad) for the multimodal distributions is higher than that of HDS by approximately 1–23%. In Figure 9, the accuracy of HDSw/BC(Irr) for unimodal distributions for is lower than that of HDS by approximately 1–11%; otherwise, the accuracies of both the methods are almost the same. On the other hand, the accuracy of HDSw/BC(Irr) for multimodal distributions is considerably higher than that of HDS by approximately 1–39%.
In summary, for the unimodal distributions, HDS with is slightly more accurate than HDSw/BC; however, for multimodal distributions, it results in a considerably lower accuracy than HDSw/BC and even BC. Because HDSw/BC(Quad) has a high accuracy, which is similar to that of HDS for a unimodal distribution, and HDSw/BC(Irr) is the most accurate for a multimodal distribution, they can be said to be superior to BC and HDS with . HDS may reduce Type I errors by selecting a lower significance level with ; it may also reduce Type II errors by selecting a higher significance level with . However, it is difficult to select an appropriate significance level without knowing the unimodality or multimodality.
In contrast, HDSw/BC method does not allow users to select the significance level but calculates an optimum significance level so that the modality of data can be easily identified. It should be noted that the performance of HDSw/BC may vary according to the types of functions that define the relationship between BC values and significance levels. Similar to the way HDS method reduces Type I or Type II errors by changing the significance levels, HDSw/BC can reduce them by using either a quadratic function or an irrational function. If the number of data points is insufficient to determine the modality of data (), either HDSw/BC(Quad) or HDSw/BC(Irr) can be used according to the type of errors that needs to be reduced, resulting in higher accuracies being achieved than those using HDS with fixed significance levels. If sufficient data are available, then the performances of HDSw/BC(Quad) and HDSw/BC(Irr) are similar for unimodal distributions and the performance of HDSw/BC(Irr) is better than that of HDSw/BC(Quad) for multimodal distributions. Thus, users are recommended to employ either HDSw/BC(Quad) or HDSw/BC(Irr) for insufficient data according to the type of errors that need to be reduced and employ HDSw/BC(Irr) for sufficient data.
5. Case Study
To demonstrate the performance of the proposed method, HDSw/BC methods are implemented to actual engineering data by comparing the results of the BC, original HDS, HDSw/BC(Quad), and HDSw/BC(Irr) methods similar to the simulation tests described in Section 4. In this study, a simulation was performed for data on temperature, wind speed, and shift factor of rubber material.
5.1. Temperature Data
Temperature can significantly affect the performance of temperaturedependent mechanical systems, components, or materials such as viscoelastic materials, as well as the health status of battery systems [4, 41]. Temperature can follow a unimodal or multimodal distribution according to different operating environments, and thus, it requires an appropriate statistical model. In this section, an assessment of multimodality is conducted for two temperature datasets collected from two sites in South Korea.
The first dataset contained hourly data collected at Daegwallyeong in South Korea from 2014 to 2018. The total number of data points is 42,864 which is assumed as a population dataset [42]. Test samples were randomly sampled from the population dataset for n = 10∼8,000 with 1,000 repetitions for each sample size. Figure 10 shows the histogram of population data of temperature where the temperature data distribution is slightly skewed to right and has a weak bimodal shape. Figure 11 shows the accuracy and error for assessing modality using the BC, HDS, and HDSw/BC methods. All methods except BC converge to the true model’s modality as increases. The estimated accuracies using HDSw/BC(Irr) increase most rapidly as increases followed by HDSw/BC(Quad) and HDS; however, BC does not correctly assess the modality as the true skewness is as small as −0.282; thus, the BC value of population is 0.4843.
Another dataset was collected hourly at Seoul, South Korea, in 2007, and the number of data points was 8,760 [42]. Samples were randomly generated for n = 10–2,000. Figure 12 shows the histogram of the temperature data at Seoul, and it has stronger bimodality than Daegwallyeong. Figure 13 summarizes the results of test for BC, HDS, and two HDSw/BC methods. The similar tendency was observed in the Seoul data, similar to the Daegwallyeong data; therefore, the original HDS and HDSw/BC methods tend to correctly identify the true modality except the increase in BC increase as n increases. However, the convergence rates for the Seoul data are considerably faster than those for Daegwallyeong data because Seoul data have stronger bimodality. Accordingly, the proposed methods are shown to be accurate, consistent, and rapidly converging to true modality than the BC and original HDS for two temperature datasets, a weak and strong bimodal distribution.
5.2. Wind Speed Data
Wind speed significantly affects the power, efficiency, and reliability of the wind turbine system. Although wind speed has been generally fitted to a unimodal distribution such as normal, Weibull, and Rayleigh distributions, it was observed that wind speed sometimes has a multimodal distribution in some areas, and wind speed was modeled by a nonparametric modeling method using the kernel density estimation (KDE) as a multimodal distribution. Hu et al. also claimed that KDE using a multimodal model for some wind speed data is more accurate than a unimodal model such as normal, Weibull, and Rayleigh models [32].
Therefore, in this study, the wind speed data at Daegwallyeong in South Korea from 2014 to 2018 were used, and the number of data points was 43,072, which was used as the population [42]. Samples were randomly generated from the population such as the temperature simulation for n = 10–8,000 with 1,000 repetitions. Figure 14 shows the histogram of wind speed data. It can be seen as a unimodal distribution as the Weibull or Rayleigh by smoothing a fluctuated frequency or as a multimodal distribution because there seem to exist explicitly different modes. Although data are sufficient to evaluate the multimodality from the histogram of the data, it is difficult for the users to quantitatively determine the unimodality or multimodality for the wind speed distribution from the histogram of the data. For this, BC, HDS, and HDSw/BC methods are used. As shown in Figure 15, HDS and two HDSw/BC methods tend to converge to the true modality (multimodality) with the increase in the number of samples. In contrast, BC converges to the unimodality unlike the temperature data because it assesses the population for wind speed data as a unimodal distribution with BC = 0.5028 because of the excess kurtosis (1.315), regardless of the large skewness (1.082). This assessment result shows that BC, which only uses the skewness and kurtosis without considering the distribution shape, has a limitation in assessing the multimodality for multimodal distribution with considerably high kurtosis.
5.3. Shift Factor of Rubber Material
Using the temperature data in Section 5.1, shift factor data related to the dynamic stiffness of rubber material can be obtained as [4]where and are the Mooney–Rivlin coefficients, is the reference temperature, and is the ambient temperature. In this study, EPDM40, which is widely used in antivibration rubber, was used, with and as and , respectively [43]. Moreover, the aforementioned hourly temperature data measured at Seoul in 2007 were used [4, 42]. Since the shift factor is directly related to the ambient temperature using equation (4), the statistical model of the shift factor needs to be modeled using the ambient temperature data. Once the modality of the shaft factor is assessed and modeled using the appropriate statistical modeling methods depending on its modality, the statistical model of the shift factor can used to analyze the dynamic stiffness of the rubber material.
The total number of temperature data is 8,760; thus, the shift factor data are calculated by using these datasets. Figure 16 shows the histogram of the shift factor of EPDM40 with 8,760 data, which is assumed as a population dataset. The shift factor has a strong bimodality like the temperature at Seoul but with different shapes. Samples are randomly sampled from the population for n = 10–2,000 with 1,000 repetitions. Figure 17 shows the test results using the BC, HDS, and both the HDSw/BC methods. A similar tendency is observed like the temperature and wind speed cases, but the convergent rate is overall quicker than the temperature and wind speed data because the shift factor has a strong bimodal distribution. As expected, the HDSw/BC(Irr) has the best accuracy for correctly identifying the bimodality of the true model, and the BC has the worst accuracy due to the low skewness of the population.
5.4. A Real Example of the Wind Turbine Power
To show the necessity of assessing multimodality, a simple reliability analysis was conducted to predict power of a wind turbine. Since the available power of the wind turbine varies according to the wind speed and air density and the air density can be calculated from the temperature, the real measured data for the temperature and wind speed at Daegwallyeong described in Sections 5.1 and 5.2 can be used to calculate the reliability of the wind turbine power. The specifications of the wind turbine are presented in Table 5 [44].where the absolute pressure (AP), specific gas constant for dry air (R_{spec}), and power coefficient() are deterministic variables and the air density (), temperature (T), and wind speed () are random variables.

The reliability analysis for the wind turbine power was performed for three cases where the distribution is assumed to have unimodality (Unimodality), known to have multimodality (Multimodality), or the modality is unknown and assessed using HDSw/BC (Estimated modality using HDSw/BC). For the temperature and wind speed, Gaussian and Weibull distributions are often used as unimodal distributions, respectively, so that they were fitted to the data of the temperature and wind speed, respectively [32]. If the multimodality is known, the PDFs are estimated using the KDE with crossvalidation bandwidth methods to represent the multimodality. In HDSw/BC, the distributions are estimated as unimodal or multimodal by obtained results of HDSw/BC. Figure 18 shows the probabilities that the available power is greater than a reference value, 50,000, which is the power at an average wind speed with 3 m/s in the site. In the multimodality case, since the true modality is used, the obtained probabilities are the most accurate as expected, and those using the unimodal distributions are the most inaccurate because of incorrect estimation of the PDFs. However, when the modality is assessed using HDSw/BC(Irr), the results are more accurate and quickly converge to the true models rather than the unimodal case as sample sizes increase. Thus, if the modality is unknown, it is indeed necessary to assess the modality of the distribution.
5.5. Summary of Case Studies
For comparison of the performance of each method, Figure 19 shows the plots of accuracy to assess the multimodality for temperature, wind speed, and shift factor of the rubber material. HDSw/BC(Quad) and HDSw/BC(Irr) are more accurate, consistent, and rapidly converge to true modality than the original HDS and BC for all of the case studies. This is because the proposed methods not only use a sample skewness and kurtosis like the BC but also a dip statistic like the original HDS. Since BC directly uses the skewness and kurtosis of the samples, it is too sensitive to the skewness and kurtosis; thus it could not be converged to the true modality for the bimodal distribution. However, since HDSw/BC methods indirectly use the skewness and kurtosis, they can be less dependent of the skewness and kurtosis. The proposed methods only use the sample skewness and kurtosis to determine the optimal significance level by modeling the relationship between the significance level and BC value for the skewness and kurtosis. Consequently, both the proposed methods well integrate the merits of the original HDS and BC; thus, these are more powerful than the existing methods, in particular, HDSw/BC(Irr) is the most powerful method for assessing multimodality of multimodal distribution.
(a)
(b)
(c)
(d)
In summary, regarding the accuracy of the HDSw/BC(Irr), to achieve over 90% accuracy, the number of temperature data at Seoul and shift factor should be at least n ≥ 500, and the number of wind speed data should be n ≥ 1,000, and the number of temperature data at Daegwallyeong needs to be n ≥ 4,000. It means that the accuracy of assessing the multimodality relies on the degree of multimodality and the quality of data. In particular, the required number of data can be considerably varied according to the degree of multimodality. Although HDSw/BC(Irr) assesses the multimodality more accurately than the BC and HDS, it still requires significantly large data. However, as a result of the simulation tests described in Section 4 and case studies described in Section 5, it is observed that the Type II error occurs more often than the Type I error in HDSw/BC methods due to the insufficient data, regardless of the unimodality or multimodality when the null distribution is a unimodal distribution. For example, if there is a lack of data whose true model is a unimodal distribution, HDSw/BC methods invariably assess them with unimodality despite the unsatisfactory data to evaluate the unimodality. If there are insufficient data whose true model is a multimodal distribution, HDSw/BC methods frequently assess them with the unimodality. However, if there is lack of data whose true model is multimodality but the data are sufficient to represent a multimodality characteristic, HDSw/BC methods assess them with multimodality. In addition, reliability analysis of the wind turbine power was performed to observe how the modality assessment using HDSw/BC affects the reliability analysis. It was found that HDSw/BC yielded accurate reliability analysis results than those obtained using the unimodality assumption when the true modality is unknown.
6. Conclusion
In this paper, bimodality coefficient (BC) and Hartigan’s dip statistic (HDS), which are the existing methods for assessing unimodality or multimodality of distributions, were studied and their limitations were discussed. To overcome the contrast limitations of the BC and HDS, HDSw/BC was proposed. The proposed HDSw/BC uses the BC values calculated from sample skewness and kurtosis to determine the optimum significance level through quadratic or irrational functions and finally identifies the modality of the data using a dip statistic of the original HDS determined by the distribution shapes. Thus, HDSw/BC methods are more reliable and stable than the BC and more useful and accurate than the original HDS. The quadratic function affects the significance level to produce the smallest significance level; thus, HDSw/BC(Quad) has the highest probability to judge the data as unimodality. However, the irrational function yields the largest significance level; thus, HDSw/BC(Irr) has the highest probability to judge the data as a multimodality.
To verify the performance of the proposed methods, the statistical simulation tests to evaluate unimodality or multimodality were carried out for assumed true unimodal, bimodal, and trimodal distributions. To verify the applicability of the proposed method to real applications, real engineering data such as temperature, wind speed, and shift factor of a rubber material were used to identify their modality. Through the simulation tests and case studies, it was shown that HDSw/BC(Quad) is more accurate, robust, reliable, and quickly converges to true modality than the BC and original HDS for both unimodal and multimodal distributions; HDSw/BC(Irr) is the most accurate method and has a high convergence rate to true modality, especially for multimodal distributions. HDSw/BC has, however, slightly lower accuracy than HDS for the unimodal distributions, but the error is very small. Accordingly, if users strongly need to avoid Type I error, i.e., incorrect identification of unimodality as multimodality, HDSw/BC(Quad) is recommended; otherwise, to avoid Type II error, i.e., incorrect identification of multimodality as unimodality, HDSw/BC(Irr) is recommended.
In this paper, the proposed methods are implemented to assess the unimodality or multimodality; however, the results of assessing the modality will affect the statistical modeling and reliability analysis. Accordingly, in future work, the effect of modality of data on statistical modeling and reliability analysis will be investigated. Furthermore, the BC, HDS, and HDSw/BC can only assess either unimodality or multimodality and they cannot evaluate whether samples have bimodal or trimodal or highly modal distributions. Therefore, a method of determining the number of modes will be developed in the future.
Appendix
A. Assessment Results of HDS according to Significance Level
In this section, results obtained using HDS are varied according to a significance level to assess a multimodality. Figure 20 shows the accuracy and error of the assessment results using HDS methods with different significance levels (α). The accuracy of HDS is varied according to the significance level; in particular, it is considerably sensitive to the significance level for multimodal distribution.
(a)
(b)
B. Comparison between Four HDSw/BC Methods
In this section, the four models (HDSw/BC with quadratic (Quad), exponential (Exp), linear (Lin), and irrational (Irr)) are compared. Figure 21 shows the results using HDSw/BC methods for a unimodal distribution with BS(50,0.4), a bimodal distribution with 0.5 N(150,10)&0.5 N(200, 10), and a trimodal distribution with 0.4 N(100, 20), 0.3 N(200, 20), & 0.3 N(300, 20). The performance of HDSw/BC(EXP) is similar to that of HDSw/BC(Quad) and that of HDSw/BC(Lin) is between those of HDSw/BC(Quad) and HDSw/BC(Irr).
(a)
(b)
(c)
C. Skewness, Kurtosis, and BC for BS Distribution
Figures 22–24 show the sample skewness, excess kurtosis, and bimodality coefficient values to clarify the varying BC values according to skewness and kurtosis as increases for BS(50,0.4) distribution. In Figures 22 and 23, the squares of skewness and excess kurtosis are 1.392 (dashdotted line in Figure 22) and 2.33 (dashdotted line in Figure 23), respectively, and the square of the sample skewness and kurtosis converge to the true values of population as n increases. In Figures 23 and 24, the convergent rate of kurtosis of the samples increases for n ≥ 30; therefore, BC values decrease and converge to the population’s BC, 0.449.
Data Availability
A registry of research data for case studies such as temperature and wind speed can be found at https://data.kma.go.kr.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2018R1A6A3A01012213) and the Human Resources Development program (no. 20164030201230) of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) grant funded by the Ministry of Trade, Industry and Energy.
References
 Y. Noh, K. K. Choi, and I. Lee, “Identification of marginal and joint CDFs using Bayesian method for RBDO,” Structural and Multidisciplinary Optimization, vol. 40, no. 1–6, pp. 35–51, 2010. View at: Publisher Site  Google Scholar
 Y.J. Kang, O.K. Lim, and Y. Noh, “Sequential statistical modeling method for distribution type identification,” Structural and Multidisciplinary Optimization, vol. 54, no. 6, pp. 1587–1607, 2016. View at: Publisher Site  Google Scholar
 J. Zhang, S. Chowdhury, A. Messac, and L. Castillo, “A multivariate and multimodal wind distribution model,” Renewable Energy, vol. 51, pp. 436–447, 2013. View at: Publisher Site  Google Scholar
 D. H. Lee and I. S. Hwang, “Analysis on the dynamic characteristics of a rubber mount considering temperature and material uncertainties,” Journal of the Computational Structural Engineering Institute of Korea, vol. 24, no. 4, pp. 383–389, 2011. View at: Google Scholar
 D. Choi, S. J. Kim, and Y. T. Oh, “Feature analysis based on beta distribution model for shaving tool condition monitoring,” Transactions of the Korean Society of Mechanical Engineers A, vol. 34, no. 1, pp. 11–18, 2010. View at: Publisher Site  Google Scholar
 J.S. Choi, S. Hong, S.B. Chi et al., “Probability distribution for the shear strength of seafloor sediment in the KR5 area for the development of manganese nodule miner,” Ocean Engineering, vol. 38, no. 1718, pp. 2033–2041, 2011. View at: Publisher Site  Google Scholar
 Y. Liu, H. Zhang, Y. Liu, Y. Deng, N. Jiang, and N. Lu, “Fatigue reliability assessment for orthotropic steel deck details under traffic flow and temperature loading,” Engineering Failure Analysis, vol. 71, pp. 179–194, 2017. View at: Publisher Site  Google Scholar
 S. Kim, S. Jun, H. Kang, Y. Park, and D.H. Lee, “Reliability based optimal design of a helicopter considering annual variation of atmospheric temperature,” Journal of Mechanical Science and Technology, vol. 25, no. 5, pp. 1095–1104, 2011. View at: Publisher Site  Google Scholar
 W. Lim, J. Jang, and T. H. Lee, “Comparative study of reliability analysis methods for discrete bimodal information,” Transactions of the Korean Society of Mechanical Engineers A, vol. 37, no. 7, pp. 883–889, 2013. View at: Publisher Site  Google Scholar
 Y.J. Kang, Y. Noh, and O.K. Lim, “Kernel density estimation with bounded data,” Structural and Multidisciplinary Optimization, vol. 57, no. 1, pp. 95–113, 2018. View at: Publisher Site  Google Scholar
 J. B. Freeman and R. Dale, “Assessing bimodality to detect the presence of a dual cognitive process,” Behavior Research Methods, vol. 45, no. 1, pp. 83–97, 2013. View at: Publisher Site  Google Scholar
 G. Giannopoulos, G. Oudatzis, G. Paterakis et al., “Red blood cell and platelet microparticles in myocardial infarction patients treated with primary angioplasty,” International Journal of Cardiology, vol. 176, no. 1, pp. 145–150, 2014. View at: Publisher Site  Google Scholar
 B. Hosenfeld, E. H. Bos, K. J. Wardenaar et al., “Major depressive disorder as a nonlinear dynamic system: bimodality in the frequency distribution of depressive symptoms over time,” BMC Psychiatry, vol. 15, no. 1, p. 222, 2015. View at: Publisher Site  Google Scholar
 P. T. Weber and D. E. Walrath, “Modeling noisy resonant system response,” Mechanical Systems and Signal Processing, vol. 84, pp. 615–624, 2017. View at: Publisher Site  Google Scholar
 J. H. Wolfe, “Pattern clustering by multivariate mixture analysis,” Multivariate Behavioral Research, vol. 5, no. 3, pp. 329–350, 1970. View at: Publisher Site  Google Scholar
 J. A. Hartigan and P. M. Hartigan, “The dip test of unimodality,” The Annals of Statistics, vol. 13, no. 1, pp. 70–84, 1985. View at: Publisher Site  Google Scholar
 SAS Institute Inc., SAS/STAT User’s Guide: Version 6, SAS Institute Inc., Cary, NC, USA, 4th edition, 1990.
 G. McLachlan and D. Peel, Finite Mixture Models, Wiley, Hoboken, NJ, USA, 2000.
 Y. J. Kang and Y. Noh, “Comparison study for the assessment of multimodality from statistical data,” in Proceeding of the Korea Society of Mechanical Engineering 2019 Spring Annual Meeting, pp. 12, Jeju, Republic of Korea, 2019. View at: Google Scholar
 Y. Lelkes, “Mass polarization: manifestations and measurements,” Public Opinion Quarterly, vol. 80, no. S1, pp. 392–410, 2016. View at: Publisher Site  Google Scholar
 S. Tempere, M.H. Schaaper, G. Lytra et al., “Molecular determinism of specific anosmia to 2bromo4methylphenol,” Flavour and Fragrance Journal, vol. 32, no. 2, pp. 126–132, 2017. View at: Publisher Site  Google Scholar
 H. Xu, S. O. Skinner, A. M. Sokac, and I. Golding, “Stochastic kinetics of nascent RNA,” Physical Review Letters, vol. 117, no. 12, Article ID 128101, 2016. View at: Publisher Site  Google Scholar
 N. G. Klaasen, C. Kos, A. Aleman, and E. M. Opmeer, “Apathy is related to reduced activation in cognitive control regions during setshifting,” Human Brain Mapping, vol. 38, no. 5, pp. 2722–2733, 2017. View at: Publisher Site  Google Scholar
 E. Zamboni, T. Ledgeway, P. V. McGraw, and D. Schluppeck, “Do perceptual biases emerge early or late in visual processing? Decisionbiases in motion perception,” Proceedings of the Royal Society B: Biological Science, vol. 283, no. 1833, pp. 1–9, 2016. View at: Publisher Site  Google Scholar
 S. G. Schreiber, U. G. Hacke, and A. Hamann, “Variation of xylem vessel diameters across a climate gradient: insight from a reciprocal transplant experiment with a widespread boreal tree,” Functional Ecology, vol. 29, no. 11, pp. 1392–1401, 2015. View at: Publisher Site  Google Scholar
 M. Zimmermann, S. Escrig, G. Lavik et al., “Substrate and electron donor limitation induce phenotypic heterogeneity in different metabolic activities in a green sulphur bacterium,” Environmental Microbiology Reports, vol. 10, no. 2, pp. 179–183, 2018. View at: Publisher Site  Google Scholar
 E. I. GalindoNava, L. D. Connor, and C. M. F. Rae, “On the prediction of the yield stress of unimodal and multimodal nickelbase superalloys,” Acta Materialia, vol. 98, pp. 377–390, 2015. View at: Publisher Site  Google Scholar
 O. ElAtwani, D. V. Quach, M. Efe et al., “Multimodal grain size distribution and high hardness in fine grained tungsten fabricated by spark plasma sintering,” Materials Science and Engineering: A, vol. 528, no. 18, pp. 5670–5677, 2011. View at: Publisher Site  Google Scholar
 B. Chen, Z. Zhong, X. Xie, and P. Lu, “Measurementbased vehicle load model for urban expressway bridges,” Mathematical Problems in Engineering, vol. 2014, Article ID 340896, 10 pages, 2014. View at: Publisher Site  Google Scholar
 G. Mei, Q. Qin, and D.J. Lin, “Bimodal renewal processes models of highway vehicle loads,” Reliability Engineering & System Safety, vol. 83, no. 3, pp. 333–339, 2004. View at: Publisher Site  Google Scholar
 Y. Ji, S. Jiang, Y. Du, and H. M. Zhang, “Estimation of bimodal urban link travel time distribution and its applications in traffic analysis,” Mathematical Problems in Engineering, vol. 2015, Article ID 615468, 11 pages, 2015. View at: Publisher Site  Google Scholar
 B. Hu, Y. Li, H. Yang, and H. Wang, “Wind speed model based on kernel density estimation and its application in reliability assessment of generating systems,” Journal of Modern Power Systems and Clean Energy, vol. 5, no. 2, pp. 220–227, 2017. View at: Publisher Site  Google Scholar
 R. Pfister, K. A. Schwarz, M. Janczyk, R. Dale, and J. Freeman, “Good things peak in pairs: a note on the bimodality coefficient,” Frontiers in Psychology, vol. 4, p. 700, 2013. View at: Publisher Site  Google Scholar
 D. N. Joanes and C. A. Gill, “Comparing measures of sample skewness and kurtosis,” Journal of the Royal Statistical Society: Series D (The Statistician), vol. 47, no. 1, pp. 183–189, 1998. View at: Publisher Site  Google Scholar
 P. M. Hartigan, “Algorithm AS 217: computation of the dip statistic to test for unimodality,” Applied Statistics, vol. 34, no. 3, pp. 320–325, 1985. View at: Publisher Site  Google Scholar
 R. P. Carver, “The case against statistical significance testing, revisited,” The Journal of Experimental Education, vol. 61, no. 4, pp. 287–292, 1993. View at: Publisher Site  Google Scholar
 J. F. Mudge, L. F. Baker, C. B. Edge, and J. E. Houlahan, “Setting an optimal α that minimizes errors in null hypothesis significance tests,” PLoS One, vol. 7, no. 2, Article ID e32734, 2012. View at: Publisher Site  Google Scholar
 D. M. Lane, “Significance level, hyperstat online contesnts,” 2004, http://davidmlane.com/hyperstat/A72117.html. View at: Google Scholar
 F. Mechler, “A direct translation into MATLAB from the original FORTRAN code of Hartigan’s subroutine DIPTST algorithm,” 2002, http://www.nicprice.net/diptest. View at: Google Scholar
 Wikipedia Contributors, “68959937 rule. Wikipedia, the free encyclopedia,” 2019, https://en.wikipedia.org/w/index.php?title=68%E2%80%9395%E2%80%9399.7_rule&oldid=891319550. View at: Google Scholar
 G.w. You, S. Park, and D. Oh, “Realtime stateofhealth estimation for electric vehicle batteries: a datadriven approach,” Applied Energy, vol. 176, pp. 92–103, 2016. View at: Publisher Site  Google Scholar
 Korea Meteorological Administration (KMA), https://data.kma.go.kr, 2019.
 D. H. Lee, S. H. Seo, Y. H. Yun, and J. G. Park, “The stiffness analysis and optimization of the rubber seat considering nonlinear behaviour,” in Proceedings of the Korean Society for Noise and Vibration Engineering Conference, pp. 244–249, Taean, Republic of Korea, 2002. View at: Google Scholar
 Npower, Wind Turbine Power Calculations RWE Npower Renewables, The Royal Academy of Engineering, London, UK, 2019, https://www.raeng.org.uk/publications/other/23windturbine.
Copyright
Copyright © 2019 YoungJin Kang and Yoojeong Noh. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.