# Untitled

a guest
Jun 19th, 2012
493
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. skew normal distribution in scipy
2. from scipy import linspace
3. from scipy import pi,sqrt,exp
4. from scipy.special import erf
5.
6. from pylab import plot,show
7.
8. def pdf(x):
9. return 1/sqrt(2*pi) * exp(-x**2/2)
10.
11. def cdf(x):
12. return (1 + erf(x/sqrt(2))) / 2
13.
14. def skew(x,e=0,w=1,a=0):
15. t = (x-e) / w
16. return 2 / w * pdf(t) * cdf(a*t)
17. # You can of course use the scipy.stats.norm versions
18. # return 2 * norm.pdf(t) * norm.cdf(a*t)
19.
20.
21. n = 2**10
22.
23. e = 1.0 # location
24. w = 2.0 # scale
25.
26. x = linspace(-10,10,n)
27.
28. for a in range(-3,4):
29. p = skew(x,e,w,a)
30. plot(x,p)
31.
32. show()
33.
34. fzz = skew(x,e,w,a) + norm.rvs(0,0.04,size=n) # fuzzy data
35.
36. def optm(l,x):
37. return skew(x,l[0],l[1],l[2]) - fzz
38.
39. print leastsq(optm,[0.5,0.5,0.5],(x,))
40.
41. (array([ 1.05206154, 1.96929465, 0.94590444]), 1)