8.MATLAB參數估計與假設檢驗-非參數參數檢驗-分布的擬合與檢驗__MATLAB資料分析與統計

分布的擬合與檢驗

更多MATLAB資料分析視頻請點擊，或者在網易雲課堂上搜尋《MATLAB資料分析與統計》 http://study.163.com/course/courseMain.htm?courseId=1003615016

 在某些統計推斷中， 通常假定總體服從一定的分布（例如常態分佈），然後在這個分布的基礎上，構造相應的統計量，根據統計量的分布做出一些統計推斷，而統計量的分布通常依賴於總體的分布假設，也就是說總體所服從的分布在統計推斷中至關重要，會影響到結果的可靠性。從這個意義上來說，由樣本觀測資料去推斷總體所服從的分布是非常必要的。這一節的就是根據樣本觀測資料擬合總體的分布，並進行分布的檢驗。

  下面以

成績為資料，推斷總成績資料所服從的分布。

在根據樣本觀測資料對總體所服從的分布做推斷時，通常需要藉助於描述性統計量和統計圖，首先從直觀上對分布形式作出判斷，然後再進行檢驗，在前面的分析中，我們已經知道總成績的資料近似服從常態分佈，下面將調用MATLAB函數（chi2gof、jbtest、kstest和lillietest）進行檢驗。

（1）卡方擬合優度檢驗

  chi2gof 函數用來做分布的卡方擬合優度檢驗，檢驗樣本是否服從指定的分布。chi2gof函數的原理是：它用若干個小區間把樣本觀測資料進行分組（預設情況下分成10個小組），使得理論上每組包含5個以上的觀測，即每組的理論頻數大於或等於5，若不滿足這個要求，可以通過合并相鄰的組來達到這個要求。

 chi2gof函數的調用格式：

  <1>h=chi2gof(x)

   進行卡方擬合優度檢驗，檢驗樣本x是否服從常態分佈，原假設樣本x服從常態分佈，其中分布參數由x進行估計。輸出參數h等於0或1，若h=0，則在顯著性水平0.05下接受原假設，認為x服從常態分佈；若為1，則在顯著性水平0.05下拒絕原假設。

 <2> [h,p]=chi2gof(.....)

   返回檢驗的p值，當p值小於或等於顯著性水平時，拒絕原假設，否則接受原假設。

 <3> [h,p,stats]=chi2gof(.....)

   返回一個結構體變數stats，它包含以下欄位：

  chi2stat:卡方檢驗統計量

  df ： 自由度

  edges：合并後各區間的邊界向量

  O： 落入每個小區間內觀測的個數，即時間頻數

 E： 每個小區間對應的理論頻數

%讀取檔案成績.xls的第一個工作表中的G2：G52中的資料，即總成績資料
score=xlsread('成績.xls','G2:G52');
%去掉總成績中的0，即缺考資料
score=score(score>0);

%進行卡方擬合優度檢驗
[h,p,stats]=chi2gof(score)

h =

     1


p =

    0.0244


stats =

    chi2stat: 9.4038
          df: 3
       edges: [49.0000 68.6000 73.5000 78.4000 83.3000 88.2000 98.0000]
           O: [4 10 6 15 4 10]
           E: [7.4844 6.9183 8.9423 9.1961 7.5245 8.9344]

由於傳回值h=1，p=0.0244<0.05,在顯著性水平0.05下，可以認為總成績不服從常態分佈。

<4> [......]=chi2gof(x,name1,val1,name2,val2,.......)

   通過可選的成對出現的參數名和參數值來控制初始分組、原假設中的分布、顯著性水平等。控制初始分組的參數與參數值如下表：

  參數名             參數值                                說明

‘nbins’             正整數，預設值為10      分組（或區間）的個數

‘ctrs’                 向量                                   指定各區間的中點

‘edges’            向量                                   指定各個區間的邊界

注意：上述三個參數不能同時指定，一次調用最多隻能指定其中的一個參數，因為後兩個參數已經潛在的指定了分組數了

%讀取檔案成績.xls的第一個工作表中的G2：G52中的資料，即總成績資料
score=xlsread('成績.xls','G2:G52');
%去掉總成績中的0，即缺考資料
score=score(score>0);

%指定各初始小區間的中點
ctrs=[50,60,70,78,85,94];
%指定ctr參數，進行卡方擬合優度檢驗
[h,p,stats]=chi2gof(score,'ctrs',ctrs)

h =

     0


p =

    0.3747


stats =

    chi2stat: 0.7879
          df: 1
       edges: [45.0000 74.0000 81.5000 89.5000 98.5000]
           O: [15 16 10 8]
           E: [15.2451 14.0220 12.3619 7.3710]

通過‘ctrs’參數控制初始分組數為6，利用ctrs向量指定了初始6個初始小區間的中點。檢驗結果表明通過合并相鄰區間，初始的6個小區間最終被合并成4個小區間。返回的h=0，p=0.3747>0.05,認為總成績服從常態分佈，該常態分佈的均值可由mean（score）計算出，標準差可由std（score）計算出。

下面通過‘nbins’參數控制初始分組數為6

%讀取檔案成績.xls的第一個工作表中的G2：G52中的資料，即總成績資料
score=xlsread('成績.xls','G2:G52');
%去掉總成績中的0，即缺考資料
score=score(score>0);

%指定‘nbins’參數進行卡方擬合優度檢驗
[h,p,stats]=chi2gof(score,'nbins',6)

h =

     0


p =

    0.3580


stats =

    chi2stat: 0.8449
          df: 1
       edges: [49.0000 73.5000 81.6667 89.8333 98.0000]
           O: [14 17 10 8]
           E: [14.4027 15.1752 12.4207 7.0014]

以上兩次調用得到的h值相同，但是p值和stats不相同，並且與chi2gof函數第一次調用的時候檢驗結論相反，這說明了卡方擬合優度檢驗對分組結果比較敏感，在使用chi2gof函數時，應使得每個分組（小區間）所包含的觀測數均在5個以上。

chi2gof函數可以利用以下參數值來控制原假設中的分布

參數名                                    參數值                                                          說明

‘cdf’                   函數名字字串、函數控制代碼、                         指定原假設中的分布，與‘expected’參數

                        由函數字串（或函數控制代碼）與函數中          不同同時出現，若為函數名字串或函數控制代碼，則x是函數的唯一

                        所含參數的參數值構成的元胞數組                 輸入參數；若是由函數名字字串（或函數控制代碼）與函數中所含有參數值

                                                                                                      構成的元胞數組，則x是函數的第一個輸入參數，其他參數為後續輸入

‘expected’                         向量                                              指定各區間的理論頻數，與‘cdf’不能同時出現

‘nparams’                      向量                                                 指定分布中待估參數的個數，它確定了卡方分布的自由度

chi2gof函數控制檢驗的其他方面的參數與參數值列表

參數名             參數值                                說明

‘emin’      非負整數，預設值5       指定一個區間對應的最小理論頻數，初始分組中，

                                                            理論頻數小於這個值的區間和相鄰區間合并。如果指定為0，

                                                            將不進行區間合并

‘frequency’   與x等長的向量               指定x中各元素出現的頻數

‘alpha’           0--1之間的數，預設值0.05        指定檢驗的顯著性水平

%讀取檔案成績.xls的第一個工作表中的G2：G52中的資料，即總成績資料
score=xlsread('成績.xls','G2:G52');
%去掉總成績中的0，即缺考資料
score=score(score>0);

%求平均值ms和標準差ss
ms=mean(score);
ss=std(score);

%參數'cdf'的值是由函數名字串與函數中所包含參數的參數值構成的元胞數組
[h1,p1,stats1]=chi2gof(score,'nbins',6,'cdf',{'normcdf',ms,ss})

%參數‘cdf’的值有函數控制代碼與函數中所包含的參數值構成的元胞數組
[h2,p2,stats2]=chi2gof(score,'nbins',6,'cdf',{@normcdf,ms,ss})

%指定初始分組數為6，檢驗總成績資料是否服從參數為ms的泊松分布
[h3,p3,stats3]=chi2gof(score,'nbins',6,'cdf',{@poisscdf,ms})

h1 =

     0


p1 =

    0.3580


stats1 =

    chi2stat: 0.8449
          df: 1
       edges: [49.0000 73.5000 81.6667 89.8333 98.0000]
           O: [14 17 10 8]
           E: [14.4027 15.1752 12.4207 7.0014]


h2 =

     0


p2 =

    0.3580


stats2 =

    chi2stat: 0.8449
          df: 1
       edges: [49.0000 73.5000 81.6667 89.8333 98.0000]
           O: [14 17 10 8]
           E: [14.4027 15.1752 12.4207 7.0014]


h3 =

     0


p3 =

    0.4871


stats3 =

    chi2stat: 1.4385
          df: 2
       edges: [49.0000 73.5000 81.6667 89.8333 98.0000]
           O: [14 17 10 8]
           E: [13.3213 16.9281 12.8698 5.8808]

從檢驗結果看出，在顯著性水平=0.05下，前兩個輸出顯示了檢驗資料服從N（ms，ss）的常態分佈，最後一個輸入顯示了檢驗資料服從參數為ms的泊松分布。    

故綜上各檢驗結果，在顯著性水平0.05下，扔認為成績資料服從常態分佈，其均值為ms，標準差為ss。

（2）  Jarque-Bera檢驗

 jbtest函數用來做Jarque-Bera檢驗，檢驗樣本是否服從常態分佈，調用該函數時不需要指定分布的均值和方差。由於常態分佈的偏度為0，峰值為3，若樣本服從常態分佈，則樣本偏度應接近於0，樣本峰度接近於3，基於此，Jarque-Bera檢驗是利用樣本偏度和峰度構造檢驗統計量。

 jbtest函數的調用格式如下：

 <1> h=jbtest（x）

    檢驗樣本x是否服從均值和方差未知的常態分佈，原假設是x服從均值常態分佈。當輸出h=1時，表示樣本在顯著性水平=0.05下拒絕原假設；當h=0時，則在顯著性水平=0.05下接受原假設。jbtest函數會把x中的NaN（不明確的數值）作為缺失資料忽略。

 <2> h=jbtest（x，alpha）

     指定顯著性水平alpha進行分布的檢驗，原假設和備擇假設同時

 <3> [h,p]=jbtest（......）

  返回檢驗的p值，當p小於或等於給定的顯著性水平alpha時，拒絕原假設，大於顯著性水平時，接受原假設

 <4>[h,p,jbstat]=jbtest（.......）

    返回檢驗統計量的觀測值jbstat

  <5>[h,p,jbstat,critval]=jbtest（.....）

      返回檢驗的臨界值critval。當jbstat>=crival時，在顯著性水平alpha下拒絕原假設

 <6> [h,p,......]=jbtest（x，alpha，mctol）

  指定一個終止容限mctol，利用蒙特卡洛類比方法計算p值的近似值

注意：jbtest函數只是基於樣本偏度和峰度進行正態性檢驗，結果受異常值的影響比較大，可能會出現比較大的偏差。

例：

randn('seed',0);  %指定隨機數產生器的初始種子為0
x=randn(10000,1); %產生10000個服從標準常態分佈的隨機數
h=jbtest(x)       %調用jbtest函數進行正態性檢驗

x(end)=5;
h1=jbtest(x)

h =

     0


h1 =

     1

從以上結果可以發現，對於一個包含10000個元素的標準常態分佈的隨機數向量，指改變它的最後一個元素，就會導致檢驗的結論完全相反，這就充分說明了jbtest函數的局限性，受異常值的影響比較大。

下面調用jbtest函數對成績資料進行正態性檢驗

%讀取檔案成績.xls的第一個工作表中的G2：G52中的資料，即總成績資料
score=xlsread('成績.xls','G2:G52');
%去掉總成績中的0，即缺考資料
score=score(score>0);

[h,p]=jbtest(score)

h =

     1


p =

    0.0193

由於傳回值h=1，p<0.05,所以在顯著性水平0.05下拒絕了原假設，認為總成績資料不服從常態分佈。不過由於jbtest函數的局限性，這個結論僅作為參考，還應結合其他函數的檢驗結果做出綜合的推斷。

（3）單樣本的Kolmogorov-Smirnov（K-S）檢驗

   kstest函數用來做單樣本的K-S檢驗：它可以作雙側檢驗，檢驗樣本是否服從指定的分布；也可做單側檢驗，檢驗樣本的分布函數是否在指定的分布函數之上或之下，這裡的分布是完全確定的，不含有未知參數。kstest函數根據樣本的經驗分布函數Fn（x）和指定的分布函數G（x）構造檢驗統計量

       KS=max（|Fn（x）-G（x）|）

  kstest函數的調用格式如下：

 <1> h=kstest（x）

   檢驗樣本x是否服從標準常態分佈，原假設是x服從標準常態分佈，備擇假設是x不服從標準常態分佈。當輸出h=1時，在顯著性水平=0.05下拒絕原假設；當h=0時，則在顯著性水平=0.05下接受原假設

 <2>h=kstest（x，CDF）

    檢驗樣本x是否服從由CDF定義的連續分布，這裡的CDF可以是包含兩列元素的矩陣，也可以是也可是是機率分布對象。當CDF是包含兩列元素的矩陣時，它的第一列表示隨機變數的可能取值，可以是樣本x中的值，也可以不是。CDF的第二列是指定分布函數的取值，如果CDF為空白，則建議樣本x是否服從標準的常態分佈。

<3> h=kstest（x，CDF，alpha）

 指定檢驗的顯著性水平alpha，預設值是0.05

 <4>h=kstest（x，CDF，alpha，type）

  用type參數指定檢驗的類型（雙側或單側）。type參數的可能取值為

  ‘unequal’：雙側檢驗，備擇假設是總體分布函數不等於指定的分布函數

   'larger':單側檢驗，備擇假設是總體分布函數大於指定的分布函數

    ‘smaller’：單側檢驗，備擇假設是總體分布函數小於指定的分布函數

 <5> [h,p,ksstat,cv]=kstest(......)

      返回檢驗的p值、檢驗統計量和觀測值ksstat和鄰界值cv

例：調用kstest函數檢驗總成績是否服從常態分佈

%讀取檔案成績.xls的第一個工作表中的G2：G52中的資料，即總成績資料
score=xlsread('成績.xls','G2:G52');
%去掉總成績中的0，即缺考資料
score=score(score>0);

%首先調用kstest檢驗是否服從標準的常態分佈
h=kstest(score)

%然後檢驗是否服從均值為ms，標準差為ss的常態分佈
ms=mean(score);
ss=std(score);
%產生cdf矩陣，用來指定分布：均值為ms，標準差ss的常態分佈
%cdf的第二列是指定分布函數的取值，由累積常態分佈函數normcdf來確定
cdf=[score,normcdf(score,ms,ss)];

%調用kstest函數，檢驗總成績是否服從由cdf指定的分布
[h1,p,ksstat,cv]=kstest(score,cdf)

h =

     1


h1 =

     0


p =

    0.5486


ksstat =

    0.1107


cv =

    0.1903

由於h=1，所以在顯著性水平0.05下，拒絕原假設（x服從標準的常態分佈）；h1=0，所以在顯著性水平=0.05下，接受原假設，認為總成績資料服從均值為ms、標準差為ss的常態分佈。

（4）雙樣本的K-S檢驗

 kstest2函數用來做兩個樣本的K-S檢驗，它可以作雙側檢驗，檢驗兩樣本是否服從相同的分布，也可以作單側檢驗，檢驗一個樣本的分布函數是否在另一個樣本的分布函數之上或之下，這裡的分布函數是完全確定的，不含未知數。kstest2函數對比兩樣本的經驗分布函數，構造檢驗統計量

   KS=max（|F1（x）-F2（x）|）

 其中F1（x）和F2（x）分別為兩樣本的經驗分布函數。

 kstest2函數的調用格式如下：

  <1>h=kstest2（x1，x2）

       檢驗樣本x1和x2是否具有相同的分布，原假設是x1與x2來自相同的連續分布，備擇假設是來自不同的連續分布。當輸出h=1時，在顯著性水平=0。05下拒絕原假設；當h=0時，則在顯著性水平=0.05下接受原假設。這裡並不要求x1與x2具有相同的長度。

   <2>h=kstest2（x1，x2，alpha）

    指定檢驗的顯著性水平alpha，預設是0.05

   <3>h=kstest2（x1，x2，alpha，type）

    用type參數指定檢驗的類型（雙側或單側）。type參數的可能取值為

   ‘unequal’：雙側檢驗，備擇假設是兩個總體的分布函數不相等

   ‘larger’ ：單側檢驗，備擇假設是第1個總體的分布函數大於第二個總體的分布函數

   ‘smaller’：單側檢驗，備擇假設是第1個總體的分布函數小於第二個總體的分布函數

   <4>[h,p]=kstest2（.....）

    返回檢驗的漸近p值，當p小於或等於顯著性水平alpha時，拒絕原假設。樣本容量越大，p值越精確，同樣要求 

                 (n1*n2)/(n1+n2)>=4

    其中，n1，n2分別為樣本x1和x2的樣本容量。

    <5>[h,p,ks2stat]=kstest2(......)

         返回檢驗統計量的觀測值ks2stat

例：

 對成績.xls中的資料，調用kstest2函數檢驗班級60101和班級60102兩個班級的總成績是否服從相同的分布。

 %讀取檔案成績.xls的第一個工作表中的班級資料即B2：B52
banji=xlsread('成績.xls','B2:B52');
%讀取檔案中的第一個工作表中的總成績資料，即G2：G52
score=xlsread('成績.xls','G2:G52');

%去除缺考成績
banji=banji(score>0);
score=score(score>0);

%分別提取60101和60102班的總成績
score1=score(banji==60101);
score2=score(banji==60102);

%調用kstest2函數檢驗兩個班級的總成績是否服從相同的分布
[h,p,ks2stat]=kstest2(score1,score2)

h =

     0


p =

    0.7597


ks2stat =

    0.1839

由於h=0，p>0.05,所以在顯著性水平0.05下，接受原假設，認為兩個班級的總成績服從相同的分布。

（5）Lillefors檢驗

 當總體均值和方差未知時，Lilliefor（1967年）提出了用樣本均值和標準差代替總體的均值和標準差，然後使用K-S檢驗，這就是所謂的Lilliefors檢驗。

 lillietest函數用來做Lilliefors檢驗，檢驗樣本是否服從指定的分布，這裡分布的參數都是未知的，根據樣本資料做出估計，可用的分布有常態分佈，指數分布，和極值分布。

 lilltest函數的調用格式如下：

 <1> h=lillietest（x）

      檢驗樣本x是否服從均值和方差未知的常態分佈，原假設是x服從常態分佈。當輸出h=1時，表示在顯著性水平=0.05下拒絕原假設；當h=0時，則在顯著性水平=0.05下接受原假設，lillietest函數會把x中的NaN（不明確的數值）作為缺失資料忽略

 <2>h=lillietest（x，alpha）

   指定顯著性水平alpha進行分布的檢驗，原假設和備擇假設同上。

 <3> h=lillietest（x，alpha，distr）

   檢驗樣本x是否服從參數distr指定的分布，distr為字串變數，可能的取值為'norm'(常態分佈，預設情況)，‘exp’（指數分布），‘ev’（極值分布）

  <4>[h,p]=lillietest(......)

   返回檢驗的p值，當p小於或等於給定的顯著性水平alpha時，拒絕原假設，大於顯著性水平，接受原假設

 <5>[h,p,kstat]=lillietest(........)

   返回檢驗統計量的觀測值kstat

  <6>[h,p,kstat,critval]=lillietest(.....)

    返回檢驗的臨界值critval。當kstat>=critval時，在顯著性水平alpha下拒絕原假設

 <7>[h,p,....]=lillietest(x,alpha,distr,mctol)

  指定一個終止容限mctol，直接利用蒙特卡洛類比法計算p值

例：利用lillietest函數對總成績資料進行正態性檢驗

%讀取檔案中的第一個工作表中的總成績資料，即G2：G52
score=xlsread('成績.xls','G2:G52');
score=score(score>0);
%調用lillietest函數進行Lilliefors檢驗，檢驗總成績是否服從常態分佈
[h,p]=lillietest(score)

%調用lillietest函數檢驗總成績是否服從整數分布
[h1,p1]=lillietest(score,0.05,'exp')

h =

     0


p =

    0.1346


h1 =

     1


p1 =

   1.0000e-03

由於h=0，所以在顯著性水平0.05下接受原假設，認為總成績服從常態分佈，由於Lilliefors檢驗用樣本均值和樣本標準差代替總體均值和標準差，因此常態分佈的均值為mean（score），標準差為std（score）；h1=1，所以在顯著性水平0.05下，拒絕原假設，認為總成績資料不服從指數分布。

（6）最終結論

  通過前面分別用chi2gof、jbtest、kstest、kstest2、lillietest函數進行檢驗，只有jbtest函數的檢驗結果認為總成績資料不符合常態分佈，其他都符號常態分佈，但由於jbtest函數的局限性（結果受異常值影響較大），所以可以認為總成績資料服從常態分佈，其均值為mean（score），標準差為std（score）

更多MATLAB資料分析視頻請點擊，或者在網易雲課堂上搜尋《MATLAB資料分析與統計》 http://study.163.com/course/courseMain.htm?courseId=1003615016