python实现常见统计量与统计分布

一. 基本概念

1. 总体

一个统计问题研究对象的全体称为总体，构成总体的每个成员称为个体。

2. 样本

一般由于总体的数量是非常巨大的，不可能对全部总体进行研究，通常会对总体进行抽样进行研究。

从总体中按一定规则抽出的一部分个体称为样本，样品的个数称为样本容量。

3. 统计量

设 $X_1,X_2{X}_3...X_n$ 为取自某总体的样本， 若样本函数中$T=(X_1,X_2{X}_3...X_n)$不含有任何未知参数， 则称$T$为统计量。

二. 常见统计量

1. 样本均值

设$X_1,X_2...X_n$为取自某总体的样本，样本均值$\overline{X}=\frac{1}{n}{\sum }_{i=1}^n{X}_i$。

import numpy as np

if __name__ == '__main__':
    data = np.array([3, 4, 5, 6, 2, 1])
    print(data.mean())

2. 样本方差与标准差

设$X_1,X_2...X_n$为取自某总体的样本，$\overline{X}$为样本均方差，则样本方差$S^2=\frac{1}{n-1}{\sum }_{i=1}^n(X_i-\overline{X}{)}^2$，$S$称为样本标准差。

import numpy as np
from scipy import stats

if __name__ == '__main__':
    data = np.array([3, 4, 5, 6, 2, 1])
    # ddof为0时自由度是n, ddof参数为1时自由度是n-1,所以我们需要指定该参数
    print(np.var(data, ddof=1))
    print(np.std(data, ddof=1))

3. 偏度

设$X_1,X_2...X_n$为取自某总体的样本，则偏度$SK=\frac{n}{(n-1)(n-2)}\ast \frac{{\sum }_{i=1}^n(X_i-\overline{X}{)}^3}{S^3}$其中$n$ 为样本量，$S$为样本标准差。

1. 偏度的计算

   import pandas as pd
   import numpy as np
   from scipy import stats
   
   
   def my_skew(data):
       n = np.size(data)
       mean = np.mean(data)
       std = np.std(data, ddof=1)
   
       tmp1 = np.sum(np.power(data - mean, 3))
       return (n * tmp1) / ((n - 1) * (n - 2) * (std ** 3))
   
   
   if __name__ == '__main__':
       # 1.使用scipy计算偏度,
       data = np.array([25.78, 15.05, 4.26, 19.14, 3.30, -35.75, 25.62, 15.15, -0.72, 17.25])
       print("scipy skewness:", stats.skew(data, bias=False))
   
       # 2.使用pandas 计算
       df = pd.DataFrame()
       df['data'] = data
       print('pandas skewness: ', df['data'].skew())
   
       # 3.自己实现偏度计算
       print("my skewness: ", my_skew(data))
   

2. 偏度与数据分布的关系

   • 如果一组数据是对称的，则偏度为0；如果偏度大于0，则数据右偏；如果偏度小于0，则数据左偏。
   • 一般偏度的绝对值大于1，数据为高度偏态分布；若绝对值在0.5 ~ 1之间，数据为中等偏态分布；绝对值越接近0，数据的偏斜程度越小。

   
   from scipy import stats
   import matplotlib.pyplot as plt
   
   plt.rcParams["font.family"] = "SimHei"  # 设置字体
   plt.rcParams["axes.unicode_minus"] = False  # 正常显示负号
   
   
   if __name__ == '__main__':
       fig, ax = plt.subplots(3, 1, figsize=(10, 8))
   
       data1 = stats.skewnorm(4).rvs(size = 10000)
       skew1 = stats.skew(data1, bias=False)
   
       data2 = stats.norm().rvs(size = 10000)
       skew2 = stats.skew(data2, bias=False)
   
       data3 = stats.skewnorm(-4).rvs(size = 10000)
       skew3 = stats.skew(data3, bias=False)
   
       ax[0].hist(data1, bins=50, label='右偏分布, skewness:{:.4f}'.format(skew1))
       ax[1].hist(data2, bins=50, label='无偏分布, skewness:{:.4f}'.format(skew2))
       ax[2].hist(data3, bins=50, label='左偏分布, skewness:{:.4f}'.format(skew3))
   
       for i in range(3):
           ax[i].legend()
   
       plt.show()
   
   

   运行结果：
   

4. 峰度

设$X_1,X_2...X_n$为取自某总体的样本，峰度$K=\frac{n(n+1){\sum }_{i=1}^n(X_i-\overline{X}{)}^4-3[{\sum }_{i=1}^n(X_i-\overline{X}{)}^2{]}^2(n-1)}{(n-1)(n-2)(n-3)S^4}$，其中$n$为样本量，$S$为样本标准差，峰度用来对数据分布的平峰或尖峰程度的测量。

1. 峰度的计算

   import pandas as pd
   import numpy as np
   from scipy import stats
   
   def my_kurtosis(data):
       n = np.size(data)
       mean = np.mean(data)
       std = np.std(data, ddof=1)
   
       tmp1 = n * (n + 1) * np.sum(np.power(data - mean, 4)) - 3 * (n - 1) * np.power(np.sum(np.power(data - mean, 2)), 2)
       tmp2 = (n - 1) * (n - 2) * (n -3) * np.power(std, 4)
       return tmp1 / tmp2
   
   
   if __name__ == '__main__':
       # 1.使用scipy计算峰度,
       data = np.array([25.78, 15.05, 4.26, 19.14, 3.30, -35.75, 25.62, 15.15, -0.72, 17.25])
       print("scipy kurtosis:", stats.kurtosis(data, bias=False))
   
       # 2.使用pandas 计算
       df = pd.DataFrame()
       df['data'] = data
       print('pandas kurtosis: ', df['data'].kurt())
   
       # 3.自己实现峰度计算
       print("my kurtosis: ", my_kurtosis(data))
   

2. 峰度与数据分布的关系

   • 峰态通常是与正太分布相比较而言，正太分布的峰度为0，峰度大于0表明数据呈现尖峰分布，峰度小于0表明数据呈现平峰分布。
     
   • 通过对不同分布进行1000000次抽样，并计算抽样值的峰度，再对比抽样值的频率分布

   from scipy import stats
   from scipy.stats import norm
   from scipy.stats import laplace
   from scipy.stats import cosine
   import matplotlib.pyplot as plt
   
   plt.rcParams["font.family"] = "SimHei"  # 设置字体
   plt.rcParams["axes.unicode_minus"] = False  # 正常显示负号
   
   
   if __name__ == '__main__':
   
       fig, ax = plt.subplots(3, 1, figsize= (10, 8))
       sample_size = 1000000
       # 从标准正太分布抽样
       data1 = norm(loc = 0, scale = 1).rvs(size = sample_size)
       kur1 = stats.kurtosis(data1, bias = False)
   
       # 从拉普拉斯分布抽样
       data2 = laplace().rvs(size = sample_size)
       kur2 = stats.kurtosis(data2, bias =False)
   
       # 从cosine分布抽样
       data3 = cosine().rvs(size = sample_size)
       kur3 = stats.kurtosis(data3, bias=False)
   
       ax[0].hist(data1, bins = 100, density=True, label = '标准正太分布, kurtosis:{:.2f}'.format(kur1))
       ax[0].legend()
   
       ax[1].hist(data2, bins = 100, density=True, label = '尖峰分布, kurtosis:{:.2f}'.format(kur2))
       ax[1].legend()
   
       ax[2].hist(data3, bins = 100, density=True, label = '平峰分布, kurtosis:{:.2f}'.format(kur3))
       ax[2].legend()
   
       plt.show()
   

   运行结果：
   

三. 常见统计分布

正太分布是统计中最常见的分布，以标准正太分布为基础构造的三个重要的统计量在实际中也有广泛的应用。

1. 卡方分布

随机变量$X_1,X_2...X_n$相互独立，且都服从标准正太分布，则随机变量$Y=X_1^2+...+X_n^2$的分布称为自由度为n的卡方分布，记为$Y\sim {\chi }^2(n)$

• 期望：$E(Y)=n$

• 方差：$D(Y)=2n$

• 画出不同参数下的卡方分布，自由度分别为$5,10,20$

  
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.stats  import chi2
  
  plt.rcParams["font.family"] = "SimHei"  # 设置字体
  plt.rcParams["axes.unicode_minus"] = False  # 正常显示负号
  
  if __name__ == '__main__':
  
      fig,ax = plt.subplots(figsize=(10, 8))
      params = [(5, 'red'), (10, 'blue'), (15, 'green')]
      x = np.linspace(0, 30, 100)
      for param in params:
          n = param[0]
          color = param[1]
          mean, var = chi2.stats(n, moments='mv')
          y = chi2(n).pdf(x)
          ax.plot(x, y, color = color, label= 'n = {:.2f}, 均值:{:.2f}, 方差: {:.4f}'.format(n, mean, var))
  
      ax.legend()
      plt.show()
  
  
  

  运行结果：
  

2. t分布

随机变量$X_1\sim N(0,1)$，$X_2\sim {\chi }^2(n)$，$X_1，X_2$相互独立，则$Y=\frac{X_1}{\sqrt{X_2/n}}$的分布称为自由度为$n$的$t$分布，记为$Y\sim t(n)$

• 期望：$E(Y)=0$

• 方差：$D(Y)=\frac{n}{n-2}(n>2)$

• 画出不同参数下的$t$分布，自由度分别为$3,6,15$

  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.stats  import t
  
  plt.rcParams["font.family"] = "SimHei"  # 设置字体
  plt.rcParams["axes.unicode_minus"] = False  # 正常显示负号
  
  if __name__ == '__main__':
  
      fig,ax = plt.subplots(figsize=(10, 8))
      params = [(3, 'red'), (6, 'blue'), (15, 'green')]
      x = np.linspace(-10, 10, 1000)
      for param in params:
          n = param[0]
          color = param[1]
          mean, var = t.stats(n, moments='mv')
          y = t(n).pdf(x)
          ax.plot(x, y, color = color, label= 'n = {:.2f}, 均值:{:.2f}, 方差: {:.4f}'.format(n, mean, var))
  
      ax.legend()
      ax.grid()
      plt.show()
  

  运行结果：
  

3. F分布

随机变量$X_1，X_2$相互独立，$X_1\sim {\chi }^2(m)，X_2\sim {\chi }^2(n)$，则$Y=\frac{X_1/m}{X_2/n}$的分布称为自由度为$m$与$n$的$F$分布，记为$Y\sim F(m,n)$

• 期望：$E(Y)=\frac{n}{n-2}(n>2)$

• 方差：$E(Y)=\frac{2n^2(m+n-2)}{m(n-2{)}^2(n-4)}(n>4)$

• 画出不同参数下的$F$分布，自由度分别为$(100,20),(100,40),(100,60)$

  
  import numpy as np
  import matplotlib.pyplot as plt
  from scipy.stats  import f
  
  plt.rcParams["font.family"] = "SimHei"  # 设置字体
  plt.rcParams["axes.unicode_minus"] = False  # 正常显示负号
  
  if __name__ == '__main__':
  
      fig,ax = plt.subplots(figsize=(10, 8))
      params = [(100,20,'red'), (100,40,'blue'), (100,60,'green')]
      x = np.linspace(0, 10, 1000)
      print(x)
      for param in params:
          m = param[0]
          n = param[1]
          color = param[2]
          mean, var = f.stats(m, n, moments='mv')
          y = f(m,n).pdf(x)
          ax.plot(x, y, color = color, label= 'm = {:.2f}, n = {:.2f}, 均值:{:.2f}, 方差: {:.4f}'.format(m, n, mean, var))
  
      ax.legend()
      ax.grid()
      plt.show()
  
  

  运行结果：