2018-09-03

数据描述性分析-数据分布

之前我们介绍了描述性分析中的一些描述统计量，它们刻画了数据的主要特征，而要对数据的总体情况作全面的描述，就要研究数据的分布。对数据分布的主要描述方法有总体分布、直方图、茎叶图等。

1 分布函数

说分布函数之前，我们要先介绍随机事件（简称事件）和随机变量。

我们做一项随机试验，随机试验的每一个结果可以看作是一个样本点，随机试验的所有可能结果的全体构成一个样本空间，我们把事件定义为样本空间的一个子集，它包含若干样本点，该事件发生当且仅当它所包含的样本点中有一个发生。但并不是说样本空间的一切子集都是事件，因为这将对给定概率带来困难，另一方面，我们要把问题中感兴趣的事件都包括进来。为此，我们的做法是，样本空间的一些子集构成一个集合，如果这个集合包含空集和全集，并且集合满足对逆、并、交、差的可列次运算封闭，那么就把这个集合叫做事件域，事件域中的元素就叫做事件。

那随机变量又是什么？它跟随机事件有什么联系呢？在随机现象中，有很大一部分问题与数值发生关系，例如在产品检验问题中，我们关心的是抽样中出现的废品数；在机票超售问题中我们关心的是某航班实到旅客数；在电话问题中我们关心的是某段时间中的话务量，它与呼叫次数及各次呼叫占用交换设备的时间长短有关。此外如测量时的误差，气体分子运动的速度，信号接收机所收到的信号，也都与数值有关。有些初看起来与数值无关的随机现象，也常常能联系数值来描述。例如在掷硬币问题中，每次出现的结果为正面或反面，与数值没有关系。但是我们能用下面方法使它与数值联系起来，当出现正面时对应数“1”，而出现反面时对应数“0”，为了计算$n$次投掷中出现的正面数就只须计算其中“1”出现的次数了。

一般地，如果$A$为某个随机事件，则一定可以通过如下示性函数使它与数值发生联系

\[1_A=\begin{cases} 1,如果A发生\\ 0,如果A不发生 \end{cases}\]

总之，这些例子中，试验的结果能用一个数来表示，这个数是随着试验结果的不同而变化的，即它是样本点的一个函数，这种量就称为随机变量。随机变量是样本空间到实数域的一个函数。对于试验结果取离散值的随机变量，我们叫它离散型随机变量，对于取连续值的随机变量，我们叫它连续型随机变量。需要强调的是，离散型随机变量与连续型随机变量没有交集，但它们的并也不是全集。

现在可以给出分布函数的定义了，设$X$是一个随机变量，对任意的实数$x$,令\[F(x)=P\left\{X\leq x\right\},x\in(-\infty,+\infty)\]则称$F(x)$为随机变量$X$的分布函数。

直观来讲，分布函数就是变量小于等于某个特定值的概率，分布函数是随机变量最重要的概率特征，它可以完整地描述随机变量的统计规律，并且决定随机变量的其他一切概率特征。

在R中，提供了计算一些分布的分布函数，在这里我们分别展示一个典型的离散型随机变量和连续型随机变量的分布函数的例子

正态分布（连续型随机变量的分布函数）

设$\mu$是均值，$\sigma^2$是方差，对于服从正态分布的随机变量$X$来说，其分布函数为\[F(x)=\frac{1}{\sqrt{(2\pi)\sigma}}\begin{matrix}\int_{-\infty}^{x}e^{-\frac{(t-\mu)^2}{2\sigma^2}}dt\end{matrix}\] 在R中的调用格式为

pnorm(x,mean=o,sd=1,lower.tail=TRUE,log.p=FALSE)

其中

x是由数值型变量构成的向量；

mean是要计算的正态分布的均值，sd是标准差，默认是0和1，可根据实际情况更改；

lower.tail是逻辑变量，当它为真时，分布函数的计算公式为$F(x)=P(X\leq x)$，当它为假时，分布函数的计算公式为$F(x)=P(X\ge x)$;

log.p是逻辑变量，当它为真时，函数的返回值不再是正态分布，而是对数正态分布。

泊松分布（离散型随机变量的分布函数）

分布函数为\[P(X=k)=\frac{\lambda e^-\lambda}{k!}\]

在R中的调用格式为

dpois(x,lambda,log=FALSE)

其中，lambda是Poisson分布的参数$\lambda$.其余参数的意义与正态分布参数意义相同。其他的分布函数也有类似的结果，不再累述。

2 直方图、经验分布图与QQ图

2.1 直方图

直方图是用来描述数据分布的常用工具之一，它是由若干矩形构成的图，在构造直方图时，先把数据值域分成若干个等长区间，之后利用这些区间把数据分组，每组作一个矩形，其高和数据落入该组的频数成比例，其底为所属区间。

在R软件中，用函数hist()画样本直方图，调用格式为hist()或者

hist(x, breaks = "Sturges",
     freq = NULL, probability = !freq,
     include.lowest = TRUE, right = TRUE,
     density = NULL, angle = 45, col = NULL, border = NULL,
     main = paste("Histogram of" , xname),
     xlim = range(breaks), ylim = NULL,
     xlab = xname, ylab,
     axes = TRUE, plot = TRUE, labels = FALSE,
     nclass = NULL, warn.unused = TRUE, ...)

在这里介绍几个常用参数：

x是由样本构成的向量；
breaks规定直方图的组距，由以下几种形式给出：

向量，给出直方图的起点、终点、组距；

数，定义直方图的组距；

字符串（默认状态）；

函数，计算组距的宽度；

freq是逻辑变量：

TRUE绘出频率直方图；

FALSE绘出密度直方图；

col表示直方图中填充的颜色；
plot是逻辑变量：

TRUE表示给出直方图（默认状态）；

FALSE表示列出直方图的各种结果但不绘图。

举一个R代码的例子如下

x<-c(1,25,43,36,8,52,10,12,18,16,15,35,12,13,15,8,19,20)
hist(x) # x是由样本构成的向量

2.2 核密度估计函数

当给定一个样本集，怎么得到该样本集的分布密度函数？解决这个问题有两种方法，一是用参数估计方法，简单来讲，就是假定样本集符合某一概率分布，然后根据样本集拟合该分布中的参数，但由于参数估计方法中需要加入主观的先验知识，往往很难拟合出真实分布的模型；二是用非参数的方法，非参数估计并不加入任何先验知识，而只是根据数据本身的特点和性质来拟合分布，这样能比参数估计方法得出更好的模型。核密度估计就是非参数估计中的一种。

我们刚才用直方图的方法来直观地展现样本的分布情况，但是用直方图表述存在一些问题：

因为我们使用有限的分组宽度代替单个数据点的位置，所以直方图会丢失信息；
直方图不是唯一的；
直观上看，直方图不够光滑；
直方图不能很好地处理一些异常点，某个远离大部分点的异常点会造成图形中间的空白。

核密度估计作为直方图配套使用的一种操作，则能很好地处理上述问题。在R中可调用的是density()函数。简单的调用格式为density(x，bw="nrd0"),x是由样本构成的向量，bw是带宽，当bw为默认值时，R软件会画出光滑的曲线。

w<-c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0)
hist(w,freq=FALSE) # x是由样本构成的向量
lines(density(w),col="blue") # 用蓝线标出核密度估计的曲线

2.3 经验分布函数

对于一个总体，如果它的分布函数未知，则可以用经验分布函数估计。我们给出经验分布函数的定义：

设$X_1,...X_n$为取自总体分布函数为$F(x)$的样本，$X_{(1)}\leq...\leq X_{(n)}$为其顺序统计量，则称\[F_n(x)=\frac{1}{n}\sum_{i=1}^n I_{\left\{X_i<x\right\}}=\begin{cases}0,x\leq X_{(1)}\\\frac{k}{n},X_{(k)}<x\leq X_{(k+1)},k=1,...,n-1\\1,x>X_{(n)}\end{cases}\]为样本$X_1,...X_n$的经验分布函数。经验分布函数具有良好的性质，比如当$n$充分大时，$F_n(x)$与$F(X)$的值在概率意义下是非常接近的。

在R中，用函数plot(ecdf())绘制样本经验分布图,调用格式为

plot(ecdf(x),...,verticals=FALSE,do.p=FALSE)

* verticals是逻辑变量，值为真时画竖线，否则不画竖线；
* do.p是逻辑变量，值为真时不画点处的记号，否则画记号。

例如

x<-c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0)
plot(ecdf(x)) # x是由样本构成的向量

2.4 QQ图

不管是直方图还是经验分布图，要从比较上鉴别样本是否近似于某种类型的分布较为困难，而QQ图可以帮助我们分析样本的分布是否近似于某种类型的分布。

为什么QQ图能做到如此地步呢，我们知道，分位数函数是分布函数的反函数，又知道一个函数与它的反函数具有相同的单调性，那么我们要看一个样本是不是近似服从正态分布，只需要看样本的顺序统计量和对应的标准正态分布的分位数序列是不是呈近似线性关系即可。（详细的理论证明在下文正态性检验中给出）

在R语言中，函数qqnorm()和qqline()提供了画正态QQ图和相应直线的方法,例如

x<-c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0)#15位同学的体重数据
qqnorm(x) # x是由样本构成的向量
qqline(x) # x是由样本构成的向量

从图中可以看出，样本点与直线基本契合，因此该样本数据可以看做是来自正态总体。

3 茎叶图、箱线图

3.1 茎叶图

茎叶图类似于直方图，只不过用数据代替矩形。在R中，用stem()函数作茎叶图，调用方法为

stem(x,scale=1,width=80,atom=1e-08)
     x是数据向量
     scale控制绘出茎叶图的长度
     width是绘图的宽度
     atom是容差

例如

x<-c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0)#15位同学的体重数据
stem(x) # x是由样本构成的向量

## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##   4 | 7
##   5 | 0779
##   6 | 22444779
##   7 | 25

在茎叶图中，数据的十位数部分表示“茎”，个位数部分作为“叶”。

3.2 箱线图

直方图和茎叶图包含了大量的样本信息，但它们并没有对样本作过多的加工和整理，而我们有时需要用少数几个统计量来对大量的原始数据进行简单的概括，而其中有五个数比较重要，分别是：最值、中位数、上（下）四分位数。把这五个数画在一个图上，在下四分位数点和上四分位数点之间画一个矩形箱子，箱子上边线是上四分位数，下边线是下四分位数，在中位数处画一横线，这幅图就叫做箱线图，箱子上下伸出的部分表示数据的散布范围，超出此范围的为异常点。

在R语言中，用boxplot()函数作箱线图。调用格式有如下两种：

boxplot(x,...outline=TRUE,horizontal=FALSE,add=FALSE)

* outline是逻辑变量，值为真时，标明异常点；
* horizontal是逻辑变量，值为真时，把箱线图绘成水平状；
* add是逻辑变量，值为真时，在原图上画图，否则替换上一张图。

或者

boxplot(formula,data=NULL,...,subset,na.action=NULL)

* formula是公式，如$y~grp$,这里y是由数据构成的数值型向量
* grp是数据的分组

举一个例子

x<-c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0)#15位同学的体重数据
boxplot(x) # x是由样本构成的向量

4 正态性检验和分布拟合检验

4.1 正态性检验

4.1.1 W检验

关于正态检验的方法有很多，这里介绍两种常用且有效的检验方法，它们是Wilk-Shapiro的W检验和D’Agostino的D检验：当样本量$n\leq50$时，用W检验；当样本量$n>50$时用D检验。以下所需要的假设检验和参数估计的基础知识在其他章节给出。

为了方便，本节始终假设$X_1,...X_n$是来自总体为$F(X)$的$IID$（相互独立）的样本，现感兴趣的假设为\[H_0:F(x)\in \left\{\phi(\frac{x-\mu}{\sigma}):\mu\in R,\sigma>0\right\}\]

W检验的检验统计量为\[W=\frac{[\sum_{k=1}^{[n/2]} a_k(n)(X_{(n+1-k)}-X_{(k)})]^2}{\sum_{i=1}^n (X_{(i)}-\bar X)^2}\]

其中$X_{(1)}\leq ...\leq X_{(n)}$为顺序统计量，$\left\{a_k(n)\right\}$是一组给定的与$n$有关的常数，可通过查W检验统计量系数表获得。W检验的拒绝域为\[\left\{W<c\right\}\]其中临界值$c$满足$P\left\{W<c|H_0\right\}=\alpha$

此检验的理由大致如下：设$X_1,...X_n$为来自$N(\mu,\sigma^2)$的IID样本，如记$Y_i=\frac{X_{(i)}-\mu}{\sigma}$,$m_i=EY_i$,$\epsilon=X_{(i)}-EX_{(i)}$，则知$m_i$是与$\mu,\sigma$无关的常数，且有\[X_{(i)}=\mu+\sigma m_i+\epsilon\]

由于当$n$较大时，$X_{(i)}\approx EX_{(i)}$，因而$\epsilon$较小。于是$n$个点$\left\{(X_{(i)},m_i):i=1,2,...n \right\}$近似在一条直线上。为判断这$n$个点是否在一条直线上，我们可采用下面的$R^2$统计量\[R^2=[\frac{\sum_{i=1}^n(X_{(i)}-\bar X)(m_i-\bar m)}{\sqrt {\sum_{i=1}^n(X_{(i)}-\bar X)^2\sum_{i=1}^n(m_i-\bar m)^2}}]^2\]

其中，$\bar m=\sum_{i=1}^nm_i/n$$

可以看出，$R^2$统计量是相关系数的平方，相关系数可以衡量两组数据是否线性相关，因此可以用$R^2$来判断$n$个点是否在一条直线上。

当$H_0$成立，即$X_i$服从$N(\mu,\sigma^2)$时，$R^2$应接近1，于是，当$R^2$远离1时，我们有理由拒绝$H_0$。由于$N(0,1)$是对称分布，所以当$H_0$成立时，$(Y_1,...Y_n)$与$(-Y_1,...,-Y_n)$有相同的联合分布，从而$Y_k$与$-Y_{n+1-k}$也有相同的分布，故$m_k=-{m_n+1-k}$,$\bar m=0$。于是，\[R^2=\frac{[\sum_{i=1}^nX_{(i)}m_i]^2}{\sum_{i=1}^n(X_{(i)}-\bar X)^2\sum_{i=1}^n m_i^2}=\frac{[\sum_{k=1}^{[n/2]}m_{n+1-k}(X_{(n+1-k)}-X_{(k)})]^2}{\sum_{i=1}^n(X_{(i)}-\bar X)^2\sum_{i=1}^nm_i^2}\] 如取\[a_k(n)=\frac{m_{n+1-k}}{\sqrt{\sum_{i=1}^nm_i^2}}\]

则知$W=R^2$，所以我们可以用$W$检验统计量来检验正态性。在R语言中，函数shapiro.test()提供W统计量和相应的p值，当$p$值小于某个显著性水平$\alpha$时，则认为样本不是来自正态分布总体；否则认为样本是来自正态分布总体。函数调用格式为shapiro.test(x)。例如

x<-c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0,53.8,45.2,56,70,57,61,63,65,54)
shapiro.test(x) # x是由样本构成的向量

## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.98305, p-value = 0.9447

假设显著性水平为0.05，那么根据结果，我们可以看出，W统计量为0.96862，十分接近于1，并且p值=0.8371>0.05,接受原假设,认为总体来自正态总体。

4.1.2 D检验

W检验虽然有效，但可惜的是它只适用于容量为3-50的样本，主要原因是，对超过50的样本容量$n$，我们很难计算其分位数值。为此，$D'Agostino$提出了D检验，D检验适用于数据量比较大的情况。D检验的检验统计量为\[D_n=\frac{\sum_{i=1}^n(i-\frac{n+1}{2})X_{(i)}}{n^{3/2}\sqrt{\sum_{i=1}^n(X_i-\bar X)^2}}\]

当它远离1时，我们有理由拒绝正态性假设，可以证明，当$H_0$成立时，$E(D_n)=0.28209479,\sqrt{Var(D_n)}=0.02998598$，故采用如下标准化随机变量\[Y_n=\frac{\sqrt n(D_n-0.028209479)}{0.02998598}\]作为检验统计量

且其拒绝域为\[\left\{Y_n<c_1\right\}\cup\left\{Y_n>c_2\right\}\]现在关于D检验的分位数也已有表可查。

在R语言中，用dago.Test()函数进行D检验，该函数在fBasics程序包中

#install.packages("fBasics")#
library(fBasics)
set.seed(1)#设置随机种子，以使得后边每次生成的随机数序列相同
x=rnorm(1000,0,1)#生成1000个服从标准正态分布的随机数
dagoTest(x) # x是由样本构成的向量

## 
## Title:
##  D'Agostino Normality Test
## 
## Test Results:
##   STATISTIC:
##     Chi2 | Omnibus: 0.0726
##     Z3  | Skewness: -0.2492
##     Z4  | Kurtosis: 0.1025
##   P VALUE:
##     Omnibus  Test: 0.9643 
##     Skewness Test: 0.8032 
##     Kurtosis Test: 0.9184 
## 
## Description:
##  Thu Jul 05 17:43:42 2018 by user: Lenovo

结果输出模型全局似然比检验结果（越小越好）、偏度系数、峰度系数以及他们各自的p值。首先p值都远大于0.05，其次，偏度系数和峰度系数都接近于0，因此接受原假设，认为总体来自正态总体。多次试验发现数据量越大，检验结果越精确。

4.2 经验分布的Kolmogorov-Smirnov检验

经验分布函数$F_n(x)$是总体分布函数$F(x)$的估计。经验分布拟合检验的方法是检验经验分布$F_n(x)$与假设的总体分布函数$F_0(x)$之间的差异。这种差异用距离来度量，Kolmogorov-Smirnov统计量是计算$F_n(x)$与$F_0(x)$的距离D:\[D=\sup_{-\infty\to\infty}|F_n(x)-F_0(x)|\]

在R语言中，函数ks.test()给出了Kolmogorov-Smirnov检验，其基本格式为ks.test(x,y)。其中x为样本构成的向量，y为原假设的数据向量或是描述原假设的字符串。例如

x<-rt(100,5) # 生成100个自由度为5的t分布随机数
ks.test(x,"pf",2,5) # 检验x是否服从自由度为（2，5）的F分布

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.50248, p-value < 2.2e-16
## alternative hypothesis: two-sided

根据p值，我们可以看出，2.2e-16<0.05,所以有理由拒绝原假设，即x不服从F(2,5)的分布。

5 参考文献

[1]薛毅. 统计建模与R软件[M]. 清华大学出版社, 2007. [2]王兆军. 数理统计教程[M]. 高等教育出版社, 2014. [3]李贤平.基础概率论[M].高等教育出版社，2010

源文件下载在线调试

返回：常用统计学, 数据描述性分析