MATLAB 資料分析方法（第2版）2.2 資料分布及其檢驗

2.2　資料分布及其檢驗

樣本資料的數字特征刻畫了資料的主要特征，而要對資料的總體情況作全面的了解，就必須研究資料的分布。上節中的資料直方圖與q-q圖等能直覺粗略地描述資料的分布，本節進一步研究如何判定資料是否服從正态分布的問題。若不服從正态分布，那麼又可能服從怎樣的分布？

2.2.1　一維資料的分布與檢驗

1.經驗分布函數

設來自總體x的容量為n的樣本x1,x2,…,xn，樣本的次序統計量為x(1),x(2),…,x(n)，對于任意實數x,定義函數

fn(x)=0,若x<x(1)

kn,若x(k)≤x<x(k+1)　（k=1,2,…,n-1）

1,若x≥x(n)（2.2.1）

稱fn(x)為經驗分布函數。

由定義可知，fn(x)表示事件{x≤x}在n次獨立重複試驗中的頻率。

1933年，格裡汶科（glivenko）證明了以下結果：對于任一實數x,當n→∞時fn(x)以機率1一緻收斂于分布函數f(x),即

plimn→∞sup-∞<x<∞fn(x)-f(x)=0=1

　　這一結論表明：對于任一實數x，當n充分大時

f(x)≈fn(x)（2.2.2）

是以可用經驗分布函數來近似代替f(x)，這一點也是由樣本推斷總體的最基本理論依據之一。

在matlab中，作經驗（累積）分布函數圖形指令cdfplot的調用格式為：

① cdfplot(x)；　　　　　　　　　　　　　　　%作樣本x（向量）的經驗（累積）分布函數圖形

② h=cdfplot(x)；% h表示曲線的環柄

③ \[h,stats\]=cdfplot(x)；%輸出stats表示樣本最小值、最大值、均值、中值與标準差

通常，将樣本的直方圖與經驗分布函數圖結合應用，對資料的分布作出推斷。

例2.2.1　生成服從标準正态分布的50個樣本點，作出樣本的經驗分布函數圖，并與理論分布函數Φ(x)比較。

解：編寫程式如下。

clear

x=normrnd (0,1,50,1);%生成服從标準正态分布的50個樣本點

\[h,stats\]=cdfplot(x);%繪制樣本的經驗分布函數圖

hold on

plot(-3:0.01:3, normcdf(-3:0.01:3,0,1), 'r')　　%繪制理論分布函數圖

legend('樣本經驗分布函數fn(x)', '理論分布函數Φ(x)' ,'location','northwest');

輸出結果:

h =

　　3.0013

stats =

　　　 min: -1.8740%樣本最小值

　　　 max: 1.6924%最大值

　　　 mean: 0.0565%平均值

　　　 median: 0.1032%中間值

　　　 std: 0.7559%樣本标準差

圖2-12表明樣本的經驗分布函數圖形與理論分布函數圖很相近。

圖2-12　n(0,1)分布函數圖及其50個樣本點的經驗分布函數圖

2.總體分布的正态性檢驗

進行參數估計和假設檢驗時，通常總是假定總體服從正态分布，雖然在許多情況下這個假定是合理的，但是當要以此為前提進行重要的參數估計或假設檢驗，或者人們對它有較大懷疑的時候，就确有必要對這個假設進行檢驗。進行總體正态性檢驗的方法有很多種，以下針對matlab統計工具箱中提供的程式，簡單介紹幾種方法。

（1）jarque-bera檢驗

jarque-bera檢驗簡稱jb檢驗，它是利用正态分布的偏度sk和峰度ku，構造一個包含sk、ku且自由度為2的卡方分布統計量

jb=n16j2+124b2~χ2(2)（2.2.3）

其中j=1n∑ni=1xi-xs3,b=1n∑ni=1xi-xs4-3。

對于顯著性水準α，當jb統計量小于χ2分布的1-α分位數χ21-α(2)時接受h0，即認為總體服從正态分布；否則拒絕h0，即認為總體不服從正态分布。這個檢驗适用于大樣本，當樣本容量n較小時需慎用。

在matlab中，jb檢驗指令jbtest的調用格式為：

h = jbtest(x,alpha)；

或

\[h,p,jbstat,cv\] = jbtest(x,alpha)；

對輸入向量x進行jarque-bera測試，顯著性水準alpha預設為0.05。輸出h為測試結果，若h=0，則不能拒絕x服從正态分布；若h=1，則可以否定x服從正态分布。輸出p為接受假設的機率值，p小于alpha，則可以拒絕是正态分布的原假設；jbstat為測試統計量的值，cv為是否拒絕原假設的臨界值，jbstat大于cv可以拒絕是正态分布的原假設。

指令jbtest一般用于大樣本，對于小樣本用指令lillietest。

（2）kolmogorov-smirnov檢驗

kolmogorov-smirnov檢驗簡稱ks檢驗，它是通過樣本的經驗分布函數與給定分布函數的比較，推斷該樣本是否來自給定分布函數的總體。設給定分布函數為g(x)，構造統計量

dn=maxn(fn(x)-g(x))（2.2.4）

即兩個分布函數之差的最大值，對于假設h0：總體服從給定的分布g(x)，及給定的α，根據dn的極限分布确定統計量關于是否接受h0的數量界限。

因為這個檢驗需要給定g(x)，是以當用于正态性檢驗時隻能做标準正态檢驗，即h0：總體服從标準正态分布n(0,1)。

在matlab中，ks檢驗指令kstest的調用格式為：

h = kstest(x)；

h = kstest(x,cdf)；

\[h,p,ksstat,cv\] = kstest(x,cdf,alpha)；

把向量x中的值與标準正态分布進行比較并傳回假設檢驗結果h。如果h=0表示不能拒絕原假設，即不能拒絕服從正态分布；若h=1，則可以否定x服從正态分布。假設的顯著水準預設值是0.05。cdf是一個兩列矩陣，矩陣的第一列包含可能的x值，第二列是假設累積分布函數g(x)的值。在可能的情況下，cdf的第一列應包含x中的值，如果第一列沒有，則用插值的方法近似。指定顯著水準alpha，傳回p值、ks檢驗統計量ksstat、截斷值cv。

（3）lilliefors檢驗

lilliefors檢驗是改進ks檢驗并用于一般的正态性檢驗，原假設h0：總體服從正态分布n(μ,σ2)，其中μ、σ2由樣本均值和方差估計。

該檢驗的matlab指令lillietest的調用格式為：

h = lillietest(x,alpha)；

 \[h,p,lstat,cv\] = lillietest(x,alpha)；

對輸入向量x進行lilliefors測試，顯著性水準alpha在0.01和0.2之間，預設時為0.05。輸出p為接受假設的機率值，lstat為測試統計量的值，cv為是否拒絕原假設的臨界值。h為測試結果，若h=0，則不能拒絕x服從正态分布；若h=1，則可以否定x服從正态分布。p小于alpha，則可以拒絕是正态分布的原假設；lstat大于cv可以拒絕是正态分布的原假設。

例2.2.2　在例2.1.5中，檢驗“中國銀行”的股票的收盤價是否服從正态分布。

解：程式如下。

a=xlsread('yhgspj.xls');　　　　　　%讀取收盤價資料

h1=jbtest(a(:,1))%jb檢驗

h2=kstest(a(:,1))%ks檢驗

h3=lillietest(a(:,1))%改進ks檢驗

qqplot(a(:,1))

程式運作結果：h1=1,h2=1,h3=1,三種檢驗都不支援“中國銀行”股票的收盤價服從正态分布，q-q圖（如圖2-13所示）也可看出收盤價不服從正态分布,進而可以認為收盤價不服從正态分布。

圖2-13　中國銀行收盤價的q-q圖

2.2.2　多元資料的正态分布檢驗

1.多元正态分布的概念與性質

設p維總體x=(x1,x2,…,xp)t的分布密度函數為

f(x1,x2,…,xp)=1(2π)p/2Σ1/2exp-12(x-μ)tΣ-1(x-μ)（2.2.5）

則稱x服從p維正态分布，記為x~np(μ,Σ)，其中

x=(x1,x2,…,xp)t,μ=(μ1,μ2,…,μp)t,Σ=σ11σ12…σ1p

σ21σ22…σ2p

…

σp1σp2…σpp

μ稱為總體均值向量，Σ稱為總體協方差矩陣。

多元正态分布具有如下性質：

1）多元正态分布的邊緣分布還是服從正态分布，但反之不真。

2）多元正态分布線上性變換下仍然服從多元正态分布。即若x~np(μ,Σ)，a為s×p階常數矩陣，d為s維常數向量，則

ax+d~ns(aμ+d,aΣat)

　　3）服從正态分布的随機變量間互相獨立與不相關等價。

2.多元正态分布的q-q圖檢驗方法

對于來自總體且由式(2.1.16)表示的樣本資料矩陣x，怎樣檢驗其是否是來自于多元正态總體呢？一般可按照以下q-q圖檢驗方法，具體的過程如下：

1）由樣本資料矩陣x計算均值向量x和協方差矩陣s。

2）計算樣本點x［t］到x的馬氏平方距離，其中x［t］表示x的第t行。（參見第4章）。

d2t=(x［t］-x)ts-1(x［t］-x)　(t=1,…,n)

　　3）對上述馬氏平方距離從小到大排序

d2(1)≤d2(2)≤…≤d2(n)

　　4）計算pt=t-0.5n(t=1,2,…,n)及χ2t，其中χ2t滿足h(χ2tp)=pt。

5）以馬氏距離為橫坐标，χ2t分位數為縱坐标作n個點(d2(t),χ2t)的平面散點圖，即分布的q-q圖。

6）考查散點圖是否在一條通過原點且斜率為1的直線上，若是則接受資料來自p維正态分布總體的假設，否則拒絕正态分布假設。

以上q-q圖檢驗方法的matlab程式實作如下。

%輸入樣本資料矩陣x

x=\[data\];

\[n,p\]=size(x);　　　　　　　　　　　　　　　　%x的行數及列數

d=mahal(x,x);%計算馬氏距離

d1=sort(d);%馬氏距離從小到大排序

pt=\[\[1:n\]-0.5\]/n;%計算分位數

x2=chi2inv(pt,p);%計算χ2t

plot(d1,x2','*',\[0:m\],\[0:m\],'-r')%作散點圖與直線y=x,其中m是正整數

以下舉例說明上述程式的應用。

例2.2.3　為了研究某種疾病，對一批60人分為三組：g1、g2、g3，同時進行4項名額的檢測，即β脂蛋白(x1)、甘油三酯(x2)、α脂蛋白(x3)、前β脂蛋白(x4)。檢測的結果列在表2-7中，現将三組檢驗資料視為一個總體，問總體是否服從四維正态分布？

表2-7　四項名額檢測資料

g1

g2

g3

x1

x2

x3

x4

260

75

40

18

310

122

30

21

320

64

39

17

200

72

34

60

35

59

37

11

240

87

45

190

27

15

360

88

28

26

170

65

225

16

295

100

36

12

270

110

24

32

205

130

23

210

82

31

380

114

69

280

67

55

42

10

46

38

20

250

117

29

107

76

73

33

125

94

103

14

330

112

13

345

127

62

22

(續)

81

120

19

119

25

57

8

370

70

135

　　

資料來源：高惠璇，應用多元統計分析，北京大學出版社，2005年第一版。

解：先将表2-7中資料按原位置作為矩陣a輸入，然後整理成樣本資料矩陣x。程式如下。

a=\[260　　75　　40　　18　　310　　122　　30　　21　　320　　64　　39　　17;

　200　　72　　34　　17　　310　　60　　35　　18　　260　　59　　37　　11;

　……

　260　　135　　39　　29　　280　　40　　37　　17　　250　　117　　36　　16\];

x=\[a(:,1:4);a(:,5:8);a(:,9:12)\];%整理成樣本資料矩陣x

\[n,p\]=size(x);

d1=sort(d);%從小到大排序

plot(d1,x2','*',\[0:12\],\[0:12\],'-r')%作圖

輸出圖形如圖2-14所示。

　圖2-14　四項檢測資料的多元

正态檢驗q-q圖

從圖2-14可以看出，資料點基本落在直線上，故不能拒絕該資料服從四維正态分布的假設。

3.多個總體協方差矩陣的相等性檢驗

（1）兩個總體協方差矩陣相等的檢驗

設從兩個總體分别抽取樣本容量為n1、n2的兩個樣本，其樣本的協方差矩陣分别為s1、s2，那麼在兩個總體協方差矩陣相等時，其總體的協方差矩陣的估計為

s=(n1-1)s1+(n2-1)s2(n1+n2-2)

若檢驗兩個總體的協方差矩陣相等，可以有如下假設檢驗：

h0:si=sh1:si≠s,(i=1,2)

檢驗統計量

qi=(ni-1)［lns-lnsi-p+tr(s-1si)］~χ2(p(p+1)/2)(i=1,2)（2.2.6）

其中·表示行列式，p是向量的維數，tr表示矩陣的迹。

對給定的α，卡方分布臨界值為λ，若qi<λ(i=1,2)則接受h0，否則拒絕h0。

例2.2.4（1989年國際數學競賽a題：蠓的分類）　蠓是一種昆蟲，分為很多類型，其中有一種名為af，是能傳播花粉的益蟲；另一種名為apf，是會傳播疾病的害蟲，這兩種類型的蠓在形态上十分相似，很難差別。現測得6隻apf和9隻af蠓蟲的觸角長度和翅膀長度資料

apf：(1.14,1.78),(1.18,1.96),(1.20,1.86),(1.26,2.00),(1.28,2.00),(1.30,1.96)；

af：(1.24,1.72)，(1.36,1.74)，(1.38,1.64)，(1.38,1.82)，(1.38,1.90)，(1.40,1.70)，(1.48,1.82),(1.54,1.82),(1.56,2.08)。

判别apf與af兩類蠓蟲的協方差矩陣是否相等。

解：編寫檢驗協方差矩陣相等的源程式如下。

apf=\[1.14,1.78; 1.18,1.96;1.20,1.86;1.26,2.;1.28,2;1.30,1.96\];

af=\[1.24,1.72;1.36,1.74;1.38,1.64;1.38,1.82;1.38,1.90;1.40,1.70;1.48,1.82;1.54,1.82;1.56,2.08\];

n1=6;n2=9;p=2;

s1=cov(apf);s2=cov(af);

s=((n1-1)*s1+(n2-1)*s2)/(n1+n2-2);　　　　　　　　　　　　　　　%計算混合樣本方差

q1=(n1-1)*(log(det(s))-log(det(s1))-p+trace(inv(s)*s1)) ;

q2=(n2-1)*(log(det(s))-log(det(s2))-p+trace(inv(s)*s2)) ;%計算檢驗統計量觀測值

輸出結果為

q1 =

　　2.5784

q2 =

　　0.7418

給定α=0.05，查表得到臨界值χ21-α(3)=7.8147（指令chi2inv(0.95,3)），由于q1=2.5784<7.8147，q2=0.7418<7.8147，故認為兩類總體協方差矩陣相同。

（2）多個總體協方差矩陣相等的檢驗

設有k個p維總體gi,i=1,2,…,k，從每個總體中分别抽取樣本容量為ni(i=1,2,…,k)的k個樣本，其樣本的協方差矩陣為s1,s2,…,sk，用s1,s2,…,sk估計Σ1,Σ2,…,Σk。其中Σi為總體gi的協方差矩陣。

原假設　　h0:Σ1=Σ2=…=Σk；

備擇假設　h1:Σ1,Σ2,…,Σk至少有一對不相等。

在h0成立時，統計量

ξ=(1-d)m~χ2(f)（2.2.7）

其中m=(n-k)lns-∑ki=1(ni-1)lnsi，s=∑ki=1(ni-1)si/(n-k)，f=p(p+1)(k-1)/2為自由度，n=n1+n2+…+nk，

d=2p2+3p-16(p+1)(k-1)∑ki=11ni-1-1n-k,ni不全等

(2p2+3p-1）(k+1)6(p+1)(n-k),ni全等（2.2.8）

　　對給定的α，計算機率p=p(ξ>χ2α(f))，若p<α則拒絕h0，否則接受h0。

以上過程和程式可用下例說明。

例2.2.5　檢驗表2-7中三個總體g1、g2、g3的協方差矩陣是否相等（α=0.1）。

a=\[data\];　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　%輸入樣本資料

g1=a(:,1:4);%提取總體1的樣本

g2=a(:,5:8);

g3=a(:,9:12);

n=size(g1,1)+ size(g2,1)+ size(g2,1);%計算總的樣本容量

\[n1,p\]= size(g1);

k=3;

f=p*(p+1)*(k-1)/2;%統計量自由度

d=(2*p^2+3*p-1)*(k+1)/(6*(p+1)*(n-k));%由式（2.2.8）計算

s1=cov(g1);%協方差矩陣

s2=cov(g2);%協方差矩陣

s3=cov(g3);%協方差矩陣

s=(n1-1)*(s1+s2+s3)/(n-k);%總體協方差矩陣估計

m=(n-k)*log(det(s))-19*(log(det(s1))+log(det(s2))+log(det(s3)));%計算式（2.2.7）中的m值

t=(1-d)*m;%計算式（2.2. 7）統計量

p0=1-chi2cdf(t,f);%卡方分布機率

輸出結果：

t= 20.3316，p0= 0.4374

由于由統計量計算得到的機率為p0=0.4374>0.1，故判定3個總體協方差矩陣相等。