Guest User

Untitled

a guest
Jun 19th, 2012
207
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)
RAW Paste Data