Detecting if an 1-dimenisional distribution is Multimodal

Question

I'm writing up some C++ code for one of my Master's coursework. What I'm actually doing at the moment isn't on the syllabus, but I wish to implement it anyway as it will allow me to produce my own memetic computing algorithm.

Onwards!

I have four different optimisation problems coded up in C++:

• Michalewicz [-512, 512] - meant to be larger
• Schwefel [-5.12, 5.12]
• Rastigin [-5.12, 5.12]
• Sphere [-5.12, 5.12]

I know from using these that the only problem which isn't multimodal is Sphere, as it only has one global minima and maxima, but the other three problems are riddled with local minima and local maxima.

My question is, where n = the dimensionality of the problem, I'm trying to produce a statistical test which tells the computer if the current optimisation problem is multimodal by taking a 1D slice (line). Depending on the output the program then selects a pool of operators (single solution or population-based optimisers) to attempt to find the global minima.

The first thing I did was to develop a distribution (in 1D) using the bounds of the problem. When I output these distributions in Excel the charts I get are what I expected. Next, I found this following paper: http://www.i3.psychologie.uni-wuerzburg.de/fileadmin/06020300/user_upload/Janczyk/pfister_et_al_2013_bimodality.pdf

They use Skewness and Kurtosis to find something called a bimodality coefficient. Now, I'm not able to reproduce these results in C++ as I suspect Matlab has some trickery in terms of how they compute Kurtosis, so my results states that every problem except Michalewicz is multimodal (Which is not correct!!)

Here's some data (1D lines from the optimisation prolem):

Sphere

25.2004 24.2064 23.2324 22.2784 21.3444 20.4304 19.5364 18.6624 17.8084 16.9744 16.1604 15.3664 14.5924 13.8384 13.1044 12.3904 11.6964 11.0224 10.3684 9.7344 9.1204 8.5264 7.9524 7.3984 6.8644 6.3504 5.8564 5.3824 4.9284 4.4944 4.0804 3.6864 3.3124 2.9584 2.6244 2.3104 2.0164 1.7424 1.4884 1.2544 1.0404 0.8464 0.6724 0.5184 0.3844 0.2704 0.1764 0.1024 0.0484 0.0144 0.0004 0.0064 0.0324 0.0784 0.1444 0.2304 0.3364 0.4624 0.6084 0.7744 0.9604 1.1664 1.3924 1.6384 1.9044 2.1904 2.4964 2.8224 3.1684 3.5344 3.9204 4.3264 4.7524 5.1984 5.6644 6.1504 6.6564 7.1824 7.7284 8.2944 8.8804 9.4864 10.1124 10.7584 11.4244 12.1104 12.8164 13.5424 14.2884 15.0544 15.8404 16.6464 17.4724 18.3184 19.1844 20.0704 20.9764 21.9024 22.8484 23.8144 24.8004 25.8064 26.8324

Rastign

15.08054031 12.44916322 5.634081898 -2.300742733 -8.269711697 -10.34815318 -8.356811024 -3.688270535 1.463148257 5.084118795 6.111823406 4.702612168 1.955348819 -0.719244438 -2.263010207 -2.371550982 -1.486566652 -0.435316748 0.069031276 -0.193614067 -0.872664091 -1.291325524 -0.871825677 0.487491228 2.285754028 3.620821814 3.631064359 1.966078457 -0.950323075 -4.013409675 -5.872922186 -5.532649855 -2.84744159 1.319464265 5.376581008 7.628965207 6.996082011 3.515915078 -1.594914797 -6.375268028 -8.888950878 -8.021360826 -3.97149892 1.776674672 7.009470734 9.652879198 8.608486302 4.214193116 -1.86474389 -7.279189126 -9.920750167 -8.757458438 -4.243997667 1.859122451 7.184423204 9.692563786 8.468277233 4.060912571 -1.759810354 -6.725172969 -8.968320054 -7.740942689 -3.66493783 1.5668076 5.90143842 7.748018971 6.575454804 3.056073443 -1.280114189 -4.713219557 -6.031660538 -4.97181358 -2.23431941 0.89973012 3.160516381 3.819244754 2.930019015 1.199675732 -0.425655394 -1.243328891 -1.11077162 -0.450071111 0.047857592 -0.142109989 -1.038342913 -2.093758866 -2.468030134 -1.508280562 0.803566029 3.68449903 5.794346702 5.824284718 3.181593178 -1.558712727 -6.695139462 -9.990991888 -9.618692643 -5.06779544 2.407550082 10.07026421 14.68369443 13.85125391 7.166887347

Schwfel

422.9184972 422.9076247 422.8927192 422.8737881 422.8508422 422.8238953 422.7929647 422.7580711 422.7192389 422.6764964 422.629876 422.5794143 422.5251525 422.4671364 422.4054168 422.34005 422.2710974 422.1986267 422.1227115 422.0434322 421.9608762 421.8751381 421.7863207 421.6945352 421.5999018 421.5025505 421.4026219 421.3002676 421.1956518 421.0889518 420.9803594 420.8700826 420.7583466 420.6453961 420.5314975 420.4169416 420.3020464 420.1871613 420.072672 419.9590062 419.8466416 419.7361161 419.6280413 419.5231219 419.4221833 419.3262139 419.2364339 419.1544183 419.0823469 419.0236428 418.985719 418.9605731 418.908803 418.841556 418.7632086 418.6763193 418.582663 418.4835824 418.3801483 418.273245 418.1636217 418.0519244 417.9387187 417.8245048 417.7097296 417.5947946 417.4800631 417.3658653 417.2525023 417.14025 417.0293614 416.9200693 416.8125885 416.707117 416.6038377 416.5029199 416.4045202 416.3087835 416.2158439 416.1258254 416.0388428 415.9550021 415.8744013 415.7971304 415.7232727 415.6529044 415.5860955 415.52291 415.4634062 415.4076369 415.3556499 415.3074884 415.2631906 415.2227907 415.1863185 415.1538 415.1252573 415.1007092 415.0801707 415.0636539 415.0511674 415.0427172 415.0383061

I'm not entirely sure on what to attempt to far. I've had a look at techniques on Wikipedia, but I'm unsure on how I'm supposed to implement them. For example, Ashman's D algorithm uses two different means and standard deviations? I can't see how this would help me with my problem.

Thanks

Answer 1

Jonathan, I am not sure if I completely understand your problem so please correct me if I am wrong.

1. You want to construct a statistical test to see if the data in a simple random sample come from a particular population distribution;

2. The reason that you wish to do this is that you want to tell your program which optimization algorithm to find a maximum or minimum for an objective function;

3. The objective function is proportional to a probability mass function and the sample is drawn by generating random numbers and plugging these into the objective function.

Is this roughly correct? You want to use a random number generator to have a peek at the objective function before deciding which algorithm to use? The goal I suspect is to minimize computation time and find the global optimum of this objective function?

One way forward that I would suggest is to define what you could call a peak or a trough in terms of what proportion of values of this sample from the objective function are higher or lower than x % of the values in the sample and to compute the distance (e.g. the euclidean distance) between the coordinates of these highest and lowest values. The coordinates are the random numbers that you used to generate the value of the objective function. If these are more than [y] units away from each other then these are different peaks or troughs. Else not.

I don't think that you need a statistical test if I understand what you are trying to do here.

Answer 2

I was able to calculate if my problem is multi-modal or not by calculating the amount of slope sign changes. What I mean by this is, if the slope goes from a positive, to a negative then this indicates the slope has changed. If this occurs more than twice, then we have a local maxima / minima.

This is the definition of a slope sign change

X = x2 - x1
Y = y2 - y1
Slope = Y / X

Answer 3

Inspired by Jonathan Brown and this, i made some code with explanations for completeness.

The goal is to find how many modes there are in a distribution

# dataIn <- "Use some data from the original question"
densHere <- density(dataIn, kernel = "gaussian", n = 1024*10)
slope = vector("numeric")
for (i in 2:length(densHere$"x")) {

    X = densHere$x[i] - densHere$x[i - 1] 
    Y = densHere$y[i] - densHere$y[i - 1]

    slope = c(slope, Y / X)
}

numberOfModes <- length(which(diff(sign(slope)) == -2))  # Now find only the slopes from positive to negative

Alternatively use splines to do the trick:

densitySpline <- smooth.spline(densHere$x, densHere$y, all.knots = TRUE, spar = 0.8)
s.0 <- predict(densitySpline, densitySpline$x, deriv = 0)  # Fit the density
s.1 <- predict(densitySpline, densitySpline$x, deriv = 1)  # Get the derivative of the fit
numberOfModes <- which(diff(sign(s.1$y)) == -2)

EDIT:

I probably came late to the party. The slope approach is unfortunately very much dependent to the smoothing done by the kernel. You will have to toy a lot with the parameters of density function.

EDIT 2:

While researching for the problem, i came across Dip Test. In R the package is called diptest. My result contains multimodal distributions and it finds nearly almost all of them, while not suffering from the usual problem of large sample size.