在数据分析项目中,变量的统计分布是各种算法选择的基础;在时序数据分析中,不同时间窗口(周期)内统计分布变化通常也是重要的特征量。例如,在风电行业中,风速分布估算是发电量估算的基础。风速的分布存在典型右偏,不符合正态分布,但到底用Weibull分布、Rayliegh分布还是截断的正态分布更好,目前还没有定论(虽然通常采用Weibull分布)。数据的统计分布拟合包括3类问题。
1)分布的非参数化拟合;
2)单个分布的拟合(参数化);
3)混合分布的拟合。
R语言的stats基础包里面有个函数——density函数,见下面代码和图3-8。
图3-8 Old Faithful泉喷发时间的概率密度
1.bw(Band width)参数
bw参数共有5种选择,见表3-2。
表3-2 bw的5种选择
在R console,直接输入函数式就可以快速看到非内置函数的源代码,以bw.nrd0函数为例:
上面的第6、7行语句使用了R语言技巧,如果条件语句不成立(中间存在NA),则执行下面的语句(lo<-hi)||(lo<-abs(x[1L]))||(lo<-1),这3个“||”表示从左到右边执行,直到第一个非空非零的表达式成立。例如,如果我们将hi设置为NA,则lo<-abs(x[1L])被首选执行,如果x[1L]不为空,lo就赋值成功。
2.核函数的选择
核函数有
不同核函数的对比见下述代码与图3-9。
3.拟合后的概率密度
可以用函数approx做插值,求新的点的概率密度。
图3-9 不同核函数的对比
R里面有fitdistrplus包,针对非剪切(non-censored)数据,其核心函数为
生存分析中常常存在数据剪切,核心函数为
1.概率分布类型的选择
可以对实际数据通过自举法估算skewness、kurtosis等指标,看它们与哪个已知的分布比较接近。图3-10所示的fitdistrplus包的groundbeef数据集,是关于一份儿童牛肉馅饼重量(单位是g)的市场调查数据,有254个样例。descdist函数给出的样本宽度是0.74(右偏),峰度略大于3,从Cullen-Frey图可看出,用Weibull、Gamma或Lognormal等右偏分布可以较好地拟合该样本分布。
2.拟合准则
度量准则有4种:最大似然估计(Maximum Likelihood Estimation,MLE)、矩匹配估计(Moment Matching Estimation,MME)、分位数匹配估计(Quantile Matching Estimation,QME)和最大拟合优度估计(Maximum Goodness-of-fit Estimation,MGE),默认采用MLE准则。
下面分别用Weibull、Gamma或Lognormal分布采用MLE准则去拟合,并绘制拟合后的概率密度图、Q-Q图、累积分布图和P-P图。
图3-10 实际数据统计量与常见分布的比对
拟合结果如图3-11所示,从概率密度图和累积分布图看,三个分布都可以较好地拟合样本数据,但从P-P图看,三个分布都没有很好地拟合出分布的中心区域,从Q-Q图看,Weibull和Gamma对右侧尾部拟合较好。在食物过量风险的业务上下文中,Weibull和Gamma可能是合适的选择。
图3-11 Weibull、Gamma和Lognormal的拟合结果
有8种拟合优度评估指标,见表3-3,其中 F ( x )是拟合模型的累积概率密度函数(Cumulative Distribution Function,CDF), F n ( x )是实际数据的累积概率密度函数(对变量数值排序后计算得到)。
表3-3 拟合优度评估指标
不同拟合优度统计量对应着不同的关注点,例如,AD平等对待尾部和中心部的拟合度。
3.fitdistr的代码
MASS包中的fitdistr函数比较典型,采用最大似然拟合方式,用stats::optim函数进行参数寻优。如图3-12所示,通过阅读,读者可以更清晰地了解背后的计算逻辑,同时加深R语言中quote、formals、match.call、eval、names(list(…))等函数用法的直观了解。
R语言有mixtools、fitdistrplus、mclust、flexmix等包,这里主要介绍mixtools,对于参数化分布,采用EM算法做极大似然估计。对于一个 r 维的样本 x i ,服从 m 个组分的混合分布,假设 x i 各个维度独立,第 k 维的概率密度函数为 f jk (·),则样本 x i 的似然度为
式中,下标 i 、 j 、 k 分别代表样本(individual)、组分(component)、变量坐标顺序(coordinate)。
以Old Faithful泉的喷发间隔数据(见图3-8)为例,演示混合概率分布的拟合。图3-13所示为拟合结果,左图中实线是Gaussian参数化拟合结果,虚线是半参数化拟合结果(bw=4);右图中实线是半参数化拟合结果(bw=1),虚线是半参数化拟合结果(bw=6)。
图3-12 fitdistr的源代码片段
图3-13 数据的拟合结果
除了本节介绍的统计分布的选择和参数拟合方法,可以进一步思考:如何度量两个概率分布的相似度或距离呢?KL散度(Kullback-Leibler divergence)、f-divergence、Wasserstein distance这些指标,读者可以自己查阅理解。