随机变量

有一些通用的概率分布类被封装在continuous random variables以及
discrete random variables中。有80多个连续性随机变量(RVs)以及10余个离散随机变量已经用
这些类建立。同样，新的程序和分布可以被用户新建（如果你构造了一个，请提供它给我们帮助发展
这个包）。

所有统计函数被放在子包scipy.stats中，且有这些函数的一个几乎完整的列表可以使用
info(stats)获得。这个列表里的随机变量也可以从stats子包的docstring中获得介绍。

在接下来的讨论中，我们着重于连续性随机变量(RVs)。几乎所有离散变量也符合下面的讨论，
尽管我们将“离散分布的特殊之处”指出它们的一些区别。

下面的示例代码我们假设scipy.stats包已被下述方式导入。

  1. >>> from scipy import stats

有些例子假设对象被这样的方式导入（不用输完整路径）了。

  1. >>> from scipy.stats import norm

获得帮助

所有分布可以使用help函数得到解释。为获得这些信息只需要使用像这样的简单调用：

  1. >>> print norm.__doc__

作为例子，我们用这种方式获取分布的上下界

  1. >>> print 'bounds of distribution lower: %s, upper: %s' % (norm.a,norm.b)
  2. bounds of distribution lower: -inf, upper: inf

我们可以通过调用dir(norm)来获得关于这个（正态）分布的所有方法和属性。应该看到，
一些方法是私有方法尽管其并没有以名称表示出来（比如它们前面没有以下划线开头），
比如veccdf就只用于内部计算（试图使用那些方法将引发警告，因为它们可能会在后续开发中被移除）

为了获得真正的主要方法，我们列举冻结分布的方法（我们将在下文解释何谓冻结）

  1. >>> rv = norm()
  2. >>> dir(rv)  # reformatted
  3.     ['__class__', '__delattr__', '__dict__', '__doc__', '__getattribute__',
  4.     '__hash__', '__init__', '__module__', '__new__', '__reduce__', '__reduce_ex__',
  5.     '__repr__', '__setattr__', '__str__', '__weakref__', 'args', 'cdf', 'dist',
  6.     'entropy', 'isf', 'kwds', 'moment', 'pdf', 'pmf', 'ppf', 'rvs', 'sf', 'stats']

最后，我们能通过内省获得所有的可用分布的信息。

  1. >>> import warnings
  2. >>> warnings.simplefilter('ignore', DeprecationWarning)
  3. >>> dist_continu = [d for d in dir(stats) if
  4. ...                 isinstance(getattr(stats,d), stats.rv_continuous)]
  5. >>> dist_discrete = [d for d in dir(stats) if
  6. ...                  isinstance(getattr(stats,d), stats.rv_discrete)]
  7. >>> print 'number of continuous distributions:', len(dist_continu)
  8. number of continuous distributions: 84
  9. >>> print 'number of discrete distributions:  ', len(dist_discrete)
  10. number of discrete distributions:   12

通用方法

连续随机变量的主要公共方法如下：

  • rvs:随机变量（就是从这个分布中抽一些样本）
  • pdf：概率密度函数。
  • cdf：累计分布函数
  • sf：残存函数（1-CDF）
  • ppf：分位点函数（CDF的逆）
  • isf：逆残存函数（sf的逆）
  • stats:返回均值，方差，（费舍尔）偏态，（费舍尔）峰度。
  • moment:分布的非中心矩。

让我们使用一个标准正态(normal)随机变量(RV)作为例子。

  1. >>> norm.cdf(0)
  2. 0.5

为了计算在一个点上的cdf，我们可以传递一个列表或一个numpy数组。

  1. >>> norm.cdf([-1., 0, 1])
  2. array([ 0.15865525,  0.5       ,  0.84134475])
  3. >>> import numpy as np
  4. >>> norm.cdf(np.array([-1., 0, 1]))
  5. array([ 0.15865525,  0.5       ,  0.84134475])

相应的，像pdf,cdf之类的简单方法可以用np.vectorize进行矢量化.

一些其他的实用通用方法:

  1. >>> norm.mean(), norm.std(), norm.var()
  2. (0.0, 1.0, 1.0)
  3. >>> norm.stats(moments = "mv")
  4. (array(0.0), array(1.0))

为了找到一个分布的中中位数，我们可以使用分位数函数ppf，它是cdf的逆。

  1. >>> norm.ppf(0.5)
  2. 0.0

为了产生一个随机变量列,使用size关键字参数。

  1. >>> norm.rvs(size=5)
  2. array([-0.35687759,  1.34347647, -0.11710531, -1.00725181, -0.51275702])

不要认为norm.rvs(5)产生了五个变量。

  1. >>> norm.rvs(5)
  2. 7.131624370075814

欲知其意，请看下一部分的内容。

偏移(Shifting)与缩放(Scaling)

所有连续分布可以操纵loc以及scale参数调整分布的location和scale属性。作为例子，
标准正态分布的location是均值而scale是标准差。

  1. >>> norm.stats(loc = 3, scale = 4, moments = "mv")
  2. (array(3.0), array(16.0))

通常经标准化的分布的随机变量X可以通过变换(X-loc)/scale获得。它们的默认值是loc=0以及scale=1.

聪明的使用loc与scale可以帮助以灵活的方式调整标准分布达到所想要的效果。
为了进一步说明缩放的效果，下面给出期望为1/λ指数分布的cdf。

F(x)=1−exp(−λx)

通过像上面那样使用scale，可以看到如何得到目标期望值。

  1. >>> from scipy.stats import expon
  2. >>> expon.mean(scale=3.)
  3. 3.0

均匀分布也是令人感兴趣的：

  1. >>> from scipy.stats import uniform
  2. >>> uniform.cdf([0, 1, 2, 3, 4, 5], loc = 1, scale = 4)
  3. array([ 0.  ,  0.  ,  0.25,  0.5 ,  0.75,  1.  ])

最后，联系起我们在前面段落中留下的norm.rvs(5)的问题。事实上，像这样调用一个分布，
其第一个参数，像之前的5，是把loc参数调到了5，让我们看：

  1. >>> np.mean(norm.rvs(5, size=500))
  2. 4.983550784784704

在这里，为解释最后一段的输出：norm.rvs(5)产生了一个正态分布变量，其期望，即loc=5.

我倾向于明确的使用loc,scale作为关键字而非像上面那样依赖参数的顺序。
因为这看上去有点令人困惑。我们在我们解释“冻结RV”的主题之前澄清这一点。

形态(shape)参数

虽然一般连续随机变量都可以通过赋予loc和scale参数进行偏移和缩放，但一些分布还需要
额外的形态参数确定其形态。作为例子，看到这个伽马分布，这是它的密度函数

γ(x,a)=λ(λx)a−1Γ(a)e−λx,

它要求一个形态参数a。注意到λ的设置可以通过设置scale关键字为1/λ进行。

让我们检查伽马分布的形态参数的名字的数量。（我们从上面知道其应该为1）

  1. >>> from scipy.stats import gamma
  2. >>> gamma.numargs
  3. 1
  4. >>> gamma.shapes
  5. 'a'

现在我们设置形态变量的值为1以变成指数分布。所以我们可以容易的比较是否得到了我们所期望的结果。

  1. >>>  gamma(1, scale=2.).stats(moments="mv")
  2. (array(2.0), array(4.0))

注意我们也可以以关键字的方式指定形态参数：

  1. >>> gamma(a=1, scale=2.).stats(moments="mv")
  2. (array(2.0), array(4.0))

冻结分布

不断地传递loc与scale关键字最终会让人厌烦。而冻结RV的概念被用来解决这个问题。

  1. >>> rv = gamma(1, scale=2.)

通过使用rv，在任何情况下我们不再需要包含scale与形态参数。显然，分布可以被多种方式使用，
我们可以通过传递所有分布参数给对方法的每次调用（像我们之前做的那样）或者可以对一个分
布对象先冻结参数。让我们看看是怎么回事：

  1. >>> rv.mean(), rv.std()
  2. (2.0, 2.0)

这正是我们应该得到的。

广播

像pdf这样的简单方法满足numpy的广播规则。作为例子，我们可以计算t分布的右尾分布的临界值
对于不同的概率值以及自由度。

  1. >>> stats.t.isf([0.1, 0.05, 0.01], [[10], [11]])
  2. array([[ 1.37218364,  1.81246112,  2.76376946],
  3.        [ 1.36343032,  1.79588482,  2.71807918]])

这里，第一行是以10自由度的临界值，而第二行是以11为自由度的临界值。所以，
广播规则与下面调用了两次isf产生的结果相同。

  1. >>> stats.t.isf([0.1, 0.05, 0.01], 10)
  2. array([ 1.37218364,  1.81246112,  2.76376946])
  3. >>> stats.t.isf([0.1, 0.05, 0.01], 11)
  4. array([ 1.36343032,  1.79588482,  2.71807918])

但是如果概率数组，如[0.1,0.05,0.01]与自由度数组,如[10,11,12]具有相同的数组形态，
则进行对应匹配，我们可以分别得到10%，5%，1%尾的临界值对于10，11,12的自由度。

  1. >>> stats.t.isf([0.1, 0.05, 0.01], [10, 11, 12])
  2. array([ 1.37218364,  1.79588482,  2.68099799])

离散分布的特殊之处

离散分布的方法的大多数与连续分布很类似。当然像pdf被更换为密度函数pmf，没有估计方法，
像fit就不能用了。而scale不是一个合法的关键字参数。Location参数，
关键字loc则仍然可以使用用于位移。

cdf的计算要求一些额外的关注。在连续分布的情况下，累积分布函数在大多数标准情况下是严格递增的，
所以有唯一的逆。而cdf在离散分布却一般是阶跃函数，所以cdf的逆，分位点函数，要求一个不同的定义：

ppf(q) = min{x : cdf(x) >= q, x integer}

为了更多信息可以看这里。

我们可以看这个超几何分布的例子

  1. >>> from scipy.stats import hypergeom
  2. >>> [M, n, N] = [20, 7, 12]

如果我们在一些整数点使用cdf，则它们的cdf值再作用ppf会回到开始的值。

  1. >>> x = np.arange(4)*2
  2. >>> x
  3. array([0, 2, 4, 6])
  4. >>> prb = hypergeom.cdf(x, M, n, N)
  5. >>> prb
  6. array([ 0.0001031991744066,  0.0521155830753351,  0.6083591331269301,
  7.         0.9897832817337386])
  8. >>> hypergeom.ppf(prb, M, n, N)
  9. array([ 0.,  2.,  4.,  6.])

如果我们使用的值不是cdf的函数值，则我们得到一个更高的值。

  1. >>> hypergeom.ppf(prb + 1e-8, M, n, N)
  2. array([ 1.,  3.,  5.,  7.])
  3. >>> hypergeom.ppf(prb - 1e-8, M, n, N)
  4. array([ 0.,  2.,  4.,  6.])

分布拟合

非冻结分布的参数估计的主要方法：

  • fit：分布参数的极大似然估计，包括location与scale
  • fit_loc_scale: 给定形态参数确定下的location和scale参数的估计
  • nnlf:负对数似然函数
  • expect:计算函数pdf或pmf的期望值。

性能问题与注意事项

分布方法的性能与运行速度根据分布的不同表现差异极大。方法的结果可以由两种方式获得，
精确计算或使用独立于各具体分布的通用算法。

精确计算一般更快。为了进行精确计算，要么直接使用解析公式，要么使用scipy.special中的
函数，对于rvs还可以使用numpy.random里的函数。

另一方面，如果不能进行精确计算，将使用通用方法进行计算。于是为了定义一个分布，
只有pdf异或cdf是必须的；通用方法使用数值积分和求根法进行求解。作为例子，
rgh = stats.gausshyper.rvs(0.5, 2, 2, 2, size=100)以这种方式创建了100个随机变量
（抽了100个值），这在我的电脑上花了19秒（译者：我花了3.5秒），
对比取一百万个标准正态分布的值只需要1秒。

遗留问题

scipy.stats里的分布最近进行了升级并且被仔细的检查过了，不过仍有一些问题存在。

  • 分布在很多参数区间上的值被测试过了，然而在一些奇葩的临界条件，仍然可能有错误的值存在。
  • fit的极大似然估计以默认值作为初始参数将不会工作的很好，用户必须指派合适的初始参数。
    并且，对于一些分布使用极大似然估计本身就不是一个好的选择。