Question

1. scipy中的stats模块中包含了很多概率分布的随机变量

   • 所有的连续随机变量都是rv_continuous的派生类的对象
   • 所有的离散随机变量都是rv_discrete的派生类的对象

1. 连续随机变量

1. 查看所有的连续随机变量：

   
   xxxxxxxxxx
   [k for k,v in stats.__dict__.items() if isinstance(v,stats.rv_continuous)]

   

2. 连续随机变量对象都有如下方法：

   • rvs(*args, **kwds)：获取该分布的一个或者一组随机值
   • pdf(x, *args, **kwds)：概率密度函数在x处的取值
   • logpdf(x, *args, **kwds) ：概率密度函数在x处的对数值
   • cdf(x, *args, **kwds)：累积分布函数在x处的取值
   • logcdf(x, *args, **kwds)：累积分布函数在x处的对数值
   • sf(x, *args, **kwds) ：生存函数在x处的取值，它等于1-cdf(x)
   • logsf(x, *args, **kwds) ：生存函数在x处的对数值
   • ppf(q, *args, **kwds)：累积分布函数的反函数
   • isf(q, *args, **kwds) ：生存函数的反函数
   • moment(n, *args, **kwds) n-th order non-central moment of distribution.
   • stats(*args, **kwds)：计算随机变量的期望值和方差值等统计量
   • entropy(*args, **kwds) ：随机变量的微分熵
   • expect([func, args, loc, scale, lb, ub, ...])：计算 $ MathJax-Element-72 $ 的期望值
   • median(*args, **kwds)：计算该分布的中值
   • mean(*args, **kwds)：计算该分布的均值
   • std(*args, **kwds)：计算该分布的标准差
   • var(*args, **kwds)：计算该分布的方差
   • interval(alpha, *args, **kwds)：Confidence interval with equal areas around the median.
   • __call__(*args, **kwds)：产生一个参数冻结的随机变量
   • fit(data, *args, **kwds) ：对一组随机取样进行拟合，找出最适合取样数据的概率密度函数的系数
   • fit_loc_scale(data, *args)：Estimate loc and scale parameters from data using 1st and 2nd moments.
   • nnlf(theta, x)：返回负的似然函数

   其中的args/kwds参数可能为（具体函数具体分析）：

   • arg1, arg2, arg3,...: array_like.The shape parameter(s) for the distribution
   • loc : array_like.location parameter (default=0)
   • scale : array_like.scale parameter (default=1)
   • size : int or tuple of ints.Defining number of random variates (default is 1).
   • random_state : None or int or np.random.RandomState instance。If int or RandomState, use it for drawing the random variates. If None, rely on self.random_state. Default is None.

   

3. 这些连续随机变量可以像函数一样调用，通过loc和scale参数可以指定随机变量的偏移和缩放系数。

   • 对于正态分布的随机变量而言，这就是期望值和标准差

2. 离散随机变量

1. 查看所有的连续随机变量：

   
   xxxxxxxxxx
     [k for k,v in stats.__dict__.items() if isinstance(v,stats.rv_discrete)]

   

2. 离散随机变量对象都有如下方法：

   • rvs(<shape(s)>, loc=0, size=1)：生成随机值
   • pmf(x, <shape(s)>, loc=0)：概率密度函数在x处的值
   • logpmf(x, <shape(s)>, loc=0)：概率密度函数在x处的对数值
   • cdf(x, <shape(s)>, loc=0)：累积分布函数在x处的取值
   • logcdf(x, <shape(s)>, loc=0) ：累积分布函数在x处的对数值
   • sf(x, <shape(s)>, loc=0) ：生存函数在x处的值
   • logsf(x, <shape(s)>, loc=0, scale=1)：生存函数在x处的对数值
   • ppf(q, <shape(s)>, loc=0) ：累积分布函数的反函数
   • isf(q, <shape(s)>, loc=0) ：生存函数的反函数
   • moment(n, <shape(s)>, loc=0)：non-central n-th moment of the distribution. May not work for array arguments.
   • stats(<shape(s)>, loc=0, moments='mv') ：计算期望方差等统计量
   • entropy(<shape(s)>, loc=0) ：计算熵
   • expect(func=None, args=(), loc=0, lb=None, ub=None, conditional=False)：Expected value of a function with respect to the distribution. Additional kwd arguments passed to integrate.quad
   • median(<shape(s)>, loc=0)：计算该分布的中值
   • mean(<shape(s)>, loc=0)：计算该分布的均值
   • std(<shape(s)>, loc=0)：计算该分布的标准差
   • var(<shape(s)>, loc=0)：计算该分布的方差
   • interval(alpha, <shape(s)>, loc=0) Interval that with alpha percent probability contains a random realization of this distribution.
   • __call__(<shape(s)>, loc=0)：产生一个参数冻结的随机变量

3. 我们也可以通过rv_discrete类自定义离散概率分布：

   
   xxxxxxxxxx
     x=range(1,7)
     p=(0.1,0.3,0.1,0.3,0.1,0.1)
     stats.rv_discrete(values=(x,p))

   只需要传入一个values关键字参数，它是一个元组。元组的第一个成员给出了随机变量的取值集合，第二个成员给出了随机变量每一个取值的概率

3. 核密度估计

1. 通常我们可以用直方图来观察随机变量的概率密度。但是直方图有个缺点：你选取的直方图区间宽度不同，直方图的形状也发生变化。核密度估计就能很好地解决这一问题。

2. 核密度估计的思想是：对于随机变量的每一个采样点 $ MathJax-Element-73 $ ，我们认为它代表一个以该点为均值、 $ MathJax-Element-78 $ 为方差的一个正态分布的密度函数 $ MathJax-Element-75 $ 。将所有这些采样点代表的密度函数叠加并归一化，则得到了核密度估计的一个概率密度函数：

   $ \frac 1N \sum_{i=1}^{N}f_i(x_i;s) $

   其中：

   • 归一化操作就是 $ MathJax-Element-76 $ ，因为每个点代表的密度函数的积分都是 1
   • $ MathJax-Element-78 $ 就是带宽参数，它代表了每个正态分布的形状

   如果采用其他的分布而不是正态分布，则得到了其他分布的核密度估计。

3. 核密度估计的原理是：如果某个样本点出现了，则它发生的概率就很高，同时跟他接近的样本点发生的概率也比较高。

4. 正态核密度估计：

   
   xxxxxxxxxx
   class scipy.stats.gaussian_kde(dataset, bw_method=None)

   参数：

   • dataset：被估计的数据集。

   • bw_method：用于设定带宽 $ MathJax-Element-78 $ 。可以为：

     • 字符串：如'scott'/'silverman'。默认为'scott'
     • 一个标量值。此时带宽是个常数
     • 一个可调用对象。该可调用对象的参数是gaussian_kde，返回一个标量值

   属性：

   • dataset：被估计的数据集
   • d：数据集的维度
   • n：数据点的个数
   • factor：带宽
   • covariance：数据集的相关矩阵

   方法：

   • evaluate(points)：估计样本点的概率密度
   • __call__(points)：估计样本点的概率密度
   • pdf(x)：估计样本的概率密度

   

5. 带宽系数对核密度估计的影响：当带宽系数越大，核密度估计曲线越平滑。

4. 常见分布

1. 二项分布：假设试验只有两种结果：成功的概率为 $ MathJax-Element-92 $ ，失败的概率为 $ MathJax-Element-80 $ 。 则二项分布描述了独立重复地进行 $ MathJax-Element-91 $ 次试验中，成功 $ MathJax-Element-108 $ 次的概率。

   • 概率质量函数：

     $ f(k;n,p)=\frac{n!}{k!(n-k)!}p^{k}(1-p)^{n-k} $

   • 期望： $ MathJax-Element-93 $

   • 方差： $ MathJax-Element-84 $

   scipy.stats.binom使用n参数指定 $ MathJax-Element-91 $ ；p参数指定 $ MathJax-Element-92 $ ；loc参数指定平移

   

2. 泊松分布：泊松分布使用 $ MathJax-Element-95 $ 描述单位时间（或者单位面积）中随机事件发生的平均次数（只知道平均次数，具体次数是个随机变量）。

   • 概率质量函数：

     $ f(k;\lambda)=\frac{e^{-\lambda}\lambda^{k}}{k!} $

     其物理意义是：单位时间内事件发生 $ MathJax-Element-108 $ 次的概率

   • 期望： $ MathJax-Element-95 $

   • 方差： $ MathJax-Element-95 $

   在二项分布中，如果 $ MathJax-Element-91 $ 很大，而 $ MathJax-Element-92 $ 很小。乘积 $ MathJax-Element-93 $ 可以视作 $ MathJax-Element-95 $ ，此时二项分布近似于泊松分布。

   泊松分布用于描述单位时间内随机事件发生的次数的分布的情况。

   scipy.stats.poisson使用mu参数来给定 $ MathJax-Element-95 $ ，使用 loc 参数来平移。

   

3. 用均匀分布模拟泊松分布：

   
   xxxxxxxxxx
     def make_poisson(lmd,tm):
       '''
       用均匀分布模拟泊松分布。 lmd为 lambda 参数； tm 为时间
       '''
       t=np.random.uniform(0,tm,size=lmd*tm) # 获取 lmd*tm 个事件发生的时刻
       count,tm_edges=np.histogram(t,bins=tm,range=(0,tm))#获取每个单位时间内，事件发生的次数
       max_k= lmd *2 # 要统计的最大次数
       dist,count_edges=np.histogram(count,bins=max_k,range=(0,max_k),density=True)
       x=count_edges[:-1]
       return x,dist,stats.poisson.pmf(x,lmd)  

   该函数首先随机性给出了 lmd*tm个事件发生的时间（时间位于区间[0,tm]）内。然后统计每个单位时间区间内，事件发生的次数。然后统计这些次数出现的频率。最后将这个频率与理论上的泊松分布的概率质量函数比较。

   

4. 指数分布：若事件服从泊松分布，则该事件前后两次发生的时间间隔服从指数分布。由于时间间隔是个浮点数，因此指数分布是连续分布。

   • 概率密度函数： $ MathJax-Element-96 $ ， $ MathJax-Element-104 $ 为时间间隔
   • 期望： $ MathJax-Element-98 $
   • 方差： $ MathJax-Element-99 $

   在scipy.stats.expon中，scale参数为 $ MathJax-Element-107 $ ；而loc用于对函数平移

   

5. 用均匀分布模拟指数分布：

   
   xxxxxxxxxx
    def make_expon(lmd,tm):
       '''
       用均匀分布模拟指数分布。 lmd为 lambda 参数； tm 为时间 
       '''
       t=np.random.uniform(0,tm,size=lmd*tm) # 获取 lmd*tm 个事件发生的时刻
       sorted_t=np.sort(t) #时刻升序排列
       delt_t=sorted_t[1:]-sorted_t[:-1] #间隔序列
       dist,edges=np.histogram(delt_t,bins="auto",density=True)
       x=edges[:-1]
       return x,dist,stats.expon.pdf(x,loc=0,scale=1/lmd) #scale 为 1/lambda

   

6. 伽玛分布：若事件服从泊松分布，则事件第 $ MathJax-Element-101 $ 次发生和第 $ MathJax-Element-102 $ 次发生的时间间隔为伽玛分布。由于时间间隔是个浮点数，因此指数分布是连续分布。

   • 概率密度函数： $ MathJax-Element-103 $ ， $ MathJax-Element-104 $ 为时间间隔
   • 期望： $ MathJax-Element-105 $
   • 方差： $ MathJax-Element-106 $

   在scipy.stats.gamma中，scale参数为 $ MathJax-Element-107 $ ；而loc用于对函数平移，参数a指定了 $ MathJax-Element-108 $

7. 用均匀分布模拟伽玛分布：

   
   xxxxxxxxxx
     def make_gamma(lmd,tm,k):
       '''
       用均匀分布模拟伽玛分布。 lmd为 lambda 参数； tm 为时间；k 为 k 参数
       '''
       t=np.random.uniform(0,tm,size=lmd*tm) # 获取 lmd*tm 个事件发生的时刻
       sorted_t=np.sort(t) #时刻升序排列
       delt_t=sorted_t[k:]-sorted_t[:-k] #间隔序列
       dist,edges=np.histogram(delt_t,bins="auto",density=True)
       x=edges[:-1]
       return x,dist,stats.gamma.pdf(x,loc=0,scale=1/lmd,a=k) #scale 为 1/lambda,a 为 k