最近太忙了,勉强利用宾馆的两个无聊的晚上,凑成了第二篇,关于正态分布的内容,正态分布是非常重要的前提,分析前需要先检验,然后看是否需要转换。
在进行真正的检验之前,我们很多时候都忽略了进行正态分布的验证。验证是否是正态分布的方法分为两类:
肉眼判断
假设检验
肉眼判断
肉眼判断当然是一句玩笑话,我们不能只凭看两眼数据就认定是否符合正态分布了,我们是通过其与正态分布的数据相比较,并不是直接做个大仙,看上一眼就解决问题了。那么怎么和正态分布的数据做比较呢,比较直观并且肉眼可断的当然是看图了,看数再厉害也要晕。
直方图和密度曲线的比较
第一种方法就是对数据做直方图和密度曲线的图,然后看数据是否基本符合正太分布的特征(也就是那个钟形曲线),这里网上关于 R 验证正态分布的时候使用直方图 hist
,以及密度曲线要用的 density
,其实二者是一回事,本质上看图还是要靠核密度估计,我们先了解一下核密度估计。
核密度估计是非参估计的一种,对我们最重要的是不需要先验知识,即可根据数据的本身特点和性质来进行拟合分布,查看数据分布的特点,最简单的是使用直方图,例如对 iris 的数据的萼片长度的分布进行查看:
layout(matrix(c(1, 2, 3, 4), 2, 2, byrow = TRUE))hist( iris$Sepal.Length, freq = FALSE, breaks = c(seq(4, 8, 0.05)), main = "breaks = 0.05", xlab = NULL)hist( iris$Sepal.Length, freq = FALSE, breaks = c(seq(4, 8, 0.3)), main = "breaks = 0.3", xlab = NULL, ylab = NULL)hist( iris$Sepal.Length, freq = FALSE, breaks = c(seq(4, 8, 0.5)), main = "breaks = 0.5", xlab = "Sepal Length")hist( iris$Sepal.Length, freq = FALSE, breaks = c(seq(4, 8, 0.8)), main = "breaks = 0.8", xlab = "Sepal Length", ylab = NULL)
也就是说,我们划分每组长度的标准(bin 的宽度)不同,肉眼看上去,直方图的形状是有变化的。而且 bin 的宽度不管怎么变化,显示的分布并非是连续的变化,要避免这些问题,就需要引入核密度估计,核密度估计的曲线根据所有点在某个特定位置的权重来计算其在整个分布曲线上的距离,当在某一区域点的数量增加时,估计的概率就会增大。
hist( iris$Sepal.Length, freq = FALSE, xlab = "Sepal Length", main = "Histogram and density of iris sepal length")lines(density(iris$Sepal.Length), col = 2, lwd = 2)
需要注意的是,不管怎样,我们很难得到非常完美的,或者理论上的分布曲线,只要曲线有一个峰值并且大致两侧对称,我们就认为是正太分布了。
python 版本的验证:
import pandas as pdimport statsmodels.api as smimport matplotlib.pyplot as pltimport seaborn as snsiris = pd.read_csv("./data/iris.csv")kwargs = dict(hist_kws={'alpha':.6}, kde_kws={'linewidth':2})plt.figure(figsize=(10,7), dpi= 80)ax = sns.distplot(iris["Sepal.Length"], label="Sepal Length histogram and density", **kwargs)ax.set_xlim([4,8])
注意,形状的变化并非是因为错误,而是因为二者默认的 Kernal Density Estimation(KDE)方法不同,但这对我们来讲,没什么影响,因为这幅图看上去的密度曲线也是近似正态分布的。另外,作图的颜色等的差别并非 R 弱一些,因为我最近有点喜欢 base plot 的语法,所以没有用 ggplot2
, 喜欢 tidyverse 语法的可以参考一下 "geom_density"。
Q-Q plot
相比密度直方图,Q-Q 图可能是更常见更容易理解的一种方式,Q 表示分位点,例如标准的正太分布曲线,0.5 的分位点是零,0.95 的分位点是 1.64。
它的实现就是使用测量数据和理论分布的数据进行作图,如果测量数据是正态分布,那么很明显,它和理论正态分布所成的散点图应该在一条直线上。
在 R 里,base plot 可以直接使用 qqnorm
和 qqline
来实现,但 car
包里 qqPlot
的实际上更贴心
library(car)layout(matrix(c(1, 2), 1, 2))qqnorm(iris$Sepal.Length, pch = 21, col = 'blue', main = 'qqplot in base plot')qqline(iris$Sepal.Length, col = "red", lwd = 3)qqPlot(iris$Sepal.Length, main="qqplot in car package", col="blue", col.lines="red")
car
会贴心的画出 95% 置信区间,以及 outlier(原始数据)。这个与上面提到的一样,不可能完全在一条直线上,大致是在一起就可以了。
python 里的实现也很简单:
import pandas as pdimport statsmodels.api as smiris = pd.read_csv("./data/iris.csv")sm.qqplot(iris["Sepal.Length"], line ='q')
其实单纯统计的细节来讲,Python 远没 R 顺手。毕竟是统计学家编的软件,不能把老本行丢了。
假设检验
上面的肉眼判断毕竟太过主观,要客观的评价,还是需要通过假设检验,用的很多,非常有名的方法为:Shapiro-Wilk 检验,它的零假设是样本来自正态分布的母本,这样的话,我们需要 p 值大于 0.05。
还是上面的数据:
shapiro.test(iris$Sepal.Length)
##
## Shapiro-Wilk normality test
##
## data: iris$Sepal.Length
## W = 0.97609, p-value = 0.01018
非常简单,一行代码搞定,结果让我们非常失望,数据不是正态分布。事实证明,上面的肉眼判断,我的看法是错误的。
python 的 scipy
中也有同样的方法:
from scipy import statsimport pandas as pdiris = pd.read_csv("./data/iris.csv")stats.shapiro(iris["Sepal.Length"])
## (0.9760899543762207, 0.010180278681218624)
尽管值略有差别,充分体现了主观的判断的不靠谱。那我们只能试试怎么将数据转换为正态分布的方法了。
常用的正态转换方法
当然,如果样本量超过了30,或者50更好,正态分布的检验并非十分重要。当然转换为正态分布的方法比较多,下面的一幅图即可概括:
来自 Applied Multivariate Statistics for the Social Sciences
方法这么多,其实也没太多特别的,就是将数据转换一下,例如我们使用最简单的对数转换看一下:
shapiro.test(log10(iris$Sepal.Length))
##
## Shapiro-Wilk normality test
##
## data: log10(iris$Sepal.Length)
## W = 0.98253, p-value = 0.05388
如果是写文章,我想这也就是死里逃生的感觉了,但毕竟比不转换强多了,上面图中列举的方法都比较简单,就不一一尝试了。
参考资料
Amrhein, Valentin, Sander Greenland, and Blake McShane. 2019. “Scientists Rise up Against Statistical Significance.” Nature. https://doi.org/10.1038/d41586-019-00857-9.
DeepAI. 2017. “What Is a Kernel Density Estimation?” https://deepai.org/machine-learning-glossary-and-terms/kernel-density-estimation.
Ed. 2017. “Transforming Data for Normality.” https://www.statisticssolutions.com/transforming-data-for-normality/.
Ford, Clay. 2017. “Understanding Q-Q Plots.” https://data.library.virginia.edu/understanding-q-q-plots/.
“Kernel Density Estimation(KDE).” 2017. https://blog.csdn.net/unixtch/article/details/78556499.
R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Wasserstein, Ronald L., Allen L. Schirm, and Nicole A. Lazar. 2019. “Moving to a World Beyond ‘P < 0.05’.” The American Statistician. https://doi.org/10.1080/00031305.2019.1583913.