# Unimodality Test for Discrete Distribution

I have several collections of discrete datasets of integer values, say, from 1 to 10, inclusive. I am interested in characterizing the various distributions in these datasets, and it is important for my purposes whether each distribution is shaped unimodally or multimodally. I am not interested in the explicit discrete mode in the sense that, for example, 8 is the most common value. Rather, I am interested in whether each distribution is shaped unimodally or multimodally. For example, what I would call a unimodal shape: And a bimodal shape: Obviously, the above two distributions exhibit very different shapes. And, as you can see, some of the datasets are very large, containing tens or hundreds of thousands of values.

The discrete nature of the datasets is somewhat problematic because the most common test for unimodality, the Hartigan-Hartigan Dip test, assumes a continuous distribution. I am not a statistician, but it appears as if this is a rather rigid assumption. I tested an R implementation of the Dip test on what appeared to be some (artificial) perfectly unimodal discrete data and detected an incredibly small p-value, suggesting that the data was actually multimodal.

The two other commonly mentioned tests for unimodality appear to be the Silverman test and the excess mass test. I know little about the latter, but this explanation of the former seemed to hint that the Silverman test applies to discrete data, although it was not said explicitly.

So, my questions:

2) What is/are the best statistical test(s) to use for my data? Is the Silverman test an appropriate choice?

3) As I am not a statistician, where, if possible, might I find an already working implementation of the above statistical test (ideally Python or R)?

• You seem to want to overlook fluctuations of a certain size when you visually characterize the "shapes," because both of these graphics are strongly multimodal: they both have (clear, significant) modes at $1$. It would therefore help to have some more guidance from you concerning what you actually mean by "shape" and how great a deviation from unimodality you are willing to overlook. – whuber Sep 13 '16 at 22:15
• @whuber You raise a good point. Even in my first example, the possible mode at 1 only represents about 1% of the data points. I admittedly lack a clear idea of what an appropriate cutoff is (perhaps that's a parameter of some statistical test?), but it seemed at intuitively to me that the second example was more clearly multimodally shaped than the first. Any extent to which I can quantify this behavior is helpful. If there were some multimodality index on which the second example had a different value than the first, that would be great. I really just want to differentiate these curves. – David Sep 14 '16 at 3:16
• The standard way to differentiate such data--they can scarcely be called "curves," because they are really just ten values that appear to be counts--is by analyzing the corresponding $2\times 10$ contingency table. With such large counts, any reasonable analysis will work well, such as a chi-squared test. – whuber Sep 14 '16 at 12:30
• @whuber Yes, you are correct. I have performed chi-squared tests, and I can substantiate the the distributions are at least different. For what it is worth, I can also somewhat differentiate these on the basis of mean, standard deviation/variance, skewness, and kurtosis. My hope was to describe the difference in some more detail (and inferentially rather than descriptively). The extent to which the datasets show some bump near the lower tail varies greatly, which I was hoping to somehow quantify. Maybe unimodality was the wrong term. Beyond what I have already done, what else might I do? – David Sep 14 '16 at 13:48

The implication of the question is that these datasets tabulate counts of values drawn independently from a discrete distribution defined on an ordered set of values such as $1,2,\ldots, 10.$ When that is the case, these counts have a multinomial distribution.

If by "mode" we mean a strict local maximum height in the graph (padding the left and right of the graph with zeros), or something like that, and if the counts are all relatively large (more than 5 or so ought to do), then an attractive method to assess the number of modes in the underlying distribution is with bootstrapping. The problem this solves is that the number of modes in the distribution might differ from the number of modes in the data. By reconstructing the experiment from the distribution defined by the data, we can see to what extent the number of modes might vary. This is "bootstrapping."

Carrying out the bootstrapping is easy: write a function to compute the number of modes in a graph and another one to repeatedly sample from the graph's data and apply that function to the sample. Tabulate its results. ExampleR code is below. When given a dataset like the second one in the question, it plots this chart of the bootstrapped mode frequencies: In 676 of 1000 bootstrap samples there were two modes; in 293 there were three; and in 31 there were four. This indicates the data are consistent with an underlying distribution with two or perhaps three modes. There is some possibility of four. The likelihood of more than four is tiny.

These results intuitively make sense, because in the dataset the frequencies of the values $8,9,10$ were close and relatively small. It is possible the true frequency of $9$ is less than those of either $8$ or $10,$ causing there to be modes at $1,8,$ and $10.$ The bootstrapping gives us a sense of how much variation in modes is likely based on the random variation implied by the assumed sampling scheme.

The results for the first set of data are always two modes. That is because the variation among counts in the thousands or tens of thousands is so small that it is extremely unlikely these data came from a distribution with any other modes besides the obvious ones at $1$ and $8.$

#
# Compute strict modes.
# Input consists of the counts in the data, in order, including any zeros.
#
n.modes <- function(x) {
n <- length(x)+1
i <- c(0, x) < c(x, 0)
sum(i[-n] & !i[-1])
}
#
# Bootstrap the mode count in a dataset.
#
n.modes.boot <- function(x, n.boot=1e3)
tabulate(apply(rmultinom(n.boot, sum(x), x), 2, n.modes), ceiling(length(x)/2+1))
#
# Plot the bootstrap results.
#
library(ggplot2)
n.modes.plot <- function(f) {
X <- data.frame(Frequency=f / sum(f))
X$Count <- factor(1:nrow(X)) X <- subset(X, Frequency > 0) ggplot(X, aes(Count, Frequency, fill=Count)) + geom_col(show.legend=FALSE) } # # Show some examples. # x <- c(70, 30,20,40,60,70,110,170,180,165) f <- n.modes.boot(x) print(n.modes.plot(f))  You could try sampling() from density() with a large bin width or using binsmooth and sampling from the function. Or just add runif(n,min=-.5,max=.5) to your data. Difference between max and min is the size of your bins. #library(diptest) uniformP<-vector() multimodeP<-vector() for(itr in 1:300){ spaled<-sample(c(1,2,3,4,5),100,replace = T) #hist(spaled,breaks = 100, plot = T) spaled2<-spaled+runif(length(spaled) ,min=-.5, max=.5) #hist(spaled2,breaks = 100, plot = T) uniformP<-append(dip.test(spaled2)$p.value,uniformP)

spaled3<-(append(sample(c(1,5),40,replace = T),spaled))
#hist(spaled3,breaks = 100, plot = T)
spaled4<-spaled3+runif(length(spaled4) ,min=-.5, max=.5)
#hist(spaled4,breaks = 100, plot = T)

multimodeP<-append(dip.test(spaled4)$p.value,multimodeP) } sum(uniformP>.1) sum(uniformP>.05) sum(uniformP>.01) sum(multimodeP>.1) sum(multimodeP>.05) sum(multimodeP>.01)   264  279  297  22  34  76 In other words at significance of .01, test rejects all but 1% of uniform distribution but misses 25% of multimodal. I conclusion it works but not super great. • Because these are discrete distributions, what meaning would this smoothing have and why would its results be appropriate for answering the question? – whuber Aug 27 '18 at 21:51 • @whuber then run Hartigan's dip.test? If there is enough data to start with runif will work. At least it stopped dip.test identifying every discreet distribution as extremely multimodal. But every time I use sample() and density() R crashes so runif it is. – ran8 Aug 28 '18 at 17:31 • AFAIK, Hartigan's Dip Test does not apply. His definition, "A distribution function$F$is unimodal with mode$m$if$F$is convex in$(-\infty,m]$and concave in$[m,\infty),\$" rules out all nontrivial discrete distributions. (J. A. Hartigan and P. M. Hartigan, The Dip Test of Unimodality, Ann. Stat. 13 No. 1 1985.) – whuber Aug 28 '18 at 17:49
• Continuous distributions can be binned. I think of the bins as touching rectangles rather than single spikes. I think HDS requires interval scale. – ran8 Aug 29 '18 at 22:43
• Certainly you can bin. But that operation forces the distribution to be discrete and destroys the convexity property required by Hartigan even to define a mode. Although it's possible the Dip test might still be adapted to such a circumstance, that would have to be established by a separate theoretical result or through simulation. You have no valid basis to assert it will work. – whuber Aug 30 '18 at 13:13