Weibull distribution and the data in the same figure (with numpy and scipy) [closed]
Asked Answered
R

1

-1

Here is my data: https://www.dropbox.com/s/xx02015pbr484es/Book2.xlsx

This is my output: enter image description here

And this is the desired output: enter image description here

As you can see,I want "the data" and the "Weibull distribution" together, (in the same figure).

This is my code:

(loc, scale) = s.exponweib.fit_loc_scale(mydata, 0.5, 0.5)
print loc, scale
x = np.linspace(mydata.min(), mydata.max(), 1000)
plt.plot(mydata, weib(mydata, loc, scale))
plt.hist(mydata, mydata.max(), normed=True)
plt.show()
Russian answered 24/11, 2013 at 4:22 Comment(3)
yes and I had deleted the duplicateRussian
my apologies, i must have had an old cache. Can you possibly upload your data as a text file? I can't read the xlsx file.Metalinguistic
@Metalinguistic there is xlrd module for excel format python-excel.orgCellaret
R
2

First, I think you want to fix location but not scale. (So both scale and shape can change). Second, I think (not 100% sure) that you can't have 0 in your data for Weibull (unless you hardcode a Weibull class yourself), so I changed your 0 to a small value 1e-8.

>>> xdata=array([1e-8,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,5,4,4,5,5,6,6,6,6,7,7,8,8,8,9,9,10,11,12,13,13,14,14,13,17,14,15,17,18,18,19,22,23,22,23,24,26,28,32,33,32,31,33,34,37,36,40,40,41,44,41,44,45,47,52,53,51,52,52,53,55,56,59,61,62,65,63,68,69,80,71,71,72,71,69,70,70,71,72,73,75,74,74,75,76,74,79,77,77,77,84,92,88,79,81,81,83,84,88,87,84,84,85,85,85,94,95,91,89,90,87,89,89,90,93,92,93,96,95,98,99,100,99,100,98,94,89,87,86,85,85,84,85,83,83,84,83,81,85,83,83,81,84,93,91,78,79,80,80,80,80,80,78,79,78,79,80,78,78,78,78,79,77,77,77,78,80,82,83,82,80,82,82,83,87,82,82,80,80,79,77,77,77,77,75,75,73,71,73,73,70,72,69,70,70,78,81,69,68,68,68,65,64,66,65,64,62,62,62,62,67,65,61,61,59,58,59,59,59,59,59,59,59,59,59,59,59,58,56,55,52,50,50,48,48,47,46,46,45,44,44,43,43,43,41,41,41,46,47,40,39,39,38,37,37,38,36,35,35,35,35,36,35,33,33,32,31,31,31,29,29,28,28,28,28,30,30,30,28,27,26,25,23,22,23,22,21,20,19,19,18,18,18,17,17,17,14,14,13,13,14,13,12,12,11,11,10,10,9,9,9,8,8,8,8,7,7,7,7,7,7,6,6,6,6,6,6,6,6,6,6,5,5,5,5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2])
>>> stats.exponweib.fit(xdata, floc=0, f0=1)
(1, 0.87120924706137459, 0, 35.884593247790207)
>>> stats.weibull_min.fit(xdata, floc=0)
(0.87120924706137459, 0, 35.884593247790036)
>>> p0, p1, p2=stats.weibull_min.fit(xdata, floc=0)
>>> ydata=stats.weibull_min.pdf(linspace(0, 120, 100), p0, p1, p2)
>>> plt.hist(xdata, 25, normed=True)
>>> plt.plot(linspace(0, 120, 100), ydata, '-')

enter image description here

The fit is actually correct. It looks ugly but it is due to a large proportion of your data is smallish.

Finally, I actually suspect that your original data is already frequency data not raw data, is that the case? (Let's assume your data is not interval-censored, that will require quite a bit of hardcode)

>>> import itertools
>>> x2data=list(itertools.chain(*[[i,]*val for i, val in enumerate(xdata)]))
>>> p0, p1, p2=stats.weibull_min.fit(x2data, floc=0)
>>> y2data=stats.weibull_min.pdf(linspace(0, 500, 100), p0, p1, p2)
>>> plt.plot(linspace(0, 500, 100), y2data, '-')
[<matplotlib.lines.Line2D object at 0x0360B6B0>]
>>> r1,r2,r3=plt.hist(x2data, bins=60, normed=True)

enter image description here Now the result looks much more reasonable. Although it still does not appears to be very closely Weibull distributed. More like http://en.wikipedia.org/wiki/Shifted_Gompertz_distribution.

Update: yes, if you have 0 in your data, you will get this when you call fit methods (scipy 0.12.0):

Warning (from warnings module):
  File "C:\Python27\lib\site-packages\scipy\optimize\optimize.py", line 438
     and numpy.max(numpy.abs(fsim[0] - fsim[1:])) <= ftol):
RuntimeWarning: invalid value encountered in subtract
Rang answered 24/11, 2013 at 16:49 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.