5.2 stat.py

#!/usr/bin/python
#  -*- coding:utf-8 -*-

import numpy as np
from scipy import stats
import math
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm


def calc_statistics(x):
    n = x.shape[0]  # 样本个数

    # 手动计算
    m = 0
    m2 = 0
    m3 = 0
    m4 = 0
    for t in x:
        m += t
        m2 += t*t
        m3 += t**3
        m4 += t**4
    m /= n
    m2 /= n
    m3 /= n
    m4 /= n

    mu = m
    sigma = np.sqrt(m2 - mu*mu)
    skew = (m3 - 3*mu*m2 + 2*mu**3) / sigma**3
    kurtosis = (m4 - 4*mu*m3 + 6*mu*mu*m2 - 4*mu**3*mu + mu**4) / sigma**4 - 3
    print ('手动计算均值、标准差、偏度、峰度：', mu, sigma, skew, kurtosis)

    # 使用系统函数验证
    mu = np.mean(x, axis=0)
    sigma = np.std(x, axis=0)
    skew = stats.skew(x)
    kurtosis = stats.kurtosis(x)
    return mu, sigma, skew, kurtosis


if __name__ == '__main__':
    d = np.random.randn(100000)
    print (d)
    mu, sigma, skew, kurtosis = calc_statistics(d)
    print ('函数库计算均值、标准差、偏度、峰度：', mu, sigma, skew, kurtosis)
    # 一维直方图
    mpl.rcParams[u'font.sans-serif'] = 'SimHei'
    mpl.rcParams[u'axes.unicode_minus'] = False
    y1, x1, dummy = plt.hist(d, bins=50, normed=True, color='g', alpha=0.75)
    t = np.arange(x1.min(), x1.max(), 0.05)
    y = np.exp(-t**2 / 2) / math.sqrt(2*math.pi)
    plt.plot(t, y, 'r-', lw=2)
    plt.title(u'高斯分布，样本个数：%d' % d.shape[0])
    plt.grid(True)
    plt.show()

    d = np.random.randn(100000, 2)
    mu, sigma, skew, kurtosis = calc_statistics(d)
    print ('函数库计算均值、标准差、偏度、峰度：', mu, sigma, skew, kurtosis)

    # 二维图像
    N = 30
    density, edges = np.histogramdd(d, bins=[N, N])
    print ('样本总数：', np.sum(density))
    density /= density.max()
    x = y = np.arange(N)
    t = np.meshgrid(x, y)
    fig = plt.figure(facecolor='w')
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(t[0], t[1], density, c='r', s=15*density, marker='o', depthshade=True)
    ax.plot_surface(t[0], t[1], density, cmap=cm.Accent, rstride=2, cstride=2, alpha=0.9, lw=0.75)
    ax.set_xlabel(u'X')
    ax.set_ylabel(u'Y')
    ax.set_zlabel(u'Z')
    plt.title(u'二元高斯分布，样本个数：%d' % d.shape[0], fontsize=20)
    plt.tight_layout(0.1)
    plt.show()

输出：

[-0.9447748   2.05723148  1.97935487 ..., -0.10145314  0.55822039
  1.02833672]
手动计算均值、标准差、偏度、峰度： 0.0039904845837 1.00298610472 0.0038894292794 -0.0130948132278
函数库计算均值、标准差、偏度、峰度： 0.0039904845837 1.00298610472 0.0038894292794037777 -0.01309481322786521

相关知识：
偏度：衡量随机变量概率分布的不对称性，是相对于均值不对称程度的度量
偏度可以为0，可以为正、负
偏度为负表示概率密度函数在左侧的尾部比右侧的长，为正则右侧尾部比左侧长

偏度为零表示数值相对均匀分布在平均值的两侧，但不一定意味着一定是对称分布
$% MathType@MTEF@2@1@+= % feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbWexLMBbXgBd9gzLbvyNv2CaeHbl7mZLdGeaGqiVu0Je9sqqr % pepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs % 0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaai % aabeqaamaabaabauaakqaabeqaaiabeo7aNnaaBaaaleaacaaIXaaa % beaakiaab2dacaqGfbGaae4waiaabIcadaWcaaqaaiaadIfacqGHsi % slcqaH8oqBaeaacqaHdpWCaaGaaeykamaaCaaaleqabaGaaG4maaaa % kiaab2facaqG9aWaaSaaaeaacaWGfbGaai4waiaacIcacaGGybGaey % OeI0IaeqiVd0MaaiykamaaCeaaleqabaGaaG4maaaakiaac2faaeaa % caGGOaGaamyraiaacUfacaGGOaGaaiiwaiabgkHiTiabeY7aTjaacM % cadaahbaWcbeqaaiaaikdaaaGccaGGDbGaaiykamaaCaaaleqabaWa % aSGaaeaacaaIZaaabaGaaGOmaaaaaaaaaOGaeyypa0ZaaSaaaeaacq % aH6oWAdaWgaaWcbaGaaG4maaqabaaakeaacqaH6oWAdaqhaaWcbaGa % aGOmaaqaamaaliaabaGaaG4maaqaaiaaikdaaaaaaaaaaOqaaiabeo % 7aNnaaBaaaleaacaaIXaaabeaakiaab2dacaqGfbGaae4waiaabIca % daWcaaqaaiaadIfacqGHsislcqaH8oqBaeaacqaHdpWCaaGaaeykam % aaCaaaleqabaGaaG4maaaakiaab2facaqG9aWaaSaaaeaacaWGfbGa % ai4waiaacIfadaahbaWcbeqaaiaaiodaaaGccaGGDbGaeyOeI0IaaG % 4maiabeY7aTjaadweacaGGBbGaaiiwamaaCeaaleqabaGaaGOmaaaa % kiaac2facqGHRaWkcaaIYaGaeqiVd02aaWbaaSqabeaacaaIYaaaaa % GcbaGaeq4Wdm3aaWbaaSqabeaacaaIZaaaaaaakiabg2da9maalaaa % baGaamyraiaacUfacaGGybWaaWraaSqabeaacaaIZaaaaOGaaiyxai % abgkHiTiaaiodacqaH8oqBcqaHdpWCdaahaaWcbeqaaiaaikdaaaGc % cqGHsislcqaH8oqBdaahaaWcbeqaaiaaiodaaaaakeaacqaHdpWCda % ahaaWcbeqaaiaaiodaaaaaaaaaaa@9C51@ \begin{align} & {{\gamma }_{1}}\text{=E }\!\![\!\!\text{ (}\frac{X-\mu }{\sigma }{{\text{)}}^{3}}\text{ }\!\!]\!\!\text{ =}\frac{E[(X-\mu ){}^{3}]}{{{(E[(X-\mu ){}^{2}])}^{{}^{3}\!\!\diagup\!\!{}_{2}\;}}}=\frac{{{\kappa }_{3}}}{\kappa _{2}^{{}^{3}\!\!\diagup\!\!{}_{2}\;}} \\ & {{\gamma }_{1}}\text{=E }\!\![\!\!\text{ (}\frac{X-\mu }{\sigma }{{\text{)}}^{3}}\text{ }\!\!]\!\!\text{ =}\frac{E[X{}^{3}]-3\mu E[X{}^{2}]+2{{\mu }^{2}}}{{{\sigma }^{3}}}=\frac{E[X{}^{3}]-3\mu {{\sigma }^{2}}-{{\mu }^{3}}}{{{\sigma }^{3}}} \\ \end{align}$

峰度$% MathType@MTEF@2@1@+= % feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbWexLMBbXgBd9gzLbvyNv2CaeHbl7mZLdGeaGqiVu0Je9sqqr % pepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9vqaqpepm0xbba9pwe9Q8fs % 0=yqaqpepae9pg0FirpepeKkFr0xfr=xfr=xb9adbaqaaeGaciGaai % aabeqaamaabaabauaakeaadaWcaaqaaiabeY7aTnaaBaaaleaacaaI % 0aaabeaaaOqaaiabeo8aZnaaCaaaleqabaGaaGinaaaaaaaaaa@4489@ \frac{{{\mu }_{4}}}{{{\sigma }^{4}}}$:
峰度是概率密度函数在均值出峰值高低的特征，通常定义四阶中心矩除以方差的平方减去3。
${{\gamma }_{2}}=\frac{{{\kappa }_{4}}}{\kappa _{2}^{2}}=\frac{{{\mu }_{4}}}{{{\sigma }^{4}}}-3=\frac{\frac{1}{n}\sum\limits_{i=1}^{n}{{{({{X}_{i}}-\bar{X})}^{4}}}}{{{(\frac{1}{n}\sum\limits_{i=1}^{n}{{{({{X}_{i}}-\bar{X})}^{2}}})}^{2}}}-3$