##
## Python code to generate density and CDF plots for the
## generalized normal distribution (requires matplotlib, numpy,
## and the special functions module of scipy).
##

import matplotlib.pyplot as plt
import numpy as np
import scipy.special as sp

def dens(X, a, b, mu):
    Y = np.abs(X-mu)/a
    Y = np.exp(-Y**b)
    Y = Y*b/(2*a)
    Y /= np.exp(sp.gammaln(1/float(b)))
    return Y

def cdf(X, a, b, mu):
    Y = (np.abs(X-mu)/a)**b
    Y = sp.gammainc(1/float(b),Y)
    Y = Y * np.sign(X-mu) / 2 + 0.5
    return Y

w = 1.5

X = np.arange(-3, 3, 0.01)

plt.clf()

colors = ['aqua', 'black', 'lime', 'deeppink', 'darkorange', 'blue']
B = [0.5, 1, 1.5, 2, 3, 8]

F = []
for c,b in zip(colors, B):
    Y = dens(X, 1, b, 0)
    f = plt.plot(X, Y, '-', color=c, lw=w)
    F.append(f)
    plt.hold(True)

s = ["$\\beta=%s$" % b for b in B]

b = plt.legend(tuple(F), tuple(s), 'upper left')
plt.ylabel("Density")
b.draw_frame(False)

plt.savefig("generalized_normal_densities.svg")
plt.savefig("generalized_normal_densities.png")

plt.clf()

F = []
for c,b in zip(colors, B):
    Y = cdf(X, 1, b, 0)
    f = plt.plot(X, Y, '-', color=c, lw=w)
    F.append(f)
    plt.hold(True)

b = plt.legend(tuple(F), tuple(s), 'upper left')
plt.ylabel("Cumulative probability")
b.draw_frame(False)
plt.ylim(0,1)

plt.savefig("generalized_normal_cdfs.svg")
plt.savefig("generalized_normal_cdfs.png")