Bimodality coefficient (BMC)

bmctest.py

# IPython log file

import sys
import numpy
numpy.set_printoptions(precision=2, suppress=True)
import scipy
import scipy.stats

print("Python version = {sys.version}".format(sys=sys))
print("NumPy version = {numpy.version.version}".format(numpy=numpy))
print("SciPy version = {scipy.version.version}".format(scipy=scipy))

SEED = 10
print("Random Seed = {SEED}".format(SEED=SEED))

numpy.random.seed(SEED)

def bmc(x):
    exc_kurt = scipy.stats.kurtosis(x)
    gamma = scipy.stats.skew(x)
    n = len(x)
    n_term = (3 * (n - 1) ** 2)/((n - 2) * (n - 3))
    return (gamma ** 2 + 1) / (exc_kurt + n_term)

n = 100
x = numpy.arange(7)
p1 = numpy.array([1, 1, 1, 1, 1, 1, 1]) / 7
p2 = numpy.array([2.5, 1, 0, 0, 0, 1, 2.5]) / 7
r1 = scipy.stats.rv_discrete(values=(x, p1))
r2 = scipy.stats.rv_discrete(values=(x, p2))
x1 = r1.rvs(size=n)
x2 = r2.rvs(size=n)
print("Sampling {n} random numbers:".format(n=n))
print("Probs = {prob}, BMC = {bmc:.2f}".format(prob=p1, bmc=bmc(x1))) 
print("Probs = {prob}, BMC = {bmc:.2f}".format(prob=p2, bmc=bmc(x2)))