QQ-plot深入理解与实现


26JUN


June 26, 2013



最近在看关于CSI(Channel State Information)相关的论文，发现论文中用到了QQ-plot。Sigh！我承认我是第一次见到这个名词，异常陌生。维基百科给出了如下定义：

“在统计学中，QQ-plot(Q代表分位数Quantile)是一种通过画出分位数来比较两个概率分布的图形方法。首先选定区间长度，点(x,y)对应于第一个分布(x轴)的分位数和第二个分布(y轴)相同的分位数。因此画出的是一条含参数的曲线，参数为区间个数。如果被比较的两个分布比较相似，则其QQ图近似地位于y=x上。如果两个分布线性相关，则QQ图上的点近似地落在一条直线上，但并不一定是y=x这条线。QQ图同样可以用来估计一个分布的位置参数。”

这段话刚开始看的时候，的确不是很清楚，难以理解。我也在网上找了一些资料，最有用的当属网上的一本在线电子书​​《Online Statistics Education: An Interactive Multimedia Course of Study》​ ​，里面的Chanpter8专门有讲解QQ-plot。本文中主要借鉴了这门书中的内容，以更浅显易懂的语言来讲清楚QQ-plot，我学习的过程中也用Matlab做了一些试验，文中将代码一并附上。

QQ-plot其实是Quantile-Quantile Plot的缩写，Quantile分位现在理解没有关系，看到最后你就会理解它的意思了。QQ-plot的目的是什么呢？是为了验证两组数据的分布是否相同或者相似，因此在实际中很多情况都会用到。为了讲清楚QQ-plot，我们先来介绍另外两种以图形的方式评价数据分布情况的方法：直方图(histogram)和 经验累积分布函数(empirical cumulative distribution function, eCDF)。

我们考虑一个随机变量X服从[0,1]区间内均匀分布，我们任取n个数据{ x1,x2...,xnx1,x2...,xn }。本例中n=100，直方图频率分布如图1所示。直方图的概率分布与bins的个数有关(bins为10,5,3)。不同的bins对应的图形也不同，图bins=10的时候还呈现锯齿状，但是bins=3的时候就趋于平稳，所以根据直方图来看累积分布不是很靠谱。随后，我们又使用eCDF对数据进行分析，如图2所示。黄色部分即为eCDF与理论CDF的误差，根据大数定理，当n取值越大，误差越小。

​​​​

图1. 直方图统计

​​​​

图2. eCDF vs 理论CDF

相关代码：




1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16


17


18


19


20


21


22


23


24


25


26


27


28




​​data = unifrnd(0,1,1,100)';​​​​%生成100个再[0,1]均匀分布的随机数​​


​​figure​​


​​subplot(3,1,1);​​


​​hist(data,10);​​​​%bins=10​​


​​xlabel(​​​​'x'​​​​);​​


​​ylabel(​​​​'Frequency'​​​​);​​


​​subplot(3,1,2);​​


​​hist(data,5);​​​​%bins=5​​


​​xlabel(​​​​'x'​​​​);​​


​​ylabel(​​​​'Frequency'​​​​);​​


​​subplot(3,1,3);​​


​​hist(data,3);​​​​%bins=3​​


​​xlabel(​​​​'x'​​​​);​​


​​ylabel(​​​​'Frequency'​​​​);​​


​​theory_y=data;​​


​​figure​​


​​[F,X]=ecdf(data);​​​​%ECDF​​


​​plot(X,F,​​​​'-k'​​​​,​​​​'LineWidth'​​​​,3);​​


​​hold on;​​


​​plot(data,theory_y,​​​​'-b'​​​​,​​​​'LineWidth'​​​​,3);​​


​​legend(​​​​'Empirical CDF'​​​​,​​​​'Theoretical CDF'​​​​,2);​​


​​hold on;​​


​​%------两曲线之间填充颜色-------​​


​​theory_y=sort(theory_y);​​


​​theory_y=[theory_y(1);theory_y];​​​​%Note:100数据进行ECDF会产生101个数的ECDF坐标，因此为了填充颜色，这里更改theory_y的维数​​


​​fill([X​​​​',fliplr(X'​​​​)],[theory_y​​​​',fliplr(F'​​​​)],​​​​'y'​​​​);​​


​​xlabel(​​​​'x'​​​​);​​


​​ylabel(​​​​'F(x)'​​​​);​​



好了，到现在开始要说QQ-plot了。我们用如下两个例子来说明:QQ-plot for 平均分布 & QQ-plot for 正态分布。

QQ-plot for 平均分布：

Sample中有n个数据， x1,x2,...,xnx1,x2,...,xn 。我们首先对数据进行排序，使之满足 x1<x2<...<xnx1<x2<...<xn 。我们将x所在区间[0,1]进行n等分。即变为 [0,1n],(1n,2n],...,(n−1n,1][0,1n],(1n,2n],...,(n−1n,1]n个自区间。为了符合平均分布，我们期望第q个数据的值坐落在第q个子区间的中间值，也就是


Eq=q−0.5nEq=q−0.5n


现在我们可以理解Quantile-Quantile(q-q) Plot了，第1个Q是Data Sample的分位即 x1,x2,...,xnx1,x2,...,xn ；第2个Q便是期望 EqEq 。因此QQ-plot其实就是n个点的集合


(q−0.5n,xq),�forq=1,2,...,n(q−0.5n,xq),�forq=1,2,...,n


因此在QQ-plot for平均分布中，当QQ点越接近y=x时，那么数据越接近平均分布。下面我们考虑表1中的5组数据，以及随机生成50个，500个，1000个数据的QQ-plot图，如图3所示。可以看出，当sample size越大，QQ-plot越接近y=x。

表1. Computing the Expected Quantile Values.

Data (x)

Rank (q)

Middle of the qth Interval

0.03

0.24

0.41

0.59

0.67

1

2

3

4

5

0.1

0.3

0.5

0.7

0.9

​​​​

图3. QQplot for uniform data

相关代码：




1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16




​​A=[0.03,0.24,0.41,0.59,0.67];​​


​​B=unifrnd(0,1,1,50);​​


​​C=unifrnd(0,1,1,500);​​


​​D=unifrnd(0,1,1,1000);​​


​​subplot(2,2,1);​​


​​gqqplot(A,​​​​'unif'​​​​);​​


​​title(​​​​'Sample size n = 5'​​​​);​​


​​subplot(2,2,2);​​


​​gqqplot(B,​​​​'unif'​​​​);​​


​​title(​​​​'Sample size n = 50'​​​​);​​


​​subplot(2,2,3);​​


​​gqqplot(C,​​​​'unif'​​​​);​​


​​title(​​​​'Sample size n = 500'​​​​);​​


​​subplot(2,2,4);​​


​​gqqplot(D,​​​​'unif'​​​​);​​


​​title(​​​​'Sample size n = 1000'​​​​);​​



QQ-plot for 正态分布：

这个就简单了，跟上面是一样的道理。我们取Z为标准的正态分布， μ=0,σ=1μ=0,σ=1 。现将n个数据进行排序，再做出相应的QQ-plot点的集合


(Φ−1(q−0.5n),xq),�forq=1,2,...,n(Φ−1(q−0.5n),xq),�forq=1,2,...,n


同样我们给出了表2,5组正态分布的数据以及其相应的期望值。为了比较，我们也随机产生了n为50,500,1000的正态分布随机数进行QQ-plot，如图4所示。

表2. Computing the expected quantile values for normal data.

Data (z)

Rank (q)

Middle of theqth Interval

Normal(q)

-1.96

-.78

.31

1.15

1.62

1

2

3

4

5

0.1

0.3

0.5

0.7

0.9

-1.28

-0.52

0.00

0.52

1.28

​​​​

图4. QQplot for normal data

代码如下：




1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16




​​E=[-1.96,-0.78,0.31,1.15,1.62];​​


​​F=randn(1,50);​​


​​G=randn(1,500);​​


​​H=randn(1,1000);​​


​​subplot(2,2,1);​​


​​gqqplot(E,​​​​'norm'​​​​);​​


​​title(​​​​'Sample size n = 5'​​​​);​​


​​subplot(2,2,2);​​


​​gqqplot(F,​​​​'norm'​​​​);​​


​​title(​​​​'Sample size n = 50'​​​​);​​


​​subplot(2,2,3);​​


​​gqqplot(G,​​​​'norm'​​​​);​​


​​title(​​​​'Sample size n = 500'​​​​);​​


​​subplot(2,2,4);​​


​​gqqplot(H,​​​​'norm'​​​​);​​


​​title(​​​​'Sample size n = 1000'​​​​);​​



好吧，到此为止，讲的差不多了，其实不难的。最后我们再来看一遍维基百科上对QQ-plot的定义：

“在统计学中，QQ-plot(Q代表分位数Quantile)是一种通过画出分位数来比较两个概率分布的图形方法。首先选定区间长度，点(x,y)对应于第一个分布(x轴)的分位数和第二个分布(y轴)相同的分位数。因此画出的是一条含参数的曲线，参数为区间个数。如果被比较的两个分布比较相似，则其QQ图近似地位于y=x上。如果两个分布线性相关，则QQ图上的点近似地落在一条直线上，但并不一定是y=x这条线。QQ图同样可以用来估计一个分布的位置参数。”

Sigh!应该理解了... 关于gqqplot函数的使用，请参考 ​​http://ackjack.com/?p=56​​


​​ Posted by sinknode ​ ​  ​​Math & Matlab​ ​   ​​matlab​​, ​​QQ-plot​​ ​ ​ 1 Comment​ ​


​​Matlab中实现QQ-plot的一个好工具gqqplot »​​