skew normal distribution in scipy

distribution, python, scipy, statistics

Solution

From the Wikipedia description,

from scipy import linspace
from scipy import pi,sqrt,exp
from scipy.special import erf

from pylab import plot,show

def pdf(x):
    return 1/sqrt(2*pi) * exp(-x**2/2)

def cdf(x):
    return (1 + erf(x/sqrt(2))) / 2

def skew(x,e=0,w=1,a=0):
    t = (x-e) / w
    return 2 / w * pdf(t) * cdf(a*t)
    # You can of course use the scipy.stats.norm versions
    # return 2 * norm.pdf(t) * norm.cdf(a*t)


n = 2**10

e = 1.0 # location
w = 2.0 # scale

x = linspace(-10,10,n) 

for a in range(-3,4):
    p = skew(x,e,w,a)
    plot(x,p)

show()

If you want to find the scale, location, and shape parameters from a dataset use `scipy.optimize.leastsq`, for example using `e=1.0`,`w=2.0` and `a=1.0`,

fzz = skew(x,e,w,a) + norm.rvs(0,0.04,size=n) # fuzzy data

def optm(l,x):
    return skew(x,l[0],l[1],l[2]) - fzz

print leastsq(optm,[0.5,0.5,0.5],(x,))

should give you something like,

(array([ 1.05206154,  1.96929465,  0.94590444]), 1)

Problem

Does anyone know how to plot a skew normal distribution with scipy? I supose that stats.norm class can be used but I just can't figure out how. Furthermore, how can I estimate the parameters describing the skew normal distribution of a unidimensional dataset?

Original source