## Code to generate density plots and CDF plots for several members
## of the generalized normal family of probability distributions.
##
## Note that this plots one of (at least) two distributions refered to
## as "generalized normal."
##
## Reuires numpy, matplotlib, and the special functions model of scipy.
import matplotlib.pyplot as plt
import numpy as np
import scipy.special as sp

def dens(X, k):
    if k!=0: Y = -np.log(1-k*X)/k
    else: Y = X
    Y = np.exp(-Y**2/2)/np.sqrt(2*np.pi)
    return Y/(1-k*X)

def cdf(X, k):
    if k!=0: Y = -np.log(1-k*X)/k
    else: Y = X
    return sp.ndtr(Y)

w = 1.5

plt.clf()

colors = ['aqua', 'lime', 'deeppink', 'darkorange', 'blue']
K = [-1, -0.5, 0, 0.5, 1]

m = 8

F = []
for c,k in zip(colors, K):
    if k==0: a,b=-m,m
    elif k>0: a,b = -m,min(m, 1/float(k))
    else: a,b = max(-m,1/float(k))+1e-8,m
    X = np.arange(a, b, 0.01)
    Y = dens(X, k)
    f = plt.plot(X, Y, '-', color=c, lw=w)
    F.append(f)
    plt.hold(True)

s = ["$\\kappa=%s$" % k for k in K]

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

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

plt.clf()

F = []
for c,k in zip(colors, K):
    if k==0: a,b=-m,m
    elif k>0: a,b = -m,min(m, 1/float(k))
    else: a,b = max(-m,1/float(k))+1e-8,m
    X = np.arange(a, b, 0.01)
    Y = cdf(X, k)
    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.xlim(-4,4)

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