用R语言进行数据分析：检验一个数据集的分布

我们可以用很多方法分析一个单变量的数据集。 最简单的办法就是直接看数字。利用函数 summary 和 fivenum 会得到 两个稍稍有点差异的汇总信息。 此外，stem (“茎叶”图)也会反映整个数据集的数字信息的。

     > attach(faithful)
     > summary(eruptions)
        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
       1.600   2.163   4.000   3.488   4.454   5.100
     > fivenum(eruptions)
     [1] 1.6000 2.1585 4.0000 4.4585 5.1000
     > stem(eruptions)
     
       The decimal point is 1 digit(s) to the left of the |
     
       16 | 070355555588
       18 | 000022233333335577777777888822335777888
       20 | 00002223378800035778
       22 | 0002335578023578
       24 | 00228
       26 | 23
       28 | 080
       30 | 7
       32 | 2337
       34 | 250077
       36 | 0000823577
       38 | 2333335582225577
       40 | 0000003357788888002233555577778
       42 | 03335555778800233333555577778
       44 | 02222335557780000000023333357778888
       46 | 0000233357700000023578
       48 | 00000022335800333
       50 | 0370

茎叶图和柱状图相似，R 用函数 hist 绘制柱状图的。

     > hist(eruptions)
     ## 让箱距缩小，绘制密度图
     > hist(eruptions, seq(1.6, 5.2, 0.2), prob=TRUE)
     > lines(density(eruptions, bw=0.1))
     > rug(eruptions) # show the actual data points

许多不错的密度图都是用 density 绘制的。在这个例子中， 我们加了一条由 density 产生的曲线。你可以用试错法（trial-and-error） 选择带宽 bw（bandwidth）因为默认的带宽值让密度曲线过于平滑 （这样做常常会让你得到非常有“意思”的密度分布）。(现在已经有 自动的带宽挑选办法，在这个例子中 bw = "SJ"给出了不错的结果。)

我们可以用函数 ecdf 得到一个数据集的 经验累积分布（empirical cumulative distribution）

     > plot(ecdf(eruptions), do.points=FALSE, verticals=TRUE)

这个分布和其他标准分布差异很大。 这种分布说明火山爆发都在3分钟以后发生。 让我们拟合一个正态分布，并且重叠上面得到的经验累积密度分布。

     > long <- eruptions[eruptions > 3]
     > plot(ecdf(long), do.points=FALSE, verticals=TRUE)
     > x <- seq(3, 5.4, 0.01)
     > lines(x, pnorm(x, mean=mean(long), sd=sqrt(var(long))), lty=3)

分位比较图（Quantile-quantile (Q-Q) plot）可以让我们更细致地研究二者的吻合程度。

     par(pty="s")       # 设置一个方形的图形区域
     qqnorm(long); qqline(long)

上述命令得到的QQ图表明二者还是比较吻合的，但右侧尾部 偏离期望的正态分布。我们可以用 t 分布 得到一些模拟数据重复上面的过程

     x <- rt(250, df = 5)
     qqnorm(x); qqline(x)

得到的QQ图常常会出现偏离正态期望的长尾部(如果是随机样本)。 我们可以用下面的命令针对特定的分布绘制 Q-Q图

     qqplot(qt(ppoints(250), df = 5), x, xlab = "Q-Q plot for t dsn")
     qqline(x)

最后，我们需要更为正规的正态性检验。 R提供了 Shapiro-Wilk 检验

     > shapiro.test(long)
     
              Shapiro-Wilk normality test
     
     data:  long
     W = 0.9793, p-value = 0.01052

和 Kolmogorov-Smirnov 检验

     > ks.test(long, "pnorm", mean = mean(long), sd = sqrt(var(long)))
     
              One-sample Kolmogorov-Smirnov test
     
     data:  long
     D = 0.0661, p-value = 0.4284
     alternative hypothesis: two.sided

(注意一般的统计分布理论（distribution theory）在这里可能无效， 因为我们用同样的一个样本对正态分布的参数进行 估计的。)