Finding minimum/maximum peaks in a n-modal distribution

I have distributions that show n-modal behavior. I need to find the values of the largest and smallest modes. For example, in the histogram below I need to find the values representing the yellow lines (first one is around 20 and the last one is around 190). The red ones are not important to me for now. One of the problems is that I cannot guarantee that any given mode has a normal distribution. In fact, I cannot guarantee any distribution at all. Also I cannot know beforehand how many modes I can find in graph.

Is there any analysis that I can do to find those values?

Figure 1: SAMPLE DISTRIBUTION HISTOGRAM

A very long time ago I learned an effective technique in the geological literature. (I apologize for not remembering the source.) It consists of studying the modes of a kernel density estimator (KDE) as the bandwidth is varied.

What happens is that with a very large bandwidth, the data look like a big lump with a single mode. This one uses a bandwidth of 60 and its mode is near 110:

As the bandwidth shrinks, the KDE outlines what the eye sees more closely and more modes appear. This one uses a bandwidth of 10 and has three obvious modes with a fourth just beginning to show near 60:

When the bandwidth shrinks too far, the KDE is too detailed. This one with a bandwidth of 1 has 36 modes:

You can explore this behavior with a "mode trace." For each bandwidth within the full range (from no detail to too detailed) it plots the modes. I have tracked the evolution of each mode and colored them accordingly. For instance, the single mode in the first figure corresponds to the central red line (shaped almost like a question mark); the four modes in the second figure correspond to the four traces rising to a height (bandwidth) of 10; the 36 modes in the third figure correspond to all 36 traces:

It's probably a good idea to use a logarithmic scale for the bandwidth, as shown here.

A glance at the mode trace will indicate how many modes to identify. I have chosen four. To determine their locations, I have found the points where the traces are the most vertical among all bandwidths smaller than the one at which all four modes first appear: at these locations the locations are stable even when the bandwidth changes. It is comforting (but not really essential) that all four locations are found using comparable bandwidths. (One really should take a bit more care in case multiple stable points appear along a trace: I would opt for the one with the largest bandwidth less than the bandwidth at which all the modes appear.)

Having located the modes, we may plot them on the original histogram:

It's then a simple matter to select the extreme modes. The mode trace will tell you how sensitive their locations are to both the number of modes you identify and to the bandwidth you use. In this example it suggests a tendency for the highest mode to grow even greater with smaller bandwidths before it splits into multiple modes, but the other three modes remain relatively stable (their traces remain nearly vertical at low bandwidths).

It doesn't much matter what shape kernel you pick. The original paper suggested using a Gaussian kernel, which I have done here. The use of a Gaussian is not tantamount to any assumption that the peaks will even approximately have Gaussian shapes. Because Gaussians are (infinitely) smooth, so is the KDE, which means you can analyze it with Calculus techniques to your heart's content.

To be perfectly clear, here is a mathematical account of the mode trace. Let the Kernel function $$K$$ have unit area and unique mode at $$0$$ and let the data be $$x_1, \ldots, x_n.$$ The KDE of the data with bandwidth $$h\ge 0$$ is the convolution

$$f(x,h) = \frac{1}{nh}\sum_{i=1}^n K\left(\frac{x-x_i}{h}\right).$$

For each $$h\ge 0,$$ let $$M(h)$$ be the set of modes of the distribution function $$x\to f(x,h).$$ The "mode trace" of the data is the union of $$M(h)$$ as $$h$$ ranges over an interval $$(0, A)$$ where $$A$$ has been chosen so large that $$M(h)$$ contains a unique element for all $$h\ge A.$$

The mode trace has additional structure: it can be decomposed (not necessarily uniquely) into the disjoint union of graphs of continuous partial functions of $$h$$ defined on intervals. This decomposition is maximal in the sense that the only points any two distinct such functions can possibly have in common are at the endpoints of their domains. I have used colors to designate these partial functions.

Apart from selecting the number of modes to use--which depends very much on your concept of the correct resolution at which to analyze your data--everything can be automated. Here is the R code I used to generate sample data, analyze them, and make the figures. Its results will be contained in a dataframe X recording the mode trace and an array modes containing information about the selected modes.

BTW, if you code your own, note that the KDE is obtained most efficiently using the Fast Fourier Transform (FFT). The most efficient method transforms the data once and then multiplies that by a sequence of transformed kernels, inverting each product to produce the KDE. To determine the range of bandwidths to search, make the largest approximately one-quarter the range of the data and the smallest perhaps 3% or 1% of that.

#
# Generate random values from a mixture distribution.
#
rmix <- function(n, mu, sigma, p) {
matrix(rnorm(length(mu)*n, mu, sigma), ncol=n)[
cbind(sample.int(length(mu), n, replace=TRUE, prob=p), 1:n)]
}
mu <- c(25, 60, 130, 190) # Means
sigma <- c(8, 13, 15, 19) # SDs
p <- c(.18, .2, .24, .28) # Relative proportions (needn't sum to 1)
n <- 1e4                  # Sample size
x <- rmix(n, mu, sigma, p)
#
# Find the modes of a KDE.
# (Quick and dirty: it assumes no mode spans more than one x value.)
#
findmodes <- function(kde) {
kde$$x[which(c(kde$$y[-1],NA) < kde$$y & kde$$y > c(NA,kde$$y[-length(kde$$y)]))]
}
#
# Compute the mode trace by varying the bandwidth within a factor of 10 of
# the default bandwidth.  Track the modes as the bandwidth is decreased from
# its largest to its smallest value.
# This calculation is fast, so we can afford a detailed search.
#
m <- mean(x)
id <- 1
bw <- density(x)$$bw * 10^seq(1,-1, length.out=101) modes.lst <- lapply(bw, function(h) { m.new <- sort(findmodes(density(x, bw=h))) # -- Associate each previous mode with a nearest new mode. if (length(m.new)==1) delta <- Inf else delta <- min(diff(m.new))/2 d <- outer(m.new, m, function(x,y) abs(x-y)) i <- apply(d, 2, which.min) g <- rep(NA_integer_, length(m.new)) g[i] <- id[1:ncol(d)] #-- Create new ids for new modes that appear. k <- is.na(g) g[k] <- (sum(!k)+1):length(g) id <<- g m <<- m.new data.frame(bw=h, Mode=m.new, id=g) }) X <- do.call(rbind, args=modes.lst) X$$id <- factor(X$$id) # # Locate the modes at the most vertical portions of their traces. # minslope <- function(x, y) { f <- splinefun(x, y) e <- diff(range(x)) * 1e-4 df2 <- function(x) ((f(x+e)-f(x-e)) / (2*e))^2 # Numerical derivative, squared v <- optimize(df2, c(min(x),max(x))) c(bw=vminimum, slope=vobjective, Mode=f(vminimum)) } # # Retain the desired modes. # n.modes <- 4 # USER SELECTED: Not automatic bw.max <- max(subset(X, id==n.modes)$$bw)
modes <- sapply(1:n.modes, function(i) {
Y <- subset(X, id==i & bw <= bw.max)
minslope(Y$$bw, Y$$Mode)
})
#
# Plot the results.
#
library(ggplot2)
ggplot(X, aes(bw, Mode)) +
geom_line(aes(col=id), size=1.2, show.legend=FALSE) +
geom_point(aes(bw, Mode), data=as.data.frame(t(modes)), size=3, col="Black", alpha=1/2) +
scale_x_log10() +
coord_flip() +
ggtitle("Mode Trace")

ggplot(data.frame(x), aes(x, ..density..)) +
geom_histogram(bins=500, fill="#2E75B2") +
geom_vline(data=as.data.frame(t(modes)),
mapping=aes(xintercept=Mode), col="#D18A4e", size=1) +
ggtitle("Histogram With Modes")

• (+1) The same broad ideas occur under the heading mode tree in mainstream statistical literature. Googling "Minotte mode tree" finds references. – Nick Cox Sep 21 '19 at 7:06
• Can you give me more references about the Mode trace. I could not understand completely. – Lin Sep 21 '19 at 18:56
• R is hard to read. – Lin Sep 22 '19 at 7:49
• Use nodes identified by a bandwidth that is appropriate for your analysis. I hesitate to recommend any automatic method, but a large vertical gap (on the log scale) between the appearance of successive modes would suggest a default cut point. In the example here, the largest gap is between 3 and 11 (approximately), corresponding to four modes. – whuber Sep 23 '19 at 21:18
• (1) "Graph" has a standard meaning in mathematics and I use it in that sense. By definition, a function $f:X\to Y$ is a subset $G_f$ of the cartesian product $X\times Y$ for which $(x,y_1)\in G_f$ and $(x,y_2)\in G_f$ implies $y_1=y_2.$ Outside of a set-theoretic context, $G_f$ is called the "graph" of $f,$ even though it is identical to $f.$ (2) Yes, when $f:X\to Y$ is a partial function, the projection $\pi_X:G_f\to X$ is not necessarily a surjection. (3) This answer already shows a simple example of a decomposition: each curve in the mode trace graphs a partial function of bw. – whuber Sep 30 '19 at 16:08

Here is some fake simulated data that are multimodal. The figure shows three kinds of plots (made with R): (1) histogram [blue], (2) tick marks below the axis, and (3) a kernel density estimator (KDE) of the data [red].

hist(x, prob=T, br=40, col="skyblue2"); rug(x)
lines(density(x), col="red")


I think the main difficulty with your plot is that it confuses local ties with overall modes. My histogram bars are also of irregular heights because I (purposely) used too many bars (by using parameter br=40). However, even with the best choice of numbers of bars, histograms are not the best kind of device for finding modes.

The 'bandwidth' of the KDE can be adjusted for finding modes. Above, I used the default bandwidth. Maybe a slightly narrower bandwidth would have worked a little better. (I used parameter adj=.5 below. You can read the R documentation of density to see how to change bandwidth and types of kernels.)

hist(x, prob=T, ylim=c(0, .01), col="skyblue2"); rug(x)


If you like, you can get a printout of around 500 heights of the KDE, then scan the numerical list left to right for increases and decreases in order to locate modes. (If scanning has to be automated, you can take differences in successive heights and see where they change sign to locate modes.)

Here is some output from density, including the first 100 heights, rounded to three places here.

DEN = density(x, adj=.5)
DEN

Call:
density.default(x = x, adjust = 0.5)

Data: x (950 obs.);     Bandwidth 'bw' = 6.928

x                y
Min.   :-13.61   Min.   :1.154e-06
1st Qu.: 54.26   1st Qu.:1.829e-03
Median :122.13   Median :3.131e-03
Mean   :122.13   Mean   :3.680e-03
3rd Qu.:190.00   3rd Qu.:5.655e-03
Max.   :257.87   Max.   :1.010e-02

round(DEN\$y[1:100],3)
[1] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[11] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[21] 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.001 0.001 0.001
[31] 0.001 0.001 0.001 0.001 0.002 0.002 0.002 0.002 0.002 0.003
[41] 0.003 0.003 0.004 0.004 0.004 0.005 0.005 0.006 0.006 0.006
[51] 0.007 0.007 0.007 0.008 0.008 0.009 0.009 0.009 0.009 0.010
[61] 0.010 0.010 0.010 0.010 0.010 0.010 0.010 0.010 0.010 0.010
[71] 0.009 0.009 0.009 0.009 0.008 0.008 0.008 0.007 0.007 0.007
[81] 0.006 0.006 0.006 0.005 0.005 0.005 0.004 0.004 0.004 0.004
[91] 0.003 0.003 0.003 0.003 0.003 0.002 0.002 0.002 0.002 0.002

• KDE's seems a very reasonable choice to my problem, anyway I would love to see/hear more ideas from SO community. Thank you BruceET. Now I just need to code my own KDE in Python (as far as I know there is no KDE in python. – Lin Sep 20 '19 at 5:32
• Wikipedia has a good article on KDE, note the Python bullet. Also, google 'density estimation Python', results of which seem promising. For very clear expositions read articles/books by Bernard Silverman. – BruceET Sep 20 '19 at 5:56

Not sure what kind of answer you're looking for, but thought I might give it a try. If you have a few of these, might be easiest to do it by hand. Cut off the sample at the trough of the left-most "bump" and find sample mode.

Since you're asking this, I'm assuming you have a bunch of these and couldn't do it by hand. Here's a pseudocode of the algorithm:

• Set i = 0
• Start with a window from 0 to i, find the mode of the data within that window.
• If the mode increases, continue increasing i and repeat the previous step.
• If the mode stop increasing, you've found your "minimum" mode.

For "maximum" mode, just do the same from the right.

• I don't think that will work, since It will stop in the first bump. For example from right to left it will stop around 240 since we have a small peak there. – Lin Sep 20 '19 at 2:50
• Ah good point. You'll want to smooth it out somewhat then. A "moving average" over its +/- 5 neighbors or something should work? – Art Sep 20 '19 at 3:18
• I cannor infer the window since the datasets can be very different. Smoothing is an alternative. But I think that there is a way to find the n-modal values, can't figureout how. – Lin Sep 20 '19 at 3:25