Skip to content

Instantly share code, notes, and snippets.

@glciampaglia

glciampaglia/bmctest.py

Created May 7, 2022 22:11
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save glciampaglia/0ab7e0a5613ec28089106facdd876ed7 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save glciampaglia/0ab7e0a5613ec28089106facdd876ed7 to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. Giovanni Luca Ciampaglia created this gist May 7, 2022.
    35 changes: 35 additions & 0 deletions bmctest.py
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,35 @@
    # 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)))