scipy.stats包简单教程

by 王宇

==生成分布==
在stats包中有许多分布，比如norm、gamma、expon等。每个分布都有缺省的参数。如果要改变这些缺省参数，那么可以给这些分布带上参数：
d = stats.norm(loc=2, scale=5)
d = stats.expon(scale=5)
d = stats.gamma(3, loc=2, scale=5)
上面每一条语句都完整的设定了一个分布的所有参数。例如expon这个分布在教科书里只有一个参数/lambda，在scipy中用scale这个参数代替了/lambda，两个参数的含义并不一样。Expon用不上loc参数，在设定时就可以忽略。
gamma这个分布有三个参数，所以需要额外设定。通用的设定模式是：
rv = generic(<shape(s)>,loc=0,scale=1)

==loc、scale参数的缺省值==
在后面你调用相应的函数，比如d.pdf时，有一条语句叫做：
args, loc, scale = self._fix_loc_scale(args, loc, scale)
_fix_loc_scale来自rv_generic类，具体做法是：
if scale is None:
scale = 1.0
if loc is None:
loc = 0.0
所以，缺省的loc是0.0，缺省的scale是1.0。

==loc、scale参数的含义是==
Stats包让几乎所有分布都有了loc和scale这两个参数，那么，这些参数的含义是什么呢？具体来看，在计算pdf时，用的语句基本上是：
x = (x-loc)*1.0/scale
return self._pdf(x) / scale
self._pdf调用的是具体分布的_pdf，如果分布是expon，那么就调用expon_gen的_pdf，它的代码是：
return exp(-x)。看到这里，想必你明白了loc和scale的功能是什么了。

所以，我们可以知道的是，如果一个分布在loc=0，scale=1参数下的mean是mu，standard variance是std，它在loc=a和scale=b的参数下的mean是mu*b+a，标准差是std*b
例如：
>>> stats.expon.stats()
(array(1.0), array(1.0))
>>> stats.expon(loc=0, scale=5).stats()
(array(5.0), array(25.0))
>>> stats.expon(loc=1, scale=5).stats()
(array(6.0), array(25.0))

>>> stats.norm.stats()
(array(0.0), array(1.0))
>>> stats.norm(loc=1, scale=4).stats()
(array(1.0), array(16.0))

>>> stats.gamma(1).stats()
(array(1.0), array(1.0))
>>> stats.gamma(1, scale=5).stats()
(array(5.0), array(25.0))
>>> stats.gamma(1, loc=6, scale=5).stats()
(array(11.0), array(25.0))

==怎么根据教科书设定参数==
这就引出一个有趣的问题，已知一个分布在教科书里的参数，我们怎么设定它的shapes、loc和scale呢？我试验了一下指数分布：

如果 =5，则令scale=0.2，这时对于任意x，pdf(x)的值为5*exp(-x*5)。

==分布的各种成员函数==
有了分布后，我们可以做一些有趣的事情。
比如d.pdf(x)告诉我们x点的概率密度：
>>> stats.norm.pdf(3)
0.0044318484119380075

d.rvs(size=n)用d这个分布做n次采样：
>>> stats.norm(1,2).rvs(10)
array([ 2.80733657, 1.0873374 , 3.97014527, 0.47112426, 3.36529344,
2.92216297, -2.11625248, 0.26656479, 1.03135824, 1.2057151 ])

==分布的参数估计方法==
另外，估计参数的方法是什么？
虽然每个分布都有fit方法，例如：
>>> stats.norm.fit([1,2])
array([ 1.50001249, 0.49996873])
返回值就是loc和scale
但是，这个fit方法大部分时候是错的，例如：

from scipy import *
import scipy.optimize as OPT
import scipy.stats as stats
def f(params):
return -sum(log(stats.expon.pdf(yList, scale=params[0])))

print "Fit:"
scale = stats.expon.fit(yList)[1]
print scale
print f((scale,))
print "fmin:"
scale = OPT.fmin(f, (0.1,))[0]
print scale
print f((scale,))
print "std:"
scale = SP.std(yList0)
print scale

yList是数据，f函数计算了-log(yList)，所以，我们用最大似然估计，也就是最小化f。以上程序运行结果为：
Fit:
1.0
20.2683027126
fmin:
Optimization terminated successfully.
Current function value: -18698.946436
Iterations: 21
Function evaluations: 42
0.0047216796875
-18698.9464363
std:
0.0411461918185
-13204.7559705

显然fit函数给出的值就是胡扯。Std函数给出的值好一点，但也不对。我们需要用fmin函数。

参考文献：http://www.rexx.com/~dkuhlman/scipy_course_01.html