skewnorm MWE script

skewnorm-mwe.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# NOTE: This script was run with NumPy version 1.7.1 (and 1.15.3)
# NOTE: This script was run with SciPy version 0.12.1 (and 1.1.0)
# NOTE: SciPy reference: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skewnorm.html (retrieved 2018-11-29)
# NOTE: Wolfram reference: https://reference.wolfram.com/language/ref/SkewNormalDistribution.html (retrieved 2018-11-29)

# Import modules ...
import matplotlib
import matplotlib.pyplot
import numpy
import scipy
import scipy.stats

# Create axis (-5 to +15 in steps of 0.01) ...
axis = numpy.linspace(-5.0, 15.0, 2001)

# Define parameters and create skewed normal distribution three ways ...
mean =  0.0
skew = 10.0
sigma = 1.0
dist1 = 2.0 * scipy.stats.norm.pdf(axis, loc = mean, scale = sigma) * scipy.stats.norm.cdf(skew * axis, loc = mean, scale = sigma)
dist2 = scipy.stats.skewnorm.pdf(axis, a = skew, loc = mean, scale = sigma)
dist3 = scipy.stats.norm.pdf(axis, loc = mean, scale = sigma) * scipy.special.erfc(-skew * (axis - mean) / (numpy.sqrt(2.0) * sigma))

# Define parameters and create skewed normal distribution three ways ...
mean = 10.0
skew = 10.0
sigma = 1.0
dist4 = 2.0 * scipy.stats.norm.pdf(axis, loc = mean, scale = sigma) * scipy.stats.norm.cdf(skew * axis, loc = mean, scale = sigma)
dist5 = scipy.stats.skewnorm.pdf(axis, a = skew, loc = mean, scale = sigma)
dist6 = scipy.stats.norm.pdf(axis, loc = mean, scale = sigma) * scipy.special.erfc(-skew * (axis - mean) / (numpy.sqrt(2.0) * sigma))

# Save data ...
with open("skewnorm-mwe.csv", "wt") as fobj:
    for i in xrange(axis.size):
        fobj.write("{0:e},{1:e},{2:e},{3:e},{4:e},{5:e},{6:e}\n".format(axis[i], dist1[i], dist2[i], dist3[i], dist4[i], dist5[i], dist6[i]))

# Save plot ...
matplotlib.pyplot.rcParams.update({"font.size" : 8})
fig, ax = matplotlib.pyplot.subplots(figsize = (9, 6))
ax.plot(axis, dist1, label = "mu = 0 (SciPy; documentation)")
ax.plot(axis, dist2, label = "mu = 0 (SciPy; in-built)")
ax.plot(axis, dist3, ".", label = "mu = 0 (Wolfram)")
ax.plot(axis, dist4, label = "mu = 10 (SciPy; documentation)")
ax.plot(axis, dist5, label = "mu = 10 (SciPy; in-built)")
ax.plot(axis, dist6, ".", label = "mu = 10 (Wolfram)")
ax.grid()
ax.legend(loc = "upper center")
fig.tight_layout()
fig.savefig("skewnorm-mwe.png", dpi = 300, format = "png", bbox_inches = "tight")