59
$\begingroup$

I understand that once we plot the values as a chart, we can identify a bimodal distribution by observing the twin-peaks, but how does one find it programmatically? (I am looking for an algorithm.)

$\endgroup$
3
  • 5
    $\begingroup$ Intriguing question. I've know nothing about this so won't attempt an answer, but googling bimodality+test gives quite a few promising leads. google.co.uk/search?q=bimodality+test $\endgroup$
    – onestop
    Jan 4, 2011 at 13:29
  • $\begingroup$ a programmatic code to test for multimodality is available at www.estima.com and is called MODES.SRC. You can easily convert it to Matlab, R or other softwares. Regards, guido. $\endgroup$
    – user3061
    Feb 5, 2011 at 10:43
  • 1
    $\begingroup$ You can't identify the number of modes (there's no limit to the number of tiny wiggles a distribution could have), but you can get a lower bound. e.g. see here projecteuclid.org/euclid.aos/1176351045 $\endgroup$
    – Glen_b
    Jan 25, 2014 at 11:05

4 Answers 4

35
$\begingroup$

Identifying a mode for a continuous distribution requires smoothing or binning the data.

Binning is typically too procrustean: the results often depend on where you place the bin cutpoints.

Kernel smoothing (specifically, in the form of kernel density estimation) is a good choice. Although many kernel shapes are possible, typically the result does not depend much on the shape. It depends on the kernel bandwidth. Thus, people either use an adaptive kernel smooth or conduct a sequence of kernel smooths for varying fixed bandwidths in order to check the stability of the modes that are identified. Although using an adaptive or "optimal" smoother is attractive, be aware that most (all?) of these are designed to achieve a balance between precision and average accuracy: they are not designed to optimize estimation of the location of modes.

As far as implementation goes, kernel smoothers locally shift and scale a predetermined function to fit the data. Provided that this basic function is differentiable--Gaussians are a good choice because you can differentiate them as many times as you like--then all you have to do is replace it by its derivative to obtain the derivative of the smooth. Then it's simply a matter of applying a standard zero-finding procedure to detect and test the critical points. (Brent's method works well.) Of course you can do the same trick with the second derivative to get a quick test of whether any critical point is a local maximum--that is, a mode.

For more details and working code (in R) please see https://stats.stackexchange.com/a/428083/919.

$\endgroup$
1
  • 3
    $\begingroup$ @venkasub, you might also look at Muller's method for root-finding. It's simpler than Brent's method and "almost" as efficient as Newton's method, but without needing to know the derivatives of the target function. $\endgroup$
    – cardinal
    Feb 5, 2011 at 20:24
27
$\begingroup$

There is a well-known paper by Silverman that deals with this issue. It employs kernel-density estimation. See

B. W. Silverman, Using kernel density estimates to investigate multimodality, J. Royal Stat. Soc. B, vol. 43, no. 1, 1981, pp. 97-99.

Note that there are some errors in the tables of the paper. This is just a starting point, but a pretty good one. It provides a well-defined algorithm to use, in the event that's what you're most looking for. You might look on Google Scholar at papers that cite it for more "modern" approaches.

$\endgroup$
1
  • 5
    $\begingroup$ +1 This paper exploits information about how the number of apparent modes depends on the kernel width (for Gaussian kernels) to select an appropriate width and then continues as described in my reply. $\endgroup$
    – whuber
    Jan 7, 2012 at 15:52
3
$\begingroup$

I came late to the party, but if you are just interested in whether it is multimodal or not, meaning you are not interested in the number of modes, you should look at diptest.

In R the package is called diptest.

$\endgroup$
2
  • 3
    $\begingroup$ The full name is Hartigan's Dip Test of Unimodality. The original paper describes gives this description: The dip test measures multimodality in a sample by the maximum difference, over all sample points, between the empirical distribution function, and the unimodal distribution function that minimizes that maximum difference. The algorithm is not trivial, but you can easily find C, Python, MATLAB, and modern Fortran translations based on the original Fortran implementation in Computation of the Dip Statistic to Test for Unimodality. $\endgroup$
    – mortehu
    Aug 28, 2019 at 0:28
  • 1
    $\begingroup$ MATLAB Code: nicprice.net/diptest and another one - gist.github.com/schluppeck/e7635dcf0e80ca54efb0. $\endgroup$
    – Royi
    Oct 17, 2019 at 1:03
0
$\begingroup$

The definition in wiki is slightly confusing to me. The probability of a continous data set having just one mode is zero. A simple way to program a bimodal distrubiton is with two seperate normal distributions centered differently. This creates two peaks or what wiki calls modes. You can actually use almost any two distributions, but one of the harder statistical opportunities is to find how the data set was formed after combining the two random data distributions.

$\endgroup$
2
  • $\begingroup$ How will you identify the centers of the two normal distributions? $\endgroup$
    – venkasub
    Jan 4, 2011 at 18:41
  • $\begingroup$ @venkasub, the center of normal distribution is its mean. If you know that bimodality comes from two normal distributions you can fit mixture model to estimate the means, which will be two modes. $\endgroup$
    – mpiktas
    Jan 5, 2011 at 9:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.