研究数据统计分布常用技巧

当我们拿到一组数值型数据,通常除了计算它们的均值、方差以外就是研究它们的统计分布特征。这里,我们将介绍一些研究统计分布的常用技巧。

排序-值曲线

拿到一组数值型数据以后,最简单的一种得到它分布特征的途径就是将这堆数据排序,然后绘制排序-值曲线(Rank value curve)。由于排序-值曲线实际上与数据的累积概率分布曲线构成反函数(后面将给出证明),所以,通过这条曲线我们便可以得知数据的分布特征。

例如,我们用Mathematica中的ParetoDistribution函数生成了一组数,代码如下:

re = RandomReal[ParetoDistribution[1, 1], 1000]

我们就可以通过排序得到它的排序-值曲线,代码如下:

gra = ListLogLogPlot[Sort[re, Greater], PlotRange -> All,  Frame -> True, FrameLabel -> {"Rank", "Value"}]

曲线图如下:

 

我们看到,这些数据中有少数的点具有非常大的数值,而大多数点的数值都很小。当点数非常多的时候,这条曲线就会非常贴近X轴。所以,为了更清楚地看出曲线,我们通常对曲线取双对数坐标,这样就会得到如下所示图:

 

这样的话,整体数据分布的情况就可以看的比较清楚了。而且,由于原数据就是一条幂律分布曲线,所以取双对数坐标后,曲线变成了一条直线。

累积概率密度曲线

另外一种常用的展示数据分布形状的图就是累积概率密度图。这个累积概率密度的定义是:

[math]\displaystyle{ C(x)=Pr\{X\gt x\}=\int_{c}^{\infty} \rho(x) dx }[/math]

因此,根据定义,我们需要扫描所有可能的x数值,然后在数据中数出来比x值大的数字所占的比例,从而得到这条累积概率密度曲线。然而,我们应该选取哪些x数值呢?一种方法是让x等间隔地取值,但这样会得到很多无效的采样点。由于有可能在某个小区间内,根本没有多少数据点。

例如,假设原始数据为{0.1,0.5,1.4,1.5},如果我们选取小区间长度为0.1,那么采样点x=0.1,0.2,0.3,0.4,0.5,...,我们会看到0.2,0.3,0.4就是无效的采样点,因为比0.2大的数据点数跟比0.3大的数据点数一样多。

一种最简单的选取数值x的方法就是用原始数据点本身。在上个例子中,我们只需要让x=0.1,0.5,1.4,1.5就可以了。这样的计算是最有效的。我们还是用帕累托分布的1000个随机数来举例说明。Mathematica的源代码如下:

sorted = Union[re];
cnt = Count[re, x_ /; x >= #] & /@ sorted;
ListPlot[Transpose[{sorted, cnt}], PlotRange -> All,  Frame -> True, FrameLabel -> {"Value", "Count"}]

这样就可以得到累积概率密度曲线(注意总坐标并没有归一化)

 

以及双对数坐标下的累积概率密度曲线

 

累积密度曲线和排序-值曲线之间的关系

仔细观察同一个分布的累积概率密度曲线与排序-值曲线就会发现,它们的形状实际上具有一定的对称性。如下图所示:

 

 

不难看出,如果将一条曲线按照y=x斜对角线对折,就可以得到另一条曲线了。如果你注意观察就会发现,累积概率密度曲线的纵坐标取值范围为[0,1],而排序-值曲线的取值范围为[0,N],这里N就是数据点的个数。但如果我们将累积概率密度曲线纵坐标放大N倍,那么它就与排序-值沿着45度对角线严格对称了。

事实上,我们可以严格地论证,只要数据点N的个数足够大,两条曲线的确互为反函数的关系。下面,我们将给予证明。

假设一组数据有N个点,它们的排序-值曲线写为V(r),其中r为数据点的排序,V(r)为该数据点的数值。当N足够大的时候,V(r)就是一条连续的曲线了。那么,我们来计算这组数据的累积概率密度函数。 根据定义,累积概率函数C(x)定义为Pr{X>x}。但是,由于我们仅有有限样本,所以我们不得不用数据的频率来代替概率。也就是说我们只能这样计算C(x):

[math]\displaystyle{ C(x)=N(\gt x)/N }[/math]

其中,N(>x)表示大于x的数据点个数。也就是说C(x)为大于x的数据点比例。这个比例可以通过V(r)曲线求出来。注意到V(r)曲线是严格递减的,所以比x数值大的数据点一定在[math]\displaystyle{ V^{-1}(x) }[/math]的左侧。这里[math]\displaystyle{ V^{-1}(x) }[/math]表示的是V(r)的逆函数在x这个数值处的取值。而[math]\displaystyle{ V^{-1}(x) }[/math]左侧的数据点个数显然就是[math]\displaystyle{ V^{-1}(x) }[/math](因为C(r)曲线在这一点的横坐标就是它)。所以,我们可以得到:

[math]\displaystyle{ NC(x)=V^{-1}(x) }[/math]

也就是说,如果将曲线C(x)的纵轴扩大N倍,那么它就是V(r)函数的反函数。

概率密度函数

虽然排序-值曲线以及累积概率密度曲线包含了数据分布的全部信息。但是在实际应用中,我们通常需要知道数据的概率密度函数曲线。但是在绘制概率密度曲线的时候,我们需要将数据可能取值的区间划分为若干等长度的小区间并计算每个小区间内数据点的个数。然而,小区间长度的选择是一个非常重要而敏感的问题。因为,如果小区间长度过大,则虽然每个区间落入的数据点比较多,但是曲线的细节信息就会被丢失;反过来,如果小区间长度过小,那么每个小区间中落入的数据点就很少,所以我们几乎看不出曲线的形状。

为了更直观地理解这个问题,让我们还是用1000个满足Pareto分布的随机数来举例说明。当我们取小区间长度dx=0.01的时候,我们得到的概率密度曲线如下:

 

而当我们选取dx=0.1的时候,我们得到的概率密度曲线就会变成如下的样子:

 

得到上面两张图的源代码如下:

dx = 1;
gra = ListLogLogPlot[  Table[{x, Count[re, y_ /; y > x && y <= x + dx]/Length[re]}, {x, 1, 1000, dx}], PlotRange -> All, Filling -> Axis, Frame ->   True, FrameLabel -> {"x", "Probability"}]

比较这两张图,我们就会发现它们的差异还是相当明显的。所以,在实际应用中,如果数据量很小,那么我们通常会试验各种可能的dx数值,以期得到最后的最好结果。

Log Bin方法

当我们研究满足幂律分布或者具有长尾特征的数据的时候,通常会采用取LogBin的方法从而完成对数据的统计。让我们还是以1000个满足Pareto分布的随机数来举例说明。当我们把dx设置为1的时候,绘制出概率密度曲线,并对其取双对数就会得到如下图形:

文件:Paretodistributiondensitycurvedxloglog0.1.png

按理说,由于数据来源于Pareto分布,因此它的概率密度函数应该是一个幂律函数,也就是说概率密度曲线在取过双对数之后应该形成一条直线。但我们看到,图中所示的曲线在偏离直线的程度还是很大的。这其中的一个很大原因在于由于数据分布满足长尾特征,所以尾部的数据非常稀疏,当我们的小区间长度取得比较小的时候,就有可能导致尾端有很多小区间中得不到足够多的数据累积,从而使得尾部(右端)的噪声较大,从而使得曲线的形状失真。

为了克服这个问题,我们通常采用LogBin方法。也就是说并不是对原始数据的可能取值做等区间的划分。而是对原始数据取对数以后做等间隔的划分,从二再次统计每个Log小区间中的数据点个数,这个时候,我们通常可以得到比较好的结果。如下图,当我们取对数后,取小区间长度为dx=0.5,我们得到的概率密度曲线如下图

文件:Paretodistributionlogbin0.5.png

我们看到,得到的曲线会更加光滑,而且曲线更接近直线。

得到上面曲线的代码如下:

dx = 0.5;
logre = Log[re];
gra = ListLogLogPlot[ Table[{Exp[x], Count[logre, y_ /; y > x && y <= x + dx]/Length[re]}, {x, Min[logre], Max[logre], dx}], PlotRange -> All, Filling -> Axis, Frame -> True, FrameLabel -> {"x", "Probability"}]