skew normal distribution in scipy

Question

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?

Answer 1

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)