第十一 章 統計與迴歸

11.1 前言

統計的技巧與資料分析常常形影不離。一般統計使用加法、累加法、平均值，中間值等等，由於處理的對象是矩陣資料，故其基本統計之技巧已經廣為應用，其觀念也會在正常之運作中出現。統計學中比較特殊應用者為機率、亂數、常態分配等，而配合應用者為其相關之圖表。

在MATLAB中，有一個統計學工具箱，內藏各種統計學上需要應用的指令，可以執行上述與統計學有關之內容。這些相關的指令大部份以M-檔案組成，所以可利用type 這個功能檢視其內容。甚至可以更改其檔案名稱與內容，增加自己需要的功能，使其成為新的指令。此外，有些指令尚搭配繪圖介面，因而可以在繪圖模式下，進行資料與圖之適配，形成具體的方程式或實驗式，以供未來研究者使用。

統計工具箱中提供約二百餘個指令檔案，其中對機率分配方面則提供廿餘種機率型態，每種均有其相關的函數，諸如：

• 機率密度函數(pdf)
• 累積分佈函數(cdf)
• 累積分佈函數之反函數
• 亂數產生器
• 均數與變異數

11.2.1 常用統計指令--平均值MEAN

11.2.1平均值MEAN

統計學上對資料處理常用趨中的處理，求取均值或中間值等，均會取中的特徵。求取一個矩陣或向量之平均值時可用指令MEAN，其格式如下：

  M=mean(A,dim)

若A為向量，其結果Ｍ為單一值，亦即向量中各元素之平均；若Ａ為矩陣，則結果Ｍ為一列向量，其中元素為各行之平均值。dim為方向性參數，其預設值為1，表示結果係行向平均，故Ｍ為列向量；若dim=2，則係沿列向平均，結果Ｍ為行向量。例如：

A = [1 2 3; 3 3 6; 4 6 8; 4 7 7] 
A =
    1     2     3
    3     3     6
    4     6     8
    4     7     7 

M1=mean(A), M2=mean(A,2) 

M1 =
   3.0000    4.5000    6.0000

M2 =
    2
    4
    6
    6 

11.2.2 常用統計指令--中間值median

11.2.2中間值median

中間值亦可利用平均值的指令型式求得，其正式指令名稱為median。但其求得之值若非正好中間值，則會以接近中間值之兩值加以平均，其結果與mean之平均值仍有不同，下面以前述之Ｂ矩陣為例，比較median 與mean兩者執行後之不同處。

M1=median(B),M2=mean(B)  

M1 =
    5.5000    5.5000    5.0000    6.0000
M2 =
    5.0000    5.5000    4.7500    6.0000  

其他的平均值則有：
幾何平均(geomean)：各元素乘積再開總數之次方。中間有零值時，其結果為零。
調諧平均(harmmean)：各元素倒數和之倒數乘以總數。
修剪平均(trimmean)：去頭去尾再平均的方式，其頭尾部份為第二參數(％)之一半比例。

下面的例子為這些平均值的比較：

X=1:2:20
locate=[geomean(X) harmmean(X) mean(X) median(X) trimmean(X,25)]  

X =
     1     3     5     7     9    11    13    15    17    19
locate =
    7.6139    4.6877   10.0000   10.0000   10.0000  

由上述這些值看來，有些的結果相同，有些則因其分離的情況而有異。中間值及修剪平均值是不管偏異的大小，一般為避免有極端的數值的影響，採用這種平均法。

11.2.3 最大值與最小值 MIN & MAX

11.2.3最大值與最小值 MIN & MAX

要求一矩陣或向量中元素之最大與最小值時，其指令之型式如下：

  C = max|min(A)
  C = max|min(A,B)
  C = max|min(A,[],dim)
  [C,I] = max(...)

若A為向量，其結果C為單一值，亦即向量中各元素之最大或最小；若Ａ為矩陣，則結果Ｃ為一列向量，其中元素為各行之最大或最小。dim為方向性參數，其預設值為1，表示結果係行向取得最大或最小，故Ｃ為列向量；若dim=2，則係沿列向操作，結果Ｍ為行向量。注意要dim之參數時，需加在第三位置。此外，在輸出項中，Ｉ表示最大或最小元素之位置，不過此項功能僅求最大值時適用。例如：

A= [7 2 3 4; 2 4 5 6; 4 6 8 5; 6 7 6 1];

[C,Index]=max(A)  

A =
     7     2     3     4
     2     4     5     6
     4     6     8     5
     6     7     6     1
C =
     7     7     8     6
Index =
     1     4     3     2  

在這個指令中，比較特殊的是兩個矩陣Ａ與Ｂ之最大或最小比較，其結果Ｃ應為與Ａ或Ｂ相同的矩陣，但比較Ａ與Ｂ中對應元素間之最大與最小。例如：

B=[1 8 7 6;4 6 2 8;7 5 6 4;8 3 4 6]
[Cab]=min(A,B)  

B =
     1     8     7     6
     4     6     2     8
     7     5     6     4
     8     3     4     6
Cab =
     1     2     3     4
     2     4     2     6
     4     5     6     4
     6     3     4     1  

11.2.4 常用統計指令--變方值VAR

11.2.4變方值VAR

變方值為各樣本品與平均值差之平方和。又稱為變值常態檢定公式，其指令型式為var(X,1)。其計算之公式如下：

VAR(X)=[(X1-r)²+(X2-r)²+...+(Xn-r)²]/(n-1)

式中，n為其總項目，若n>1，則採用n-1以獲得無偏差之變方值。相關指令如下：

　　Y=var(X)
　　Y=var(X,1)
　　Y=var(X, w)
　　Y=var(X, w, dim)

Y=var(X)為執行上式計算之無偏差變方，Y=var(X,1)則採用n為除數。另w則為加權值向量，此權值必須為正值，且長度與Ｘ相同。

若Ｘ為矩陣時，通常預設為行向計算，但可以利用dim=2參數改為以列向為計算基礎，其結果為行向量。var指令會將其元素除以總和，因此權值總和為一。若w值為零，其結果如var(X)；若為１則如var(X,1)。

範例

>> x=1:9
x =

    1     2     3     4     5     6     7     8     9

>> var(x) %除以(n-1)
ans =

   7.5000

>> var(x,1) %除以n

ans =
   6.6667

>> w=[1 1 1 1 1 0.5 0.5 0.5 0.5] %設為權值，向量與x同
w =

   1.0000    1.0000    1.0000    1.0000    1.0000    0.5000    0.5000    0.5000    0.5000

>> var(x,w)
ans =

   5.9184

11.2.5 標準差STD

11.2.5標準差STD

標準誤差為各樣本品與平均值間之常態差，其值實際上為上述變方var執行結果之開方值，其計算公式如下：




　其中 x或r為其平均值。

第一式為正常化之型式，由於平均值佔去一個自由度，故其樣本數應減一，其結果稱為非偏差變方值，通常採用var(X)之指令型式。當樣本數增多時，可以使用第二式計算，

其指令格式如下：

 s = std(X)
 s = std(X,flag)
 s = std(X,flag,dim)

上述之指令均會回應一個標準差值s。若X為矩陣，則會產生一個列向量；令dim=2則會產生行向量之標準差值。

如前述變方之情況，正常化之樣本數應減去一。此為指令之預設值，或flag=0。若要採用分母N為母數時，可令flag=1。下面為利用隨機函數產生Ｘ再計算變方、行向及列向標準差：

>> X=rand(3,4)*10,V=var(X),X1=std(X),X2=std(X,1),X3=std(X,0,2)

X =
  8.3812    3.7948    7.0947    1.8965
  0.1964    8.3180    4.2889    1.9343
  6.8128    5.0281    3.0462    6.8222
V =
 18.8712    5.4672    4.3013    8.0259 %變方值
X1 =
  4.3441    2.3382    2.0739    2.8330　％標準差(/N-1)
X2 =
  3.5469    1.9091    1.6934    2.3131　％標準差(/N)
X3 =
  2.9757
  3.5149
  1.7976

11.2.6 共方差COV

11.2.6 共方差COV

共方差為兩向量之觀察值與其平均差之乘積和，其計算之函數定義如下：

 cov(x1,x2)=E[(x1-u1)(x2-u2)]
 　其中　ui=Exi 

根據上式，其相關指令格式為：

C = cov(X)
C = cov(x,y)

在COV之指令，若 X為向量，其回應值應為變方值。若其為矩陣，則各列為觀察值，各行則成為變數，而COV(X)則為共方矩陣。其對角線元素 DIAG(COV(X))即為每行之變方差向量。若將之排序後，即SQRT(DIAG(COV(X)))，其結果為標準差之向量。以下為例：

>> x=[4    -2     1;  9     5     7]
C1=cov(x)
M=diag(cov(x))
sort(diag(cov(x)))  

x =
     4    -2     1
     9     5     7
C1 =
   12.5000   17.5000   15.0000
   17.5000   24.5000   21.0000
   15.0000   21.0000   18.0000
M =
   12.5000
   24.5000
   18.0000
ans =
   12.5000
   18.0000
   24.5000  

COV(X,Y)則是求取 X與 Y兩等長度之向量之共方差， X與 Y兩向量即使為列向量亦會自動改為行向量，其效果等於COV( [X(:) Y(:)] )。這兩個指令均設法加以常態化，故母數除以N-1，以消除偏差。若要維持使用N為母數，則可增加參數１，即採用 COV(X,1) 或 COV(X,Y,1)指令之型式。

11.2.7 相關係數 corrcoef

11.2.7相關係數 corrcoef

兩個變數相關性可由相關係數求得。其指令型式如下：

R = corrcoef(X)
R = corrcoef(x,y)
[R,P]=corrcoef(...)
[R,P,RLO,RUP]=corrcoef(...)
[...]=corrcoef(...,'param1',val1,'param2',val2,...)

基本上，R=CORRCOEF(X)在於計算一個R矩陣，其內有Ｘ陣列行間之相關係數。而CORRCOEF(X,Y)則計算 X 與 Y兩行向量之相關係數，其意義與CORRCOEF([X Y])相同。

假設 C為共方矩陣，且 C = COV(X)，則R=CORRCOEF(X)之定義為：

   R(i,j) = C(i,j)/sqrt{C(i,i)C(j,j)}

輸入項中除ＸＹ等資料矩陣外，尚可輸入其他特定變數與常數。這些可以用 'PARAM1',VAL1成對表示，其項目包括：

     參數值：
      'alpha'    顯著水準，預設值為0.05(即95%信任區間)
      'rows'     使用 'all' (預設值)表示使用所有列值；
                 'complete'表示使用沒有含NaN 值之列；
                 'pairwise'表示計算R(i,j)時使用不含
                 NaN值之 i行或 j行。

輸出值中， P表示檢驗無關係假設之Ｐ值矩陣。每一個Ｐ值代表隨機可以觀察得到之最大值域。若 P(i,j)值很小，例如小於 0.05，則R(i,j) 之關係甚為顯著。

此外，有RLO與RUP代表95%信任水準之下限與上限矩陣，其大小與Ｒ相同。

例一

>> x=1:5
x =

     1     2     3     4     5

>> y=x.^3
y =

     1     8    27    64   125

>> r=corrcoef(x,y)
r =

    1.0000    0.9431
    0.9431    1.0000

答案中r之值愈接近於1，其相關性愈高。此例中，對角線為自己對自己(即x對x；y對y)故其相關性為1，其餘x對y或y對x，兩者相關性一樣，其數值為0.9431，也相當高。

例二

利用常態分配亂數指令randn產生30X4大小之資料，開始時先利用第四行建立與其他行間之關係，以橫向加總於第四行。其後以corrcoef求相關係數r及機率p。就機率而言，p值愈小，表示兩者之相異性更強，其結果可利用find指令找出小於0.05以下之機率項目。

       x = randn(20,4);       % uncorrelated data
       x(:,4) = sum(x,2)     % introduce correlation
       [r,p] = corrcoef(x)   
               % compute sample correlation and p-values
       [i,j] = find(p<0.05); 
               % find significant correlations [i,j] 

x =
   0.0828   -0.5703   -0.0716   -0.5850
   0.7662   -1.4986   -2.4146   -4.2576
   2.2368   -0.0503   -0.6943    2.2430
   0.3269    0.5530   -1.3914   -0.0113
   0.8633    0.0835    0.3296    0.7592
   0.6794    1.5775    0.5985    2.2962
   0.5548   -0.3308    0.1472   -0.3822
   1.0016    0.7952   -0.1014    2.6212
   1.2594   -0.7848   -2.6350   -2.4089
   0.0442   -1.2631    0.0281   -1.3408
  -0.3141    0.6667   -0.8763   -1.7822
   0.2267   -1.3926   -0.2655   -1.1188
   0.9967   -1.3006   -0.3276    2.0588
   1.2159   -0.6050   -1.1582   -0.2577
  -0.5427   -1.4886    0.5801   -2.8740
   0.9122    0.5585    0.2398    1.9573
  -0.1721   -0.2774   -0.3509   -2.2362
  -0.3360   -1.2937    0.8921   -0.5890
   0.5415   -0.8884    1.5783   -0.4617
   0.9321   -0.9865   -1.1082   -0.4434
r =
   1.0000    0.1950   -0.3475    0.5143
   0.1950    1.0000    0.0929    0.5785
  -0.3475    0.0929    1.0000    0.3822
   0.5143    0.5785    0.3822    1.0000
p =
   1.0000    0.4100    0.1333    0.0203
   0.4100    1.0000    0.6969    0.0075
   0.1333    0.6969    1.0000    0.0963
   0.0203    0.0075    0.0963    1.0000
ans =
    4     1
    4     2
    1     4
    2     4 

11.2.8 群組函數grpstats

11.2.8群組函數grpstats

前面討論到之平均值求法，通常應用於整個陣列之值，若要應用到比較複雜的分組平均問題，則必須使用不同的函數才能達成。此項指令之格式如下：

means = grpstats(X, group)
[means, sem, counts, name] = grpstats(X, group, whichstats)
grpstats(x, group, alpha)

輸入參數中Ｘ為求平均值之對象，Ｘ可為多行，其平均結果也會多行。group則為與Ｘ同列長之陣列，可能由多項分組之向量組成，其內容可為字串列或細胞陣列之文字，如{G1 G2 G3}。若Ｘ中之元素同屬分組中之一項，則其平均值會出現在該項下。

在輸出項中，第一項means為群組平均，sem為組內標準差，counts為各組之項數，name則為各組之名稱。上述項目並非一成不變，亦可以在輸入參數whichstats內依自己之需要進行設定，這個設定有特定的名稱，其名稱必須使用細胞陣列。項目包括：

• 'mean' 組平均
• 'sem'　組標準差
• 'numel'　各組之數目
• 'gname'　各組之名稱
• 'std'　標準差
• 'var'　變異值
• 'meanci'　平均值之95%上下範圍
• 'predci'　新值預測之95%信任範圍

輸入參數中有alpha，可改變其顯著水準，其預設值為0.05，或為95%之信任水準。

輸出項中，means 即為各分組項之平均值，sem為該分組項之標準差，counts為該分組下之觀察值數目，而name則為該分組之名稱。

範例一：

　   x = 1     2     3     4     5     6     7     8     9    10
　group= 1     1     1     1     1     2     2     2     2     2;

>> grpstats(x,group)

ans =

3
8

上述結果為分兩組的平均，前五項為一組，後五項為一組。結果第一組平均為3，第二組為8。
組別間，其項數並不一定要相同，例如：

範例二：

x =
  1     2     3     4     5     6     7     8     9

>> group=[1 1 1 1 2 2 2 2 2]
group =

  1     1     1     1     2     2     2     2     2

>> [m,s,c]=grpstats(x,group)

m =

 2.5000
 7.0000


s =

 0.6455
 0.7071

c =

  4
  5

其輸出之第一項為平均值，第二項為標準差，第三項為各組之項數。故即使各組之樣本數不同也可以得到對應組之統計資料。

範例三：

設有200個觀測值分成四小組，每一觀測值分成五項，其平均範圍由1-5。為製造這樣的數據，下面之例子實際上應用了許多特定的函數：

• unidrnd(4,100,1)　平均製造一個100X1的陣列，其中之數值分配為1:4的整數範圍，以每項分別以1,2,3,4隨機出現。
• normrnd(true_mean,1)　常態分配之亂數函數，其平均值為true_mean，其標準差為1。
• true_mean((ones(100,1),:)　利用原來設定之(ones(100,1),:)陣列，重覆100次。

執行此程式後，由於n為細胞陣列，故全改為字串才能同時顯現其結果，其結果如下：

group = unidrnd(4,100,1);%create a 100X1 matrix in random [ 1,2,3,4 ]
true_mean = 1:5;
true_mean = true_mean(ones(100,1),:); %100X5 matrix
x = normrnd(true_mean,1); %randomize
[m, s, c,n] = grpstats(x,group);

[n num2cell(m)]
ans =
'1'    [0.9584]    [1.8200]    [2.8412]    [4.1669]    [5.0220]
'2'    [0.8972]    [1.8393]    [2.9423]    [4.0044]    [4.9437]
'3'    [0.9768]    [2.1093]    [3.1565]    [3.9860]    [5.0585]
'4'    [1.1164]    [2.2249]    [2.8920]    [4.1323]    [5.3251]

範例四：

利用matlab所附的carsmall.mat示範檔案，其中參數項目包括重量(Weight)、年份(Model_Year)等資料，利用該項資料求其年份下之平均車重、預測值、年份名稱及各年份下之數量。最後並利用errorbar繪出其範圍。

% cargroup.m
load carsmall
[Weight Model_Year]'
[means,p,year,count] = grpstats(Weight,Model_Year,...
 {'mean','predci','gname','numel'})
n = length(means);
errorbar((1:n)',means,p(:,2)-means)
set(gca,'xtick',1:n,'xticklabel',year)
title('95% prediction intervals for mean weight by year')

先將上述程式存為cargroup.m檔案，執行後應有許多資料，其中僅選擇本題所需要者。其過程如下：

>> cargroup

ans =

Columns 1 through 7

     3504        3693        3436        3433        3449        4341        4354
       70          70          70          70          70          70          70

Columns 8 through 14

     4312        4425        3850        3090        4142        4034        4166
       70          70          70          70          70          70          70

＝　＝＝　＝＝　＝＝　＝＝　＝＝　＝＝　＝＝　＝＝　＝＝　＝＝　＝＝
Columns 92 through 98

     2835        2665        2370        2950        2790        2130        2295
       82          82          82          82          82          82          82

Columns 99 through 100

     2625        2720
       82          82


means =

1.0e+003 *

 3.4413
 3.0787
 2.4535


p =

1.0e+003 *

 1.7770    5.1056
 1.3832    4.7742
 1.7184    3.1887


year =

 '70'
 '76'
 '82'


count =

 35
 34
 31

11.2.9 表列函數tabulate

11.2.9表列函數tabulate

TABLE = tabulate(x)

這個指令可將一向量Ｘ之觀測值製成一表格，其第一行為Ｘ向量中之相同數值，第二行為該數值出現之次數，最後一行為該值出現之百分比。若Ｘ值為含有文字串之陣列或細胞陣列，則第一行為陣列內之獨一的名稱，其餘兩行則相同。下面為利用rand之隨機函數取100個值作比較。

tabulate(round(rand(100,1)*6))  
  Value    Count   Percent
      0        8      8.00%
      1       16     16.00%
      2       22     22.00%
      3       21     21.00%
      4       10     10.00%
      5       16     16.00%
      6        7      7.00%  


tabulate(fix(rand(100,1)*6))  
  Value    Count   Percent
      0       14     14.00%
      1       14     14.00%
      2       21     21.00%
      3       15     15.00%
      4       22     22.00%
      5       14     14.00%  

由上面之結果比較可以看出：利用四捨五入法之round函數，其在兩端之機率顯然少很多。
利用前節之汽車carsmall.mat資料，亦可以tabulate指令作簡單之統計：

load carsmall
>> tabulate(Origin)
    Value    Count   Percent
      USA       69     69.00%
   France        4      4.00%
    Japan       15     15.00%
  Germany        9      9.00%
   Sweden        2      2.00%
    Italy        1      1.00%

看起來美國地區還是美國車多。

上述tablulate指令之左邊若不給予參數，則會自動產生上述之格式，分成三行，即名稱、數量及百分比。若結果給予一個參數時，其內容會轉為細胞陣列。因此必要時，須利用cell2mat函數轉換為數值陣列。以上述資料為例，下面的型式會有不同的結果：

>> a=tabulate(Origin)

a = 

    'USA'        [69]    [69]
    'France'     [ 4]    [ 4]
    'Japan'      [15]    [15]
    'Germany'    [ 9]    [ 9]
    'Sweden'     [ 2]    [ 2]
    'Italy'      [ 1]    [ 1]

顯然其內容均是細胞陣列，若要作成餅圖，則必須稍微變換：

>> pie(cell2mat(a(:,2)),a(:,1))

11.2.10 百分值prctile

11.2.10百分值prctile

在一般大量樣本之情況下，可以利用百分值去確定樣本之合理對應值，由此百分比與對應值之關係可以瞭解資料之外形、位置以及擴散度。其指令格式如下：

Y = prctile(X, p)

此指令計算Ｘ之樣本值中一個大於p%部份之對應值位置，此值並不一定是原有之觀測值，只求其比例位置。輸入參數 p 必須落在[0 100]間，可為常數或向量。若 p = 50% 時，則Ｙ值應對應Ｘ之中間值(median)。Ｘ之資料可為向量或矩陣，而 p 則可能為一向量或其中之常數。下表說明可能之之四種狀況：

X      p               Y= prctile(X, p)
====   ========   ===========================================
向量   數值常數   數值等於X中第p級百分值
矩陣   數值常數   列向量含X中每一行內第p級百分值，其中
                      y(i)=prctile(X(:,j),p)
向量   向量       Ｙ向量與p同，其第i項為Ｘ中第p(i)級百分值
                      y(i)=prctile(X,p(i))
矩陣   向量       Y矩陣中，其第(i,j)項為Ｘ中第j行第p(i)級百分值
                      y(i,j)=prctile(X(:,j),p(i))
====   ========   ===========================================

執行例：

>> x=randn(100,1);
>> Y=prctile(x,[10:10:90])
Y =

-1.2258 -1.0402 -0.6829 -0.4320 -0.2074  0.0287  0.3215  0.6855  1.2649

%Y與p之向量相同，與Ｘ為列或行向量無關。

%X2若為矩陣，則p先與Ｘ之行向量作百段分值。

X2 =
  0     1     2     3     4
  1     2     3     4     5
  2     3     4     5     6

>>Y2=prctile(X2,50)
Y2 =
  1     2     3     4     5

>>X2=[0:4;1:5; 2:6],Y2=prctile(X2,[20 30 50])

X2 =
  0     1     2     3     4
  1     2     3     4     5
  2     3     4     5     6

>>Y2 =
 0.1000    1.1000    2.1000    3.1000    4.1000
 0.4000    1.4000    2.4000    3.4000    4.4000
 1.0000    2.0000    3.0000    4.0000    5.0000

X2若為矩陣，p為向量時，每個p(i)先與Ｘ之行向量作百段分值，並Ｙ依序完成行向量之值。

11.2.11細分值quantile

11.2.11細分值quantile

細分值與百分值之意義類似，但其區間為小數，介於[0 1]之間，以配合累積密度函數之使用，其指令格式如下：

Y = quantile(X, p)
Y = quantile(X, p, dim)

其輸出值Ｙ為Ｘ觀測值中傳回值，p為數值或累積機率值之向量。當Ｘ為向量時，Ｙ之大小與p相同，而Y(i)則是第p(i)之細分對應值。當Ｘ為矩陣，則Ｙ之第i列為第p(i)對Ｘ之行向量之細分值。

其細分方向亦可利用dim設定，但Ｙ在dim指定的方向長度應與p之長度相同。

細分值常應用於累加機率，故其範圍為[0 1]。若Ｘ為一具Ｎ元素之向量，則QUANTILE依下列方式運算：

• 1)Ｘ值經過排序，並細分為(0.5/N), (1.5/n), ..., ((N-0.5)/N) 細分段
• 2) 線性內插法用於計算 (0.5/N) 與 ((N-0.5)/N)間機率之細分段。
• 3) Ｘ中之最大值與最小值作為機率外圍之細段。

X=1:10;y=quantile(X,.50)  
y =
   5.5000 

y = quantile(X,[.025 .25 .50 .75 .975]) 
y =
   1.0000    3.0000    5.5000    8.0000   10.0000 

11.2.12堆疊矩陣REPMAT

堆疊矩陣之使用，前面也曾述及，其相關語法如下：

B = repmat(A,m,n)
B = repmat(A,[m n])
B = repmat(A,[m n p...])

這是一個處理大矩陣且內容有重複時使用之。其功能是以Ａ之內容堆疊在一(M x N)的矩陣B中。B矩陣之大小由MXN及A矩陣之內容決定。例如：

>>B=repmat( [1 2;3 4],2,3)

B =
   1     2     1     2     1     2
   3     4     3     4     3     4
   1     2     1     2     1     2
   3     4     3     4     3     4

其結果變為4X6。Ａ也可以置放文字串，如：

>>C=repmat(' Long live the king!', 2,2)
C =
Long live the king! Long live the king!
Long live the king! Long live the king!

也可置放特定常數：

>> D=repmat(NaN,2,5)

D =
 NaN   NaN   NaN   NaN   NaN
 NaN   NaN   NaN   NaN   NaN

11.2.13中央慣性矩MOMENT

11.2.13中央慣性矩MOMENT

統計分析之技巧中，通常以觀察值之平均值為參考標，經處理後再進行判定。MOMENT的指令就是採用相差距離之次方作為計算之基準，如：

mk = E(x-u)k

式中之m值變成以不同k值產生之中央慣性矩，也是各種差異度之期望值。k值為與均值差之次方，又稱階數(order)，其值必須為正。利用觀察資料依各值對中央均值之值距作階數之乘方而累加之和，再除以樣本數。當k=1時稱為第一慣性矩，所得結果應為零，亦即為均值之位置；k=2時即為第二慣性矩，即為上述討論之變方，只是此時之除數為n。其指令型式如下：

m = moment(X,order)

舉例：

>> X=randn(4,6)
m1=moment(X,1)
m2=moment(X,2)
m3=moment(X,3)
m4=moment(X,4)

X =
-1.7258    0.1387    1.1508   -0.3735   -1.5731   -0.2433
0.8132   -0.8595   -0.6080   -0.8320    2.0157    0.1733
1.4419   -0.7523    0.8062    0.2869   -0.0720    0.9232
0.6723    1.2296    0.2171   -1.8189    2.6289   -0.1786
m1 =
 0     0     0     0     0     0
m2 =
1.4524    0.7053    0.4445    0.5872    2.8011    0.2149
m3 =
-1.6611    0.3293   -0.1237   -0.1293   -1.1069    0.0795
m4 =
4.6600    0.8526    0.3402    0.6391   11.1516    0.0919

上述之結果可知，m1均為零，m2與m4因為偶數次方故均為正數；m3則有正負。

11.2.14歪度SKEWNESS

11.2.14歪度SKEWNESS

歪度為第三中央慣性矩除以標準差之三次方，由此可以測出分佈的偏向，或稱為費雪(Fisher)歪度公式。其計算方式如下：

  s = E(x-μ)3/σ3

其中μ為x之平均值，σ為標準差，E表示括號內之期望值。y即為所謂之歪度。其指令之型式如下：

 s = skewness(X)
 s = skewness(X,flag)

輸入參數flag為校正偏差用，flag=1時不作校正，為其預設值；flag=0時則有校正。通常少量之樣本代表一群體時，所得歪度會產生偏差，故視樣本大小需進行校正。因此，修正時，必須執行skewness(X,0)指令。

若x為向量，則s 會得到樣本歪度，或為單一值。當Ｘ為矩陣時，s歪度將以列向量表示。若其結果s等於零，表示對稱分佈；s>0則為正偏(即高峰偏左)；s<0為負偏(高峰在右)。進行判斷時，若其絕對值小於1.96，則屬常態分配，大於1.96則為非常態。例如：

>> X=randn(4,5)
s=skewness(X)

X =
 2.0941    1.6820   -0.1586   -0.5266    0.2486
 0.0802    0.5936    0.8709   -0.6855    0.1025
-0.9373    0.7902   -0.1948   -0.2684   -0.0410
 0.6357    0.1053    0.0755   -1.1883   -2.2476
s =
 0.2794    0.4979    0.9682   -0.4978   -1.1201

由歪度結果得知，這此數值雖屬常態分佈，但也有中央偏向問題。第一、二、三行之高峰偏左，第四及第五行則高峰偏右。

11.2.15峰度KURTOSIS

11.2.15峰度KURTOSIS

峰度是採用第四中央慣性矩除以標準差的四次方而得。其公式如下：

  k =E(x-μ)4/σ4

此處k值稱為kurtosis，其指令格式如下：

k = kurtosis(X)
k = kurtosis(X,flag)

其中flag之功能與skewness相同，為校正偏差用，flag=1時為其預設值，不作校正；flag=0作校正。kurtosis(X)之意義在於表示峰度之趨勢分佈。在常態分配中，峰度值為3。大於3時其峰度將高於常態峰度；小於3時則低於常態峰度。以此可以作為山峰高度之判斷。
例如：

X=randn(100,5);k=kurtosis(X) 

k =
   3.1023    2.6613    2.7877    3.5796    2.7580 

這種常態分佈狀態下，當樣本增多時，峰度值將趨近於3。上述五行之資料當中，第一及四行之山峰較高，其餘之峰頂則較低。（峰度值3即表示標準的山形的意識，也比較容易記）

11.3 統計指令一覽表

在大部份之機率分佈計算中，有些指令尚提供參數及信任水準的計算方法。為獲得相關指令，可以使用help stats進行查詢。例如：

>> help stats
  Statistics Toolbox
  Version 5.0 (R14) 05-May-2004
 
  
Distributions.
   Parameter estimation.
    betafit     - Beta parameter estimation.
    binofit     - Binomial parameter estimation.
    dfittool    - Distribution fitting tool.
    evfit       - Extreme value parameter estimation.
    expfit      - Exponential parameter estimation.
    gamfit      - Gamma parameter estimation.
    lognfit     - Lognormal parameter estimation.
    mle         - Maximum likelihood estimation (MLE).
    mlecov      - Asymptotic covariance matrix of MLE.
    nbinfit     - Negative binomial parameter estimation.
    normfit     - Normal parameter estimation.
    poissfit    - Poisson parameter estimation.
    raylfit     - Rayleigh parameter estimation.
    unifit      - Uniform parameter estimation.
    wblfit      - Weibull parameter estimation.
 
   
Probability density functions (pdf).
    betapdf     - Beta density.
    binopdf     - Binomial density.
    chi2pdf     - Chi square density.
    evpdf       - Extreme value density.
    exppdf      - Exponential density.
    fpdf        - F density.
    gampdf      - Gamma density.
    geopdf      - Geometric density.
    hygepdf     - Hypergeometric density.
    lognpdf     - Lognormal density.
    mvnpdf      - Multivariate normal density.
    nbinpdf     - Negative binomial density.
    ncfpdf      - Noncentral F density.
    nctpdf      - Noncentral t density.
    ncx2pdf     - Noncentral Chi-square density.
    normpdf     - Normal (Gaussian) density.
    pdf         - Density function for a specified distribution.
    poisspdf    - Poisson density.
    raylpdf     - Rayleigh density.
    tpdf        - T density.
    unidpdf     - Discrete uniform density.
    unifpdf     - Uniform density.
    wblpdf      - Weibull density.
 
   
Cumulative Distribution functions (cdf).
    betacdf     - Beta cdf.
    binocdf     - Binomial cdf.
    cdf         - Specified cumulative distribution function.
    chi2cdf     - Chi square cdf.
    ecdf        - Empirical cdf (Kaplan-Meier estimate).
    evcdf       - Extreme value cumulative distribution function.
    expcdf      - Exponential cdf.
    fcdf        - F cdf.
    gamcdf      - Gamma cdf.
    geocdf      - Geometric cdf.
    hygecdf     - Hypergeometric cdf.
    logncdf     - Lognormal cdf.
    nbincdf     - Negative binomial cdf.
    ncfcdf      - Noncentral F cdf.
    nctcdf      - Noncentral t cdf.
    ncx2cdf     - Noncentral Chi-square cdf.
    normcdf     - Normal (Gaussian) cdf.
    poisscdf    - Poisson cdf.
    raylcdf     - Rayleigh cdf.
    tcdf        - T cdf.
    unidcdf     - Discrete uniform cdf.
    unifcdf     - Uniform cdf.
    wblcdf      - Weibull cdf.
 
   
Critical Values of Distribution functions.
    betainv     - Beta inverse cumulative distribution function.
    binoinv     - Binomial inverse cumulative distribution function.
    chi2inv     - Chi square inverse cumulative distribution function.
    evinv       - Extreme value inverse cumulative distribution function.
    expinv      - Exponential inverse cumulative distribution function.
    finv        - F inverse cumulative distribution function.
    gaminv      - Gamma inverse cumulative distribution function.
    geoinv      - Geometric inverse cumulative distribution function.
    hygeinv     - Hypergeometric inverse cumulative distribution function.
    icdf        - Specified inverse cdf.
    logninv     - Lognormal inverse cumulative distribution function.
    nbininv     - Negative binomial inverse distribution function.
    ncfinv      - Noncentral F inverse cumulative distribution function.
    nctinv      - Noncentral t inverse cumulative distribution function.
    ncx2inv     - Noncentral Chi-square inverse distribution function.
    norminv     - Normal (Gaussian) inverse cumulative distribution function.
    poissinv    - Poisson inverse cumulative distribution function.
    raylinv     - Rayleigh inverse cumulative distribution function.
    tinv        - T inverse cumulative distribution function.
    unidinv     - Discrete uniform inverse cumulative distribution function.
    unifinv     - Uniform inverse cumulative distribution function.
    wblinv      - Weibull inverse cumulative distribution function.
 
   
Random Number Generators.
    betarnd     - Beta random numbers.
    binornd     - Binomial random numbers.
    chi2rnd     - Chi square random numbers.
    evrnd       - Extreme value random numbers.
    exprnd      - Exponential random numbers.
    frnd        - F random numbers.
    gamrnd      - Gamma random numbers.
    geornd      - Geometric random numbers.
    hygernd     - Hypergeometric random numbers.
    iwishrnd    - Inverse Wishart random matrix.
    lognrnd     - Lognormal random numbers.
    mvnrnd      - Multivariate normal random numbers.
    mvtrnd      - Multivariate t random numbers.
    nbinrnd     - Negative binomial random numbers.
    ncfrnd      - Noncentral F random numbers.
    nctrnd      - Noncentral t random numbers.
    ncx2rnd     - Noncentral Chi-square random numbers.
    normrnd     - Normal (Gaussian) random numbers.
    poissrnd    - Poisson random numbers.
    randg       - Gamma random numbers (unit scale).
    random      - Random numbers from specified distribution.
    randsample  - Random sample from finite population.
    raylrnd     - Rayleigh random numbers.
    trnd        - T random numbers.
    unidrnd     - Discrete uniform random numbers.
    unifrnd     - Uniform random numbers.
    wblrnd      - Weibull random numbers.
    wishrnd     - Wishart random matrix.
 
   
Statistics.
    betastat    - Beta mean and variance.
    binostat    - Binomial mean and variance.
    chi2stat    - Chi square mean and variance.
    evstat      - Extreme value mean and variance.
    expstat     - Exponential mean and variance.
    fstat       - F mean and variance.
    gamstat     - Gamma mean and variance.
    geostat     - Geometric mean and variance.
    hygestat    - Hypergeometric mean and variance.
    lognstat    - Lognormal mean and variance.
    nbinstat    - Negative binomial mean and variance.
    ncfstat     - Noncentral F mean and variance.
    nctstat     - Noncentral t mean and variance.
    ncx2stat    - Noncentral Chi-square mean and variance.
    normstat    - Normal (Gaussian) mean and variance.
    poisstat    - Poisson mean and variance.
    raylstat    - Rayleigh mean and variance.
    tstat       - T mean and variance.
    unidstat    - Discrete uniform mean and variance.
    unifstat    - Uniform mean and variance.
    wblstat     - Weibull mean and variance.
 
   Likelihood functions.
    betalike    - Negative beta log-likelihood.
    evlike      - Negative extreme value log-likelihood.
    explike     - Negative exponential log-likelihood.
    gamlike     - Negative gamma log-likelihood.
    lognlike    - Negative lognormal log-likelihood.
    nbinlike    - Negative likelihood for negative binomial distribution.
    normlike    - Negative normal likelihood.
    wbllike     - Negative Weibull log-likelihood.
 
  
Descriptive Statistics.
    bootstrp    - Bootstrap statistics for any function.
    corr        - Linear or rank correlation coefficient.
    corrcoef    - Linear correlation coefficient with confidence intervals.
    cov         - Covariance.
    crosstab    - Cross tabulation.
    geomean     - Geometric mean.
    grpstats    - Summary statistics by group.
    harmmean    - Harmonic mean.
    iqr         - Interquartile range.
    kurtosis    - Kurtosis.
    mad         - Median Absolute Deviation.
    mean        - Sample average (in MATLAB toolbox).
    median      - 50th percentile of a sample.
    moment      - Moments of a sample.
    nanmax      - Maximum ignoring NaNs.
    nanmean     - Mean ignoring NaNs.
    nanmedian   - Median ignoring NaNs.
    nanmin      - Minimum ignoring NaNs.
    nanstd      - Standard deviation ignoring NaNs.
    nansum      - Sum ignoring NaNs.
    nanvar      - Variance ignoring NaNs.
    prctile     - Percentiles.
    quantile    - Quantiles.
    range       - Range.
    skewness    - Skewness.
    std         - Standard deviation (in MATLAB toolbox).
    tabulate    - Frequency table.
    trimmean    - Trimmed mean.
    var         - Variance (in MATLAB toolbox).
 
  
Linear Models.
    addedvarplot - Created added-variable plot for stepwise regression.
    anova1      - One-way analysis of variance.
    anova2      - Two-way analysis of variance.
    anovan      - n-way analysis of variance.
    aoctool     - Interactive tool for analysis of covariance.
    dummyvar    - Dummy-variable coding.
    friedman    - Friedman's test (nonparametric two-way anova).
    glmfit      - Generalized linear model fitting.
    glmval      - Evaluate fitted values for generalized linear model.
    kruskalwallis - Kruskal-Wallis test (nonparametric one-way anova).
    leverage    - Regression diagnostic.
    lscov       - Least-squares estimates with known covariance matrix.
    lsqnonneg   - Non-negative least-squares.
    manova1     - One-way multivariate analysis of variance.
    manovacluster - Draw clusters of group means for manova1.
    multcompare - Multiple comparisons of means and other estimates.
    polyconf    - Polynomial evaluation and confidence interval estimation.
    polyfit     - Least-squares polynomial fitting.
    polyval     - Predicted values for polynomial functions.
    rcoplot     - Residuals case order plot.
    regress     - Multivariate linear regression.
    regstats    - Regression diagnostics.
    ridge       - Ridge regression.
    robustfit   - Robust regression model fitting.
    rstool      - Multidimensional response surface visualization (RSM).
    stepwise    - Interactive tool for stepwise regression.
    stepwisefit - Non-interactive stepwise regression.
    x2fx        - Factor settings matrix (x) to design matrix (fx).
 
  
Nonlinear Models.
    nlinfit     - Nonlinear least-squares data fitting.
    nlintool    - Interactive graphical tool for prediction in nonlinear models.
    nlpredci    - Confidence intervals for prediction.
    nlparci     - Confidence intervals for parameters.
 
  Design of Experiments (DOE).
    bbdesign    - Box-Behnken design.
    candexch    - D-optimal design (row exchange algorithm for candidate set).
    candgen     - Candidates set for D-optimal design generation.
    ccdesign    - Central composite design.
    cordexch    - D-optimal design (coordinate exchange algorithm).
    daugment    - Augment D-optimal design.
    dcovary     - D-optimal design with fixed covariates.
    ff2n        - Two-level full-factorial design.
    fracfact    - Two-level fractional factorial design.
    fullfact    - Mixed-level full-factorial design.
    hadamard    - Hadamard matrices (orthogonal arrays).
    lhsdesign   - Latin hypercube sampling design.
    lhsnorm     - Latin hypercube multivariate normal sample.
    rowexch     - D-optimal design (row exchange algorithm).
 
  
Statistical Process Control (SPC).
    capable     - Capability indices.
    capaplot    - Capability plot.
    ewmaplot    - Exponentially weighted moving average plot.
    histfit     - Histogram with superimposed normal density.
    normspec    - Plot normal density between specification limits.
    schart      - S chart for monitoring variability.
    xbarplot    - Xbar chart for monitoring the mean.
 
  
Multivariate Statistics.
   Cluster Analysis.
    cophenet    - Cophenetic coefficient.
    cluster     - Construct clusters from LINKAGE output.
    clusterdata - Construct clusters from data.
    dendrogram  - Generate dendrogram plot.
    inconsistent - Inconsistent values of a cluster tree.
    kmeans      - k-means clustering.
    linkage     - Hierarchical cluster information.
    pdist       - Pairwise distance between observations.
    silhouette  - Silhouette plot of clustered data.
    squareform  - Square matrix formatted distance.
 
   
Dimension Reduction Techniques.
    factoran    - Factor analysis.
    pcacov      - Principal components from covariance matrix.
    pcares      - Residuals from principal components.
    princomp    - Principal components analysis from raw data.
    rotatefactors - Rotation of FA or PCA loadings.
 
   
Plotting.
    andrewsplot - Andrews plot for multivariate data.
    biplot      - Biplot of variable/factor coefficients and scores.
    glyphplot   - Plot stars or Chernoff faces for multivariate data.
    gplotmatrix - Matrix of scatter plots grouped by a common variable.
    parallelcoords - Parallel coordinates plot for multivariate data.
 
   
Other Multivariate Methods.
    barttest    - Bartlett's test for dimensionality.
    canoncorr   - Cannonical correlation analysis.
    cmdscale    - Classical multidimensional scaling.
    classify    - Linear discriminant analysis.
    mahal       - Mahalanobis distance.
    manova1     - One-way multivariate analysis of variance.
    mdscale     - Metric and non-metric multidimensional scaling.
    procrustes  - Procrustes analysis.
 
  
Decision Tree Techniques.
    treedisp    - Display decision tree.
    treefit     - Fit data using a classification or regression tree.
    treeprune   - Prune decision tree or creating optimal pruning sequence.
    treetest    - Estimate error for decision tree.
    treeval     - Compute fitted values using decision tree.
 
  
Hypothesis Tests.
    ranksum     - Wilcoxon rank sum test (independent samples).
    signrank    - Wilcoxon sign rank test (paired samples).
    signtest    - Sign test (paired samples).
    ztest       - Z test.
    ttest       - One sample t test.
    ttest2      - Two sample t test.
 
  
Distribution Testing.
    jbtest      - Jarque-Bera test of normality
    kstest      - Kolmogorov-Smirnov test for one sample
    kstest2     - Kolmogorov-Smirnov test for two samples
    lillietest  - Lilliefors test of normality
 
  
Nonparametric Functions.
    friedman    - Friedman's test (nonparametric two-way anova).
    kruskalwallis - Kruskal-Wallis test (nonparametric one-way anova).
    ksdensity   - Kernel smoothing density estimation.
    ranksum     - Wilcoxon rank sum test (independent samples).
    signrank    - Wilcoxon sign rank test (paired samples).
    signtest    - Sign test (paired samples).
 
  
Hidden Markov Models.
    hmmdecode   - Calculate HMM posterior state probabilities.
    hmmestimate - Estimate HMM parameters given state information.
    hmmgenerate - Generate random sequence for HMM.
    hmmtrain    - Calculate maximum likelihood estimates for HMM parameters.
    hmmviterbi  - Calculate most probable state path for HMM sequence.
 
  
Statistical Plotting.
    andrewsplot - Andrews plot for multivariate data.
    biplot      - Biplot of variable/factor coefficients and scores.
    boxplot     - Boxplots of a data matrix (one per column).
    cdfplot     - Plot of empirical cumulative distribution function.
    ecdfhist    - Histogram calculated from empirical cdf.
    fsurfht     - Interactive contour plot of a function.
    gline       - Point, drag and click line drawing on figures.
    glyphplot   - Plot stars or Chernoff faces for multivariate data.
    gname       - Interactive point labeling in x-y plots.
    gplotmatrix - Matrix of scatter plots grouped by a common variable.
    gscatter    - Scatter plot of two variables grouped by a third.
    hist        - Histogram (in MATLAB toolbox).
    hist3       - Three-dimensional histogram of bivariate data.
    lsline      - Add least-square fit line to scatter plot.
    normplot    - Normal probability plot.
    parallelcoords - Parallel coordinates plot for multivariate data.
    probplot    - Probability plot.
    qqplot      - Quantile-Quantile plot.
    refcurve    - Reference polynomial curve.
    refline     - Reference line.
    surfht      - Interactive contour plot of a data grid.
    wblplot     - Weibull probability plot.
 
  
Statistics Demos.
    aoctool     - Interactive tool for analysis of covariance.
    disttool    - GUI tool for exploring probability distribution functions.
    polytool    - Interactive graph for prediction of fitted polynomials.
    randtool    - GUI tool for generating random numbers.
    rsmdemo     - Reaction simulation (DOE, RSM, nonlinear curve fitting).
    robustdemo  - Interactive tool to compare robust and least squares fits.
 
  
File Based I/O.
    tblread     - Read in data in tabular format.
    tblwrite    - Write out data in tabular format to file.
    tdfread     - Read in text and numeric data from tab-delimitted file.
    caseread    - Read in case names.
    casewrite   - Write out case names to file.
 
  
Utility Functions.
    combnk      - Enumeration of all combinations of n objects k at a time.
    grp2idx     - Convert grouping variable to indices and array of names.
    hougen      - Prediction function for Hougen model (nonlinear example).
    statget     - Get STATS options parameter value.
    statset     - Set STATS options parameter value.
    tiedrank    - Compute ranks of sample, adjusting for ties.
    zscore      - Normalize matrix columns to mean 0, variance 1.

11.4 隨機亂數指令

事件發生的機率雖然可以估計，但是必須由過去發生的事實去統計，並以此預測未來將發生的事件。這種運作的方式在模擬的作業上有實際的困難，故僅能依靠亂數產生器模擬事件之發生。亂數的特性是發生時間或次數均屬隨機性，每次發生的機率均依循其特定的模式，常用者為均勻分佈的亂數型態，每一事件出現之機率是相等的。利用亂數產生之數值可以模擬機械之故障率、顧客來訪的時機、旅客搭機的安排，甚或穀物乾燥中心作業的動線規劃等。

在MATLAB中常用的亂數產生指令為rand，這是個依均勻機率分配出現之原則產生一個或一組在[0, 1]間之亂數，每次呼叫其值均不一樣。這些數值雖有其範圍，但產生之後可以作適當調整，以符合實際所需。所產生之數值可為單一數字，亦可為特定之矩陣，其大小由其後之參數界定之。若僅輸入一個參數(n)，其所產生之結果為一nxn之方矩陣；若為(m, n)則得 m x n之矩陣。

例如：

>>rand  
ans =
    0.9501  

>>rand(2,4)    
ans =
    0.2311    0.4860    0.7621    0.0185
    0.6068    0.8913    0.4565    0.8214  

除均勻機率分配產生者外，亦可使用randn之指令，其產生之亂數則係依常態分態的型式。例如：

>> randn(5,6)

ans =

    0.8956    0.5689   -0.2340    0.6232    0.2379    0.3899
    0.7310   -0.2556    0.1184    0.7990   -1.0078    0.0880
    0.5779   -0.3775    0.3148    0.9409   -0.7420   -0.6355
    0.0403   -0.2959    1.4435   -0.9921    1.0823   -0.5596
    0.6771   -1.4751   -0.3510    0.2120   -0.1315    0.4437

利用常態分佈之指令randn可以產生正負值，其值也可能大於一。其他分配所需之亂數產生器可以參考如下之各項指令。注意每項指令之後可能需要相關之參數，讀者可以利用help之輔助指令獲得所需之資訊。

亂數產生器(Random Number Generators)

依據實際之用途，可以有不同之亂數指令可以應用，

    均勻分佈亂數rand
    常態分佈亂數randn
　  常態亂數(Normal (Gaussian))-normrnd.
    貝他亂數(Beta)-betarnd.
    二項式亂數(Binomial)- binornd.
    X's亂數(Chi square)-chi2rnd.
    極值亂數(Extreme value)-evrnd.
    指數亂數(Exponential)-exprnd.
    F亂數(F)-frnd.
    伽瑪亂數(Gamma)-gamrnd.
    伽瑪亂數( Gamma (unit scale))-randg.
    幾何亂數(Geometric)-geornd.
    超幾何亂數(Hypergeometric)-hygernd.
    反威夏亂數(Inverse Wishart)-iwishrnd.
    常態對數亂數(Lognormal)-lognrnd.
    多變異亂數(Multivariate normal)-mvnrnd.
    多變T亂數( Multivariate t)-mvtrnd.
    負二項式亂數(Negative binomial)-nbinrnd .
    非中央F亂數(Noncentral F)-ncfrnd .
    非中央亂數(Noncentral)-nctrnd.
    非中央X's亂數(Noncentral Chi-square)-ncx2rnd.
    波義松亂數(Poisson)-poissrnd.
    定群樣本亂數(finite population)-randsample .
    魏利亂數( Rayleigh )-raylrnd .
    T亂數(T)- trnd .
    離散均勻亂數( Discrete uniform)-unidrnd .
    均勻亂數(Uniform)-unifrnd .
    魏伯亂數(Weibull)-wblrnd  .
    威夏亂數矩陣( Wishart)-wishrnd .

上述諸多指令，可以實際統計上之需要實際選取應用。例如：以亂數的理念，產生一群數字作排例時，可用randperm(n)指令，以產生不同之排列組合：

>> randperm(7)
ans =
     1     3     7     6     2     5     4

>> randperm(7)
ans =
     7     6     3     5     4     2     1

>> randperm(7)
ans =
     4     1     7     2     6     3     5

其結果每次顯示不同組合，均為整數。此指令因此可以模擬擲骰子，若改用randperm(6)，則成為每趟擲六次，顯現排列數據。這個指令與rand指令不同，若使用rand則必須先化成整數，然後利用程式進行模擬擲骰之動作，其過程比較複雜。

MATLAB亦提供一示範指令，randtool。此指令可以示範各種型態之分配機率，利用下拉式清單可以選擇所需的分配型式，並直接設定所需之範圍與變異量。其結果以對應之質方圖顯示。

>> randtool

由於不同之分佈，有不同的指令，不容易記清楚。不過，MATLAB提供一個通用的指令，其型式如下：

random('name', a1,a2,m,n) 

這個指令則可在name處置放各種指令名稱，其後可依需要再加上對應之參數，其中a1、a2為該函數所需之參數，而最後之[m,n]則為最後之矩陣大小，變成一個相當有彈性的指令，例如：

random('normal',0,1,3,3)  

ans =
   -0.1867    2.1832    1.0668
    0.7258   -0.1364    0.0593
   -0.5883    0.1139   -0.0956  

其中之'name'名稱可以為 'beta', 'binomial', 'chi-square' 或 'chi2', 'extreme value' 或 'ev', 'exponential', 'f', 'gamma', 'geometric', 'hypergeometric' 或 'hyge', 'lognormal', 'negative binomial' 或 'nbin', 'noncentral f' 或 'ncf', 'noncentral t' 或 'nct', 'noncentral chi-square'　或 'ncx2', 'normal', 'poisson', 'rayleigh', 't', 'discrete uniform' 或 'unid', 'uniform', 'weibull' 或 'wbl'等等。部份字母符合亦可，且大小寫不拘。如此可以免去背記上述指令之苦。

11.5常態分配

常態分配之亂數中，不像均勻分佈的型態，基本上它會集中在某區域，故最常發生的事件會集中在平均值附近。在均勻分佈型態應用上，常需設定上下界，讓其出現限定在特定範圍內；常態分配則無明顯的上下界限，且由於集中存在平均值附近，故與均值間之距離會有正負差。其表示方法如下：

y = μ+σx

其中μ為平均值，而σ為其標準差，亦即常態分配之亂數值。在MATLAB中以randn之指令函數產生標準差值σ。此指令係以標準差為１，而平均值為0作成亂數結果，故若有某一群組分佈平均為10，標準差為5，且樣本為200點，則產生後必須乘以所需之標準差5再加上平均值10，其做法如下表示：

y=10+5*randn(1, 200)
>>hist(10+5*randn(1,200))    

上例利用randn製造以10為均值，5為誤差值之隨機樣本的方法，亦可直接利用下面的指令來完成：

y = normrnd(10,5,1,200)

此指令中，前二項參數分別為指定之均值及標準誤差，後二項則為產生之矩陣大小設定。

11.6機率分配

這是一個交談式的指令，可以直接變換參數並得其結果：

>> disttool

這是一個展示各種分配及機率分佈之指令。沒有任何參數，但可以在圖中更改所要的參數。開始時可選分佈種類，如常態、貝他、二項式、指數、二極值、Ｆ、伽瑪、幾何、常態對數等等。選定其一後立即會顯示。其後選擇所需函數名稱。僅有機率分配(pdf)、累積分配函數(cdf)等二項。上兩項選定後，即可輸入Ｘ值，求得對應之機率，或選定機率求得對應之Ｘ值。亦可利用滑鼠直接在圖上點出位置，顯示其對應之Ｘ值與機率。此外，在其左下方有參數設定值，可以利用滑尺為之。

11.7 變異數分析ANOVA1

單向變異分析(ANOVA1)旨在尋找不同類別下之資料是否具有相同均值，亦即決定不同類別下，各量測值是否具有差異特性。在線性模式下，單向變異分析是最簡單的狀況。任何量測值可以矩陣表示如下：

yij = αjεij

其中，yij中每行j代表不同類別，並具有類別內之平均值 。在工具箱中存有hogg.mat資料，下載後可以作為執行anova1函數之用，其過程如下：

>> load hogg
>> hogg

hogg =
 24    14    11     7    19
 15     7     9     7    24
 21    12     7     4    19
 27    17    13     7    15
 33    14    12    12    10
 23    16    18    18    20

ANOVA指令之型式如下：

P = ANOVA1(X,GROUP,DISPLAYOPT)

此指令會將Ｘ組合之樣本矩陣中之各行視為比較之組別，以此決定組別間之平均值是否會相等。GROUP為分組向量，其長度應與X之行數相同，其內容即為各組之名稱。指令左邊為結果之機率，其值愈小表示相等之機率愈小，亦即顯著性之差異愈大。GROUP向量可以容許使用字串或細胞陣列，但必須以行表示。若不使用組別名稱則可以空格。

DISPLAYOPT為顯示開關，可用 'on' (預設值)或 'off'決定是否顯示結果圖。若要文字輸出ANOVA表，則可在左邊加一參數如ANOVATAB：

[P,ANOVATAB, stats] = ANOVA1(...)

參數ANOVATAB為細胞陣列，stats則為一些統計參數。

例：

>> [p,tbl,stats]=anova1(hogg)
p =
0.00011971
tbl =
'Source'     'SS'        'df'    'MS'        'F'         'Prob>F' 
'Columns'    [   803]    [ 4]    [200.75]    [9.0076]    [0.00011971]
'Error'      [557.17]    [25]    [22.287]          []              []
'Total'      [1360.2]    [29]          []          []              []
stats =
 gnames: [5x1 char]
      n: [6 6 6 6 6]
 source: 'anova1'
  means: [23.833 13.333 11.667 9.1667 17.833]
     df: 25
      s: 4.7209

在本例中，hogg之資料實際上各行代表不同處理之細菌繁殖情形。利用ANOVA變異分析表之結果亦可做 F-test，以證明不同處理是否仍然具有相同的結果。此例所得之p值約為 0.0001，這是一個相當小的值。換句話說，此種結果強烈顯示不同的處理其結果之差異是顯著的。若就或然率考慮，每10,000次實驗當中只有一次有結果相同的機會。對於研究者而言，這是一個很大的鼓舞，因為至少由統計分析結果可證明：不同處理比較下，應有很大的差異性的。當然，上述p值要正確，其基本假設是各項變異應獨立，且屬常態分配，其變異值也是固定。這些差異性由圖也可以看出端倪。



上述之指令，若要加上各行之名稱，則可另以如下之指令為之：

>> [p,tbl,stats]=anova1(hogg,{'A1','A2','A3','A4','A5'}')

其執行結果與前述無異，但在最後之組別圖中會出現分組之組名。

11.7.1多重變數比較

11.7.1Multiple Comparisons

有些時候光比較眾處理之間有無差異仍然不足，有時必須知道那一對有顯著的差異。當然此時也可以藉助一連串的t-test逐對比較，但這種方式仍有缺點。在做t-test時，通常會先計算t的過程中與一特定值比較。此值之選定通常是認為各處理之均值超過此值很小（例如5%）。當均值有顯著差異時，統計的數據超過該值的機率反而變大，容易造成誤判。利用MULTCOMPARE指令可以解決此問題，並可執行平均值之多重比較。其語法如下：

   COMPARISON = MULTCOMPARE(stats)

參數stats為統計資料，可以直接得自 anova1, anova2, anovan, aoctool等函數之輸出值，其輸出參數COMPARISON則是一個矩陣，每一列中具有五行，第 1-2行是被比較之兩樣本編號，第 3-5 則分別為其差異之最低值、估計值及最高值。

>>COMPARISON = MULTCOMPARE(stats) 
COMPARISON =
           1            2       2.4953         10.5       18.505
           1            3       4.1619       12.167       20.171
           1            4       6.6619       14.667       22.671
           1            5      -2.0047            6       14.005
           2            3      -6.3381       1.6667       9.6714
           2            4      -3.8381       4.1667       12.171
           2            5      -12.505         -4.5       3.5047
           3            4      -5.5047          2.5       10.505
           3            5      -14.171      -6.1667       1.8381
           4            5      -16.671      -8.6667     -0.66193

以第二列為例，其值為 [ 1　 3　 4.1619 　12.167　 20.171]，表示第一組之平均減去第三組之平均為12.167，而在95%信任水準下，其上下差距為[4.1619 　20.171]。這個區間不包括零，故可以認定這兩組之平均有顯著差異。實際上由圖比較，可以更清楚分辨。由圖所示第一組與第二、三、四組有非常明顯的差別，若比較第四組與第五組，其差異也相當明顯。其他組之間則沒有太顯著的差異。在實際圖的應用上，若用滑鼠點上某一組，則具有顯著差異的會顯示出來，不顯著者則呈灰色，如此可以清楚分辨。

除上述之stats參數外，實際上後面亦可依需要另外附加，其格式如下：

COMPARISON = MULTCOMPARE(STATS, 'PARAM1',val1, 'PARAM2',val2,...)

其後添加之'PARAM1'及val1等參數可以是如下之設定：

      'alpha' ：顯著水準（預設值0.05）
     'displayopt'：估計圖顯示開關 'on' (預設)或  'off'。
     'ctype' ：使用之臨界值，如： 'tukey-kramer' (預設), 'dunn-sidak', 'bonferroni','scheffe'等，可以輸入兩個以上，以空白分隔。
     'dimension'   預設值1。如 [1 3]為第1及第3之預測值。
     'estimate'    選擇估計法進行比較。如：
         anova1: 比較類別均值
         anova2:  'column' (預設)或 'row' 平均。
         anovan:  比較樣本邊界平均值。
         aoctool:  'slope'（預設）, 'intercept',或 'pmm' 。
         kruskalwallis:  欄位平均值。
         friedman:  欄位均值。

輸出參數：

   [COMPARISON,MEANS,H] = MULTCOMPARE(...)

左邊除前述之comparison參數外，可加入 MEANS，等於欄位估計量及標準差，屬矩陣型式，其後 H則為圖握把。

11.7.2雙變異分析

單變異與雙變異分析上不同的地方是後者具有兩組類別，其所定義的特性不相同。例如某汽車公司有兩家工廠，每家工廠均生產相同的三車種。此時生產之汽車之里程數作合理的比較時，其差異可能有兩項：其一是甲乙兩工廠間之差異，其二是車型間之差異。 此時工廠別與車型別兩者均會影響汽車之里程數。其間差異可能來自工廠間之製造，也可能是因車種設計或規格上的問題，後者之問題可能與工廠無關。 此外，某工廠也可能對生產某一車種有獨到的製造能力（可能因為有較佳的生產線吧！），卻製造其他車種則與另一工廠不相上下，這種效應稱為加減性，或兩項類別間之交互作用所產生之影響。

雙變異分析(ANOVA)是線性模式之特殊狀況。就汽車之里程數可以表示如下：

　　yijk=μ+α.j+βi.+γij+εijk

其中，y_ijk為里程數之觀察資料（列指標為i，行指標為j，重複標為k）。μ整個矩陣之平均值。α_.j為由於車型j與μ差異之均值，β_i.則是工廠別i與μ差異之均值，γ_ij為交互作用項，其在行向之和或列向之和為零。ε_ijk 為整個隨機之變異量。

由上述汽車之製造例可知，其觀察之里程數變成一種矩陣的型式，具有j行與i列的組別項，此時由行與列構成之組別或稱為處理(Treatment)，對應於行列方向之交叉點即為細胞(Cell)，每個細胞位置必須重複置放樣本觀察數，或稱為重複數( repetition)。若以矩陣表示，此重複數必須置於k方向。

ANOVA2為分析雙變異數之指令。其格式如下：

　　ANOVA2(X, REPS, DISPLAYOPT)

輸入參數Ｘ為觀察之資料矩陣，其行與列均需二項以上，以作為比較之基準。在各行的資料代表一類別；各列者則為另一類別。若每一對應細胞有多個觀察資料，則由REPS指定。若REPS=3，代表每個細胞位置有三筆資料，這些資料依列之類別依序置於列中。因此若REPS=3，則在列中每三筆資料屬其中一組，下三筆屬第二組，如此類推，其總列數因而應為三(REPS)的倍數。另外參數DISPLAYOPT之定義與前述指令之用法相同。

輸出參數則表示如下：

[P, TABLE, stats] = ANOVA2(...)

其中 P為p值之向量，而 TABLE則為細胞陣列，包括ANOVA表的內容。stats則為分析後之統計資料，可以作為其他指令如MULTCOMPARE函數之輸入。

若有非平衡變異資料則可採用ANOVAN指令進行分析，詳細情形可以參考手冊或輔助資料。
茲以工具箱中存有mileage.mat之車型生產資料為例，可以作為本節之分析示範：

>>load mileage
>>mileage
mileage =
        33.3         34.5         37.4
        33.4         34.8         36.8
        32.9         33.8         37.6
       32.6         33.4         36.6
       32.5         33.7         37
　        33         33.9         36.7 

此為三種車型（行）及二家工廠製造之汽車每加侖里程資料，其reps=3，故第1-3列代表第一家工廠製造之車輛，第4-6列代表第二家工廠。將相關參數代入anova2指令，可得結果為：

>>[p,tbl,stats] = anova2(mileage,3)
p =
2.4278e-010    0.0039095      0.84106

tbl =
'Source'         'SS'        'df'    'MS'         'F'          'Prob>F'
'Columns'     [53.351]    [ 2]    [ 26.676]    [ 234.22]    [2.4278e-010]
'Rows'        [ 1.445]    [ 1]    [  1.445]    [ 12.688]    [  0.0039095]
'Interaction' [  0.04]    [ 2]    [   0.02]    [0.17561]    [    0.84106]
'Error'       [1.3667]    [12]    [0.11389]           []               []
'Total'       [56.203]    [17]           []           []               []
stats =
  source: 'anova2'
 sigmasq: 0.11389
colmeans: [32.95 34.017 37.017]
    coln: 6
rowmeans: [34.944 34.378]
    rown: 9
   inter: 1
    pval: 0.84106
      df: 12

除文字及統計數據外，尚有一個標準的 ANOVA表，其中第一行為行，其平方和(SS)、自由度(df)及均方值 (MS=SS/df), F 檢定及或然率 p-值等。 F檢定可以檢驗車型間、工廠間及車型X工廠之交互作用（經過調整後之增加效應）產生之里程數是否相同。由於 p-值在車型間之效應幾乎近於零(至小數四位)，表示車型對里程數之變異甚大，有顯著之影響。經由 F檢定結果，因車型而發生同平均值之情況，其機率為萬分之一。若使用MULTCOMPARE指令比較：

>>COMPARISON = MULTCOMPARE(stats)

Note:  Your model includes an interaction term.  A test of main
effects can be difficult to interpret when the model includes
interactions.
COMPARISON =
        1            2      -1.5865      -1.0667     -0.54686
        1            3      -4.5865      -4.0667      -3.5469
        2            3      -3.5198           -3      -2.4802

利用multcompare函數比較結果，三車型之差異有明顯之不同。而 p-值對工廠間之機率值為 0.0039，這也是一個非常顯著之差異。顯然，某一工廠製造的汽車，其里程數是比另一廠為高。由其 p-值得知，只有千分之四的機率兩工廠製造的汽車之里程數才會相同。就工廠Ｘ車型間交互作用之影響而言，則不顯著。其 p-值僅 0.8411，亦即結果中，可能百分八四機率會出現無交互作用之影響。當然，此處顯示之 p-值若要正確，基本上整個樣本之分佈必須獨立、常態分配並有固定的變異常數。

11.7.3多變異分析ANOVAN

11.7.3多變異分析ANOVAN

在變異分析的過程中，有些資料並非刻意形成，而是因調查或分析產生的特定資訊，無法事先規劃作出等量分組或先經實驗分組。例如，有一些汽車銷售資訊，在銷售單中可能有不同車種、型號、排氣量，甚至出產國別等。若就其分類有些數量並不一定對等，或處於非平衡狀態。多變異分析可以同時處理平衡或非平衡資料，其功能與anova1、anova2略有不同。此模式之架構如下:


其中有相乘之兩項或三項者，即為其間之交互作用(interaction)產生之效應。在所需之Ｘ資料向量中，應有其對應之類別資訊。依其類別資訊可以區分為若干層次，此時可將要設為層次之行歸納在一個細胞陣列 group中。其內可有Ｎ層之名單，以此指認Ｘ資料之歸屬。在group細胞中之名單可以是一向量、文字陣列或文字細胞陣列，以大括符集合起來，惟其個數需與Ｘ中之項目相符。其相關語法如下：

p = anovan(x,group)
p = anovan(x,group,'Param1',val1,'Param2',val2,...)
[p,table] = anovan(...)
[p,table,stats] = anovan(...)
[p,table,stats,terms] = anovan(...)

其參數定義大體上與anova1、anova2等指令相同，每個 GROUP中之變數則必須具有與Ｘ等量之元素。在分析後之變異數分析表中，將列出GROUP中之變數名稱。另外增加之參數'Paramx',valx，則可說明如下：

'Paramx'	valx
'sstype'	1, 2, 或 3（預設值）平方和之型式
'varnames'	組別名稱預設為 'X1', 'X2', 'X3', ..., 'XN' ，但可利用此參數定義新名稱
'display'	顯示ANOVA表，'on' | 'off'
'random'	以一向量指定群組參數是否為隨機。
'alpha'	信任水準（預設值為0.05或95%）
'model'	決定採用何種函數，如 'linear'其p值僅計算Ｎ主效應（預設值）； 'interaction'除主效應外尚計算兩者間之交互效應； 'full'計算所有主效應及交互效應。若為整數k(k<=N)，則表示計算至k層級。如k=3，表示主效應＋兩者間及三者間之交互效應。 若為零壹值矩陣之定義，如[1 0 0] A主效應、[0 1 0] B主效應、[0 0 1 ] C主效應、[1 1 0] AB交互效應、[1 0 1] AC交互效應、[0 1 1] BC交互效應、[1 1 1] ABC交互效應等。

在輸出項當中，TERMS矩陣可作為下一次執行ANOVA時，其 MODEL輸入之引數。若MODEL參數不採用隨機型態，其變異分析表T 中之欄位將為TERMS 、ＳＳ、ＤＦ、奇偶數、均方、Ｆ檢測及Ｐ值等。若具隨機功能時，則會顯示TERMS之型態（即固定或隨機）、期望均方、Ｆ值之除數均方、Ｆ值之除數自由度、除數定義、變異項估計值、低限值及高限值等。

輸出參數STATS 的結構包括使用MULTCOMPARE函數之參數及如下項：

coeffs      估計係數
coeffnames  各係數之名稱
vars        組別參數值矩陣
resid      適配模式下之餘數。

具隨機效應時，則另有下列欄位：

ems         期望均方值。
denom       分母定義。
rtnames     隨機項名稱。
varest      變異項估計值（每一隨機項一個）。
varci       各變異項之信任區間

實際上anova2與anovan也可以通用，只是其安排的語法略有不同而已。anova2是將重複的觀察值置於列向，用reps來計算其每組之重複次數，而anovan則必須另外建相關的行作為指示其階數之用。例如下面的例子，假設屬anova2之資料，且是二重複，即reps=2：

X = [25 17 21; 28 18 65; 45 5 60; 42 10 95]

X =
25    17    21
28    18    65
45     5    60
42    10    95
p=anova2(X,2)
p =
0.0175    0.1930    0.2321

上述資料格式若要改為適用anovan，則必須建立對應之行向量。首先建立一個對應三欄的位置m1；其次再建立對應重複的群組m2：

m1=repmat(1:3,4,1)
m2=[ones(2,3);2*ones(2,3)]

m1 =
1     2     3
1     2     3
1     2     3
1     2     3
m2 =
1     1     1
1     1     1
2     2     2
2     2     2

其次，再將Ｘ，m1，m2以行向對應起來，如：

X=X(:);m1=m1(:);m2=m2(:);
[X m1 m2]
ans =
25     1     1
28     1     1
45     1     2
42     1     2
17     2     1
18     2     1
5     2     2
10     2     2
21     3     1
65     3     1
60     3     2
95     3     2

此時就可以利用anovan執行分析了。

anovan(X, {m1 m2},2)
ans =
0.0175
0.1930
0.2321

其結果與上述anova1相同。看起來使用naovan好像比較複雜，那是因為轉換為anova1之型式複雜的綠故，若準備資料時就朝naovan之需要處理，則會簡單許多。指令中之參數2即為'model'=2的意思，表示分析結果包括交互作用之影響。若不考慮交互作用，則可以免去，其結果將完全與anova1相同。

下面是一個具三方變異分析之例子，設ｙ為觀察之資料，其餘三組g1、g2、g3為對應之組合項，各組內僅有兩種內容，g1為[1 2]、g2['hi' 'lo']、g3['june' 'may']。由三組內容可知，其項目可為數值，可為文字，但文字必須使用文字字串。

y = [52.7 57.5 45.9 44.5 53.0 57.0 45.9 44.0]';
g1 = [1 2 1 2 1 2 1 2];
g2 = {'hi';'hi';'lo';'lo';'hi';'hi';'lo';'lo'};
g3 = {'may'; 'may'; 'may'; 'may'; 'june'; 'june'; 'june'; 'june'};

利用上述資料進行三方變異分析如下：

>>p = anovan(y, {g1 g2 g3})
p =
0.4174
0.0028
0.9140

輸出項Ｐ值包括三個主項所產生之效應，由其變異分析表可知分為x1、x2、x3對應g1、g2、g3等三組。由於沒有第三項參數，故僅列示主效率部份。若需其交互效應則必須加入2之參數在最後一項，或令'model' = 2。三個Ｐ值則對應個別零假設(Null Hypothesis)Ho1、Ho2、Ho3之機率值。基本上Ｐ值趨近零，對於零假設的成分愈小。例如，Ho2之對應Ｐ值0.0028已足夠小到可以認定該組合項下之某一樣本平均值無法與另一組有顯著之差別。在進行此項研斷之前，通常必須先選定Ｐ值之顯著水準門檻值，常用的的有0.005或 0.001，依問題的特性而定。

上述的結果並未考慮交互效應，若要增加此項功能，必須在參數方面設定，即令'mode'=2，或直接使用2即可。例如：

>>[p,table]=anovan(y, {g1 g2 g3},'model',2)

p =
0.0347
0.0048
0.2578
0.0158
0.1444
0.5000
table =
'Source'  'Sum Sq.' 'd.f.' 'Singular?' 'Mean Sq.'    'F'       'Prob>F'
'X1'     [  3.7812] [   1] [   0]    [  3.7812]   336.1111]    [0.0347]
'X2'     [199.0013] [   1] [   0]    [199.0013]  [1.7689e+004] [0.0048]
'X3'     [  0.0612] [   1] [   0]    [  0.0612]  [     5.4444] [0.2578]
'X1*X2'  [ 18.3012] [   1] [   0]    [ 18.3012]  [1.6268e+003] [0.0158]
'X1*X3'  [  0.2112] [   1] [   0]    [  0.2112]  [    18.7778] [0.1444]
'X2*X3'  [  0.0113] [   1] [   0]    [  0.0113]  [          1] [0.5000]
'Error'  [  0.0113] [   1] [   0]    [  0.0113]          []    []
'Total'  [221.3788] [   7] [   0]            []          []    []

結果會有三項值，其內容大致相同。Ｐ值之首三項為主效應，其餘三項為交互效應。在名稱上均預設為x1、x2、x3，此名稱亦可自行設定：

>>p= anovan(y,{g1 g2 g3},'varnames',{'time' 'limit' 'month'},'model',2);

>> aoctool
Note: 6 observations with missing values have been removed.

AOCTOOL 適合單向變異分析(ANOCVA)及繪製變量圖。其執行格式如下：

AOCTOOL(X,Y,G,ALPHA)

上式之輸入參數：X,Y為對應之資料，其中Ｘ為已知值，Ｙ為回應值。Ｇ則是作為分組之變數。而ALPHA則是信任水準，其預設值為0.05，亦即其信任水準為95%。

執行這個程式後，將出現三個圖表，其中一個為交談曲線圖，可以調整參數使其產生預測值變化，其二為為變異分析表(ANOVA)，另一為參數值預測。

除上述輸入參數外， 各參數尚可指定名稱，其圖表亦可開關。其指令型式如下：

AOCTOOL(X,Y,G,ALPHA,XNAME,YNAME,GNAME,DISPLAYOPT,MODEL)

其中XNAME,YNAME,GNAME分別為Ｘ軸、Ｙ軸及組別名稱。 DISPLAYOPT 則控制圖形是否顯示（'on'為預設值，或否'off'）， MODEL 則控制最初適配時之選項，包括：

 1    'same mean'　（共同均值）
 2    'separate means'　（分開均值）
 3    'same line'　（共同預測線）
 4    'parallel lines'　（平行預測線）
 5    'separate lines'　（分離預測線）

此外，此一指令亦容許有輸出，其型式如下：

[H,ATAB,CTAB,STATS] = AOCTOOL(...)

其中 H代表圖中各線之握把， ATAB 與 CTAB則為ANOVA表及係數表之內容，而 STATS則為多重均值比較之結果，亦即為 MULTCOMPARE 函數執行的結果。

下面的例子中為利用load carsmall下載最近十餘年（1970-1982）的汽車相關資料。其中包括車重(Weight)、加侖里程數(MPG)及年份(Model_Year)。其資料內容如下，請特別注意在加侖里程之資料中有六點之資料從缺。

>> Weight'
 3504     3693     3436     3433     3449     4341     4354
 4312     4425     3850     3090     4142     4034     4166
 3850     3563     3609     3353     3761     3086     2372
 2833     2774     2587     2130     1835     2672     2430
 2375     2234     2648     4615     4376     4382     4732
 2464     2220     2572     2255     2202     4215     4190
 3962     4215     3233     3353     3012     3085     2035
 2164     1937     1795     3651     3574     3645     3193
 1825     1990     2155     2565     3150     3940     3270
 2930     3820     4380     4055     3870     3755     2605
 2640     2395     2575     2525     2735     2865     3035
 1980     2025     1970     2125     2125     2160     2205
 2245     1965     1965     1995     2945     3015     2585
 2835     2665     2370     2950     2790     2130     2295
 2625     2720


>> MPG'
18.0000   15.0000   18.0000   16.0000   17.0000   15.0000   14.0000   14.0000
14.0000   15.0000       NaN       NaN       NaN       NaN       NaN   15.0000
14.0000       NaN   15.0000   14.0000   24.0000   22.0000   18.0000   21.0000
27.0000   26.0000   25.0000   24.0000   25.0000   26.0000   21.0000   10.0000
10.0000   11.0000    9.0000   28.0000   25.0000   25.0000   26.0000   27.0000
17.5000   16.0000   15.5000   14.5000   22.0000   22.0000   24.0000   22.5000
29.0000   24.5000   29.0000   33.0000   20.0000   18.0000   18.5000   17.5000
29.5000   32.0000   28.0000   26.5000   20.0000   13.0000   19.0000   19.0000
16.5000   16.5000   13.0000   13.0000   13.0000   28.0000   27.0000   34.0000
31.0000   29.0000   27.0000   24.0000   23.0000   36.0000   37.0000   31.0000
38.0000   36.0000   36.0000   36.0000   34.0000   38.0000   32.0000   38.0000
25.0000   38.0000   26.0000   22.0000   32.0000   36.0000   27.0000   27.0000
44.0000   32.0000   28.0000   31.0000


>> Model_Year'
70    70    70    70    70    70    70    70    70    70    70    70    70    70
70    70    70    70    70    70    70    70    70    70    70    70    70    70
70    70    70    70    70    70    70    76    76    76    76    76    76    76
76    76    76    76    76    76    76    76    76    76    76    76    76    76
76    76    76    76    76    76    76    76    76    76    76    76    76    82
82    82    82    82    82    82    82    82    82    82    82    82    82    82
82    82    82    82    82    82    82    82    82    82    82    82    82    82
82    82

利用AOCTOOL分析時，設車重為輸入值，加侖里程為對應值，年份作為組別，其指令格式如下：

aoctool(Weight,MPG,Model_Year,0.05)

actool

預測圖

　

上圖為依年份所得之預測值，共有三個不同年份，可以依年份之類別選擇顯示。由此圖可知車重增加時，每加侖之里程數將減少。且年份愈近者，其每加侖之里程數愈高。

在模式(Model)中，則可選不同的預測線處理方式。



此表則表示三條預測線之截距與斜率。就截距而言，總截距為45.9798，70年份則比此值略低-8.5805，亦即45.9798¬-8.5805=37.3993。總的斜率為-0.0078，70年份則比此值略高0.002，或為-0.0058。如此類推，即可看出統計上之預測值。其對應之標準差及Ｔ試驗值則附於其後。各組之公式如下：








就變異分析表(ANOVA)分析，年份及重量均相當顯著，其交叉效果也在0.0072之水準，亦在顯著之範圍，此可由圖中82年份與其他年份在車重4000公斤有交叉的現象有關。
若使用參數輸出的方式，則可執行如下指令：

[h,a,c,s] = aoctool(Weight,MPG,Model_Year,0.05)

Note: 6 observations with missing values have been removed.
h =
152.0028
153.0013
154.0013
204.0013
205.0013
206.0013
207.0013
208.0013
a =
Columns 1 through 5
'Source'           'd.f.'    'Sum Sq'         'Mean Sq'        'F'
'Model_Year'      [   2]    [   807.6896]    [   403.8448]   [ 51.9762]
'Weight'          [   1]    [2.0502e+003]    [2.0502e+003]    [263.8679]
'Model_Year*Weight'  [   2]  [    81.2188]   [    40.6094]    [  5.2266]
'Error'            [  88]    [   683.7420]   [     7.7698]            []
Column 6
'Prob>F'
[1.2212e-015]
[          0]
[     0.0072]
        []

c =
'Term'         'Estimate'    'Std. Err.'      'T'           'Prob>|T|'
'Intercept'    [ 45.9798]    [     1.5208]    [ 30.2330]    [2.9266e-048]
'  70'         [ -8.5805]    [     1.9619]    [ -4.3737]    [3.3417e-005]
'  76'         [ -3.8902]    [     1.8686]    [ -2.0818]    [     0.0403]
'  82'         [ 12.4707]    [     2.5568]    [  4.8775]    [4.7391e-006]
'Slope'        [ -0.0078]    [5.5717e-004]    [-14.0031]    [4.0955e-024]
'  70'         [  0.0020]    [6.6152e-004]    [  2.9605]    [     0.0039]
'  76'         [  0.0011]    [6.5328e-004]    [  1.7425]    [     0.0849]
'  82'         [ -0.0031]    [9.9914e-004]    [ -3.0994]    [     0.0026]
s =
 source: 'aoctool'
 gnames: {3x1 cell}
      n: [3x1 double]
     df: 88
      s: 2.7874
  model: 5
 slopes: [3x1 double]
slopecov: [3x3 double]
intercepts: [3x1 double]
intercov: [3x3 double]
    pmm: [3x1 double]
 pmmcov: [3x3 double]

11.8內插法

內插技巧是資料分析時常用的工具，尤其在尋求等距之對應資料作為比對時，常需利用內插法尋求實際對應之值。當然內插法亦可用來尋求資料間之中間對應值。在MATLAB中有下列指令處理內插法：

指令名稱	說明
interp1,interp1q	一維內插法，有如查表法。interp1q與interp1功能相同，但執行速度較快。
interp2	二維內插指令
interp3, interpn	三維及多維內插指令
spline	曲尺內插

以東方人體型之身高與體重標準為例，其對應資料如下：

身高  157.5  160 162.5 165  167.5  170   172.5  175   177.5 180
小體型 52.5  54   55    57   58    59.5  61.5   63.5   66   67.5
中體型 56    57.5 59    60   61.5  63.5  66     67.5   69   71
大體型 60    61.5 63.5  65   67    68.5  70.5   72.5   74   76.5

由上表可直接查得160cm對應之中體型重為57.5kg。但若要知道中間值，例如，161cm身高的人的體重則必須用內插法。所以內插法之使用有一個特徵是其相關值是連續增加或減少或具有相關性的。因此利用內插法較不易失敗。尋求內插值可以使用上述幾個指令，但若屬一維的對應，則仍以使用interp1為多，其呼叫之型式如下：

YI = INTERP1(X,Y,XI)

這個指令類似查表法，前面X,Y其實是一個對應表，後面之XI則是要查之進入值，而所得之結果則置於左邊之YI。此處，X,Y必須是一個對應的向量，XI則可為向量或常數值，其個數亦不必與X,Y相同。只是所得之YI的個數亦會與XI相同。以上述之身高體重為例，若要知道161、163、171cm等高度之對應中體型體重時，則可呼叫如下：

height=[157.5 160 162.5 165 167.5 170 172.5 175 177.5 180];
M_wt=[56 57.5 59 60 61.5 63.5 66 67.5 69 71];
YI=interp1(height,M_wt,[161 163 171])

YI =
58.1000   59.2000   64.5000

利用上項功能，實際上也可以重新定域新範圍，並獲得新的對應值。例如希望能獲得150,155,160,165,179,175,180之對應體重，並由此作圖，則可以執行如下：

H=[160:5:180];
W=interp1(height,M_wt, H);
plot(height,M_wt,H,W,'o');xlabel('Height, cm');ylabel('Weight, kg')

使用同樣的資料，若想查尋58、62、70kg人的對應身高時，亦可利用interp1指令查詢之，只是其順序必須稍加更動：

XI=interp1(M_wt,height,[58 62 70])
XI =
160.8333  168.1250  178.7500

上述之指令中，實際上它內部採用之內插技巧是內設為線性內插，亦即兩點間均以直線作為考量。為求精確，有些內插必須採用特定的方法。在interp1的指令中，實際上亦提供此項選擇，其參數位置如下：

 YI = INTERP1(X,Y,XI,'method')

不同的技巧則在'method'處表示。可以選擇如下幾種：

   'linear'   - 線性內插（預設值）。
   'spline'   - 區段立方曲線內插。
   'nearest'  - 以最靠近的點進行內插。
   'pchip'    - 保留外形區段立方內插（Hermite型）。
   'cubic'    - 與 'pchip'同。

一般之內插法預設為線性，故若需繪圖，並取得較為平滑曲線，則可選用 'spline'，或稱為立方曲線。下面為線性與立方曲線所顯示出不同的效果：

     x = 0:10; y = sin(x); xi = 0:.25:10;
yi = interp1(x,y,xi); plot(x,y,'o',xi,yi)

由於原來之分點之區隔較粗，故僅能顯示直線的曲線。若採用 'spline'的功能，則會產生比較平滑的曲線：

yi = interp1(x,y,xi,'spline'); plot(x,y,'o',xi,yi)

指令中之X參數可以省略，但此時指令用X=[1 2 …length(Y)]為其參數。與rnterp1指令類似的有rnterp1q，後者之功能大體與前者相同，只是後者執行較快，適用於不等距之間隔，因為其輸入值並未經嚴格檢查。但此指令使用時，Ｘ之值必須為單一對應值，且須遞增排列；而Ｙ則需為行向量，若為矩陣則其長度應與X相同。

11.9 外插法

如果要查的資料超高範圍會如何？以前項之身高體重為例：

XI=interp1(M_wt,height,80)
XI =
NaN

顯然沒有資料，因為它已超出範圍。此時若要有資料可以採用外插法，即在後面之參數加上'extrap'，並且要表示外插法所使用之方法，例如'linear'是。

XI=interp1(M_wt,height,80,'linear','extrap')
XI =
191.2500

只是一個體重80kg的人，他的標準身高 191.2cm是否合理，那就不是外插法所能管得了的了。

內插法亦可應用於多維資料中。例如：

x = [1:10]'; y = [ x.^2, x.^3, x.^4 ];
xi = [ 1.5, 2.5; 5.5, 8.5], yi = interp1(x,y,xi)

xi =
1.5000    2.5000
5.5000    8.5000
yi(:,:,1) =
2.5000    6.5000
30.5000   72.5000
yi(:,:,2) =
4.5000   17.5000
170.5000  620.5000
yi(:,:,3) =
1.0e+003 *
0.0085    0.0485
0.9605    5.3285

由於輸入值為2x2之向量，所得之yi具有三個函數，故會產生大小為2x2x3之對應值。

INTERP2 指令

處理函數如z=f(x,y)這種二維的內插問題時，可以使用interp2指令。其呼叫型式如下：

ZI = INTERP2(X,Y,Z,XI,YI)

指令中，X,Y,Z為對應值，其中X,Y為自變數，Ｚ為對應函數值。進行內插時，必須設定XI,YI值，其後結果ZI置於等號之左邊。在輸入參數後邊，亦可加入前述之'method'項，以獲得不同的內插效果。其預設項仍為'linear'，即線性內插。下面的例子為先設定XY方向之方格，建立Ｚ函數；然後以較細的內插點繪出對應點，將其Ｚ座標提高15個單位，以與原函數比較。

[X,Y] = meshgrid(-3:.5:3);
Z = peaks(X,Y);
[XI,YI] = meshgrid(-3:.25:3);
ZI = interp2(X,Y,Z,XI,YI,'spline');
mesh(X,Y,Z), hold, mesh(XI,YI,ZI+15)
hold off
axis([-3 3 -3 3 -5 20])

在interp2指令之參數中，若Ｚ前面之Ｘ、Ｙ兩項省略，如：

ZI = INTERP2(Z,XI,YI)

則表示內部預設 X=1:N，Y=1:M，其中 [M,N]=SIZE(Z)。如此可以省掉一些裝飾的工夫。interp2也可以使用外插法，其設定方式與interp1相同。所有ＸＹ之對應必須單一的，其間隔可以不是均勻分割的。

[x,y]=meshgrid(-5:5,-5:.5:5);
z=x.^2+y.^2;
[xi,yi] = meshgrid(-3:.2:3,-3:.2:3);
zi = interp2(x,y,z,xi,yi);
mesh(xi,yi,zi+50);hold on;
mesh(x,y,z);hold off

11.10 spline指令

由於立方曲線之內插可產生較為平滑的曲線，故常為內插法的選項。為此，MATLAB特別另建一指令處理立方曲線之內插。立方曲線內插係以三次方多項式表示，其基本型式如下：

　yi(x)=c1(x-xi)3+c2(x-xi)²+c3(x-xi)+c4

其中i代表每一待內插之區間， ，若總共有n點，則 其中之係數 則必須利用內插之求法求解。其邊界條界如下：

 此多項點必須通過區段之兩端 。
 相鄰兩個多項式在共同端點之斜率（第一導數）必須一致。
 相鄰兩個多項式在共同端點之曲率（第二導數）必須一致。

要求出上項係數 ，可以使用sprine指令之另一種型式：

　　pp=spline(x,y)

這是一個中間結果的型式，因為實際要知道其係數值的不多，大部份都是直接得到內插對應值即可。上述指令之(x,y)為基本資料輸入，pp則為結構變數。其內有form、breaks、coefs等變數。比較重要的是其breaks之Ｘ值基點及對應之係數coefs。依結構變數之呼叫法即為pp.breaks及pp.coefs。以下為執行例：

x=[10  12 14 17];y=[120 140 150 175];
pp=spline(x,y)
pp.breaks
pp.coefs

pp =
    form: 'pp'
  breaks: [10 12 14 17]
   coefs: [3x4 double]
  pieces: 3
   order: 4
     dim: 1
ans =
  10    12    14    17
ans =
  0.2738   -2.8929   14.6905  120.0000
  0.2738   -1.2500    6.4048  140.0000
  0.2738    0.3929    4.6905  150.0000

因此，由上述得到之係數可以代入內插多項式：



上述之資訊對內插之結果並沒什麼幫助，但有時需要這些係數以作為其他指令之用途。但所得之結果亦可利用y=ppval(pp,x)指令得到內插之結果。在spline(x,y)之指令中，若 Y之元素項比Ｘ多二項時，其最初及最後項將用為設定立方曲線之端點斜率之用。若 Y為一向量，亦即：

      f(X) = Y(2:end-1),  Df(min(X))=Y(1),    Df(max(X))=Y(end)

例如：

     x = -4:4; y = [0 .15 1.12 2.36 2.36 1.46 .49 .06 0];
     cs = spline(x,[0 y 0]);
     xx = linspace(-4,4,101);
     plot(x,y,'o',xx,ppval(cs,xx),'-');

基本上，spline指令亦可直接進行內插，輸入項及可以呼叫型式如下：

 YI = SPLINE(X,Y,XI)

例如：

      x = 0:10;  y = sin(x);
      xx = 0:.25:10;
      yy = spline(x,y,xx);
      plot(x,y,'o',xx,yy)

實際上YI = SPLINE(X,Y, XI)等同 YI = PPVAL(SPLINE(X, Y), XI) 。前面所述PP=SPLINE(X,Y)之指令，也可以使用 PP = UNMKPP(BREAKS, COEFS)或 PP = UNMKPP(SPLINE(X,Y))獲得同樣的效果。其中BREAKS為ｘ之向量，而 COEFS則為三次方多項式之係數。

11.11 多項式適配polytool

多項式繪圖工具polytool

指令polytool函數可以利用交談的方式配合多項式曲線進行繪圖及預測。其語法如下：

polytool(x,y)
polytool(x,y,n)
polytool(x,y,n,alpha)
polytool(x,y,n,alpha,xname,yname)
h = polytool(...)

內容包括一般直線及多項式線 ，其中n代表多項式之階數，其預設值為1；x及y為對應之資料陣列。alpha為信任區間，其預設值為0.05。這個指令配合繪圖，可以將預測之曲線直接繪出。其參數可以自動產生（在無任何輸入參數時），亦可依需要輸入xy之對應觀測值，其他參數之定義與polyfit相同。另外亦可加上Ｘ軸與Ｙ軸之名稱 xname, yname。左邊之h值為其握把。其產生多項式之型式如下：

y=β0+β1x+β2x2+...+βnn+ε

　　

在沒有參數時，此指令自動用下面之函數產生，例如：

x = (1:10)';
y = 50 + 4*x - 0.75*x.^2 + randn(10,1);
>> polytool(x,y,3)

指令polytool之交談圖

在圖中可以改變階數，並可改變X值以求其對應Ｙ值及變化範圍。其相關係數可以輸出至變數空間中作進一步處理。

多項式適配例

在工具箱中，MATLAB有存一筆資料polydata.mat，可以用load指令下載再用此指令執行，進行預測。此指令以適配多項式，預測中間值。

>>load polydata;
>>whos

x         1x43                      344  double array
x1        1x101                     808  double array
y         1x43                      344  double array
y1        1x101                     808  double array

這個檔案中含有四個變數，其中x、y為觀測值，而x1、y1則為真實函數對應值。為使用ploytool進行適配，函數指令如下：

polytool(x,y)

其所得之資料點及適配線如下圖：


由圖中所示，顯然使用預計線性適配之結果並不很好。大部份之資料點均集中在左側，而右側較遠處則有兩點在左側方向。此時若利用其交談功能，將位階改為三次方項，結果為


此時適配之資料顯然較佳，其上下限也較為接近。當然要求更接近之預測值也許可以試著增加其階數。但有時並不一定能到十分美意的地步，這就有待讀者自行酙酌。

由此工具之協助，最後可由左下角之輸出按鈕輸出相關之適配方程式係數至工作空間內，讀者可自行摸索，本例之對應範圍為：

　
y=49.348+3.7361x-0.6149x2-0.010235x3

　　
其變化範圍上下限為：

　
Lower limit:
y=45.134+0.57715x-1.2665x2-0.049307x3

Higher limit:
y=53.562+6.895x+0.03667x2-0.028838x3

　

11.11.2 多項式迴歸polyfit & polyval

前文已介紹ployfit多項式適配之用法，現在再介紹其詳細的應用。polyfit指令係以多項式函數尋找配合Ｘ資料之曲線，主要在求其適配之多項式係數：

y(x)=p1xn+p2xn-1+...+pnx+pn+1

求係數與求預測值之語法分別如下：

[p,S] = polyfit(x,y,n)
Y = polyval(p,X)
[Y,DELTA] = polyval(p,X,S)
[Y,DELTA] = polyconf(p,X,S)
[Y,DELTA] = polyconf(p,X,S,alpha)

其中Ｘ為觀測值，n為模擬之階數。輸出參數p則為多項式之係數，S為特性矩陣，可用於polyval指令函數以產生預測值。若y資料為獨立常態性且定值變方，則利用polyval或polyconf所產生之誤差區間，或Ｙ DELTA至少包含50%之預測點。polyconf之輸入參數有alpha項，此為顯著水準，其預設值為0.05，此為設定信任水準求得相關對應值之指令。

例如：要適配一組隨機之變數樣本，利用normrnd來完成。normrnd之函數指令，其功能與randn指令相同，只是normrnd(m,s,M,N)之參數中，m為設定之均值，s為設定之標準差。

[p,S]=polyfit(1:20,[1:20]+normrnd(0,1,1,20),2)
p =
-0.0083    1.1649   -0.7796
S =
  R: [3x3 double]
 df: 17
normr: 4.0758

得到的方程式為

y(x)=-0.0083x²+1.1649x - 0.7796

使用上述得到之參數，可代入polyval進行預測：

X=2:0.5:15;
[Y, D]=polyval(p,X,S);
plot(X,Y,'r-',X,Y+D,'b:',X,Y-D,'b:'),
xlabel('X');ylabel('Y');grid

例二：汽車之資料包括車重，選樣測定其里程數，其中有部份不合格者，試求不合格之比例：

% A set of car weights
weight = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]';
% The number of cars tested at each weight
tested = [48 42 31 34 31 21 23 23 21 16 17 21]';
% The number of cars failing the test at each weight
failed = [1 2 0 3 8 8 14 17 19 15 17 21]';
% The proportion of cars failing for each weight
proportion = failed ./ tested;
plot(weight,proportion,'s')
xlabel('Weight'); ylabel('Proportion');
% We can try fitting a straight line to these data.
linearCoef = polyfit(weight,proportion,3);
linearFit = polyval(linearCoef,weight);
plot(weight,proportion,'s', weight,linearFit,'r-', ...
[2000 4500],[0 0],'k:', [2000 4500],[1 1],'k:')
xlabel('Weight'); ylabel('Proportion');

11.12通用線性模式

11.12.1通化線適配Glmfit

常用線性迴歸 (generalized linear model)可以用於直線之適配或任何能與其參數構成線性之函數，其資料之取得假設為常態分佈。雖然這個方法並不十分實際，但仍然算是最簡單的迴歸方式。通用線性模式將線性模式加以擴充，其一是引入連結函數，使參數之線性可以放宽；其二是誤差分佈即使非常態分配的狀態也可使用。

B = glmfit(X, Y, 'distr')
B = glmfit(X, Y, 'distr', 'link', 'estdisp', offset, pwts, 'const')
[B,dev,stats] = glmfit(...)

通化線模 GLMFIT指令係以預測矩陣Ｘ及迴應值Ｙ來配對，其採用的分配曲線由'distr'指定。執行結果置於Ｂ，是為估計之係數向量。 分配曲線'distr'可為 'normal'、 'binomial'、 'poisson'、 'gamma'與 'inverse gaussian'等。分配參數可適配至Ｘ中之各行，預設值是自行採用權威性連結。在大部份狀況下，Ｙ為回應之量測值；二項式分配曲線中，Ｙ為兩行資料，第一行為量測數目，第二行為試驗之數目（二項式中之Ｎ參數）。Ｘ之列數與Ｙ相同，包括每一觀察點之預測值。

參數'link'亦提供選擇性連結，以此提供f(mu) = xb關係式之修正，使分配參數mu與預測子xb可以產生較佳的適配性。'link'可為下列情況中之一種：

文字串 ：其內容及代表意義如下表：

'link'	意義	預設連結
'identity'	μ=xb	'normal'
'log'	log(μ)=xb	'poisson'
'logit'	log(μ/(1-μ))=xb	'bionomial'
'probit'	norminv(μ)=xb
'comploglog'	log(-log(1-μ))=xb
'logloglink'	log(-log(μ))=xb
'reciprocal'	1/μ=xb	'gamma'
p（一個數字）	μp=xb	'inverse gsussian'(p= -2)

-一個細胞陣列的型式如 {@FL @FD @FI} ，此三函數分別定義該連結 (FL)、連結之導數 (FD)及反連結 (FI)。
- 一個細胞陣列的型式包括三個 inline 函數以定義連結、導數及反連結函數。

上述後兩項連結中，其細胞陣列可使用外加函數或以inline指令設定之三函數──連結、其導數及反函數。例如：

本函數：　　FL=inline('x.^-0.5')
導數函數：　FD=inline('-0.5*x.^-1.5')
反函數：　　FI=inline('x.^-2')

然後在細胞陣列位置填入　{FL FD FI}
上述函數定義亦可改寫為：　{@FL @FD @FI}，只要將此三檔案定義在M-檔案中。或用匿函數如 FL=@(x) x.^-0.5

導數函數：　FD=@(x)  -0.5*x.^-1.5
反函數：　　FI=@(x) x.^-2

然後在細胞陣列位置填入　{FL FD FI}即可。

其他， 'estdisp'為離散參數之開關，設定為'on'時可使二項式或波義松分佈曲線之離散參數估計值( dispersion parameter)及標準差；設定為 'off' 時則使用理論離散值。對某些分配而言，'on'常為預設值。

輸入參數offset則是一個額外的預測向量，惟其係數固定為1.0。pwts則為先前權重之向量，例如每一對Ｘ與Ｙ的觀察頻率。參數 'const'可為 'on' (預設值) 以包含常數項，常數項省略時則設定為 'off' 。常數項為Ｂ之第一元素（不要直接在Ｘ矩陣中的第一行輸入１）。

輸出方面，dev為正解之偏異值，stats則為統計資料結構，包括下面之欄位： stats.dfe (誤差之自由度)、 stats.s (理論與離散參數估計值), stats.sfit (離散參數估計值)、stats.estdisp(=1有離散估計值；=0固定)stats.beta(Ｂ係數之向量) stats.se (B係數之標準差), stats.coeffcorr (B相關係數矩陣), stats.t (B之t 檢驗)、 stats.p (B之p-值), stats.resid (殘值)、stats.residp (Pearson 殘值)、stats.residd (偏異殘值)、stats.resida (Anscombe殘值)。

範例：汽車之資料中，依重量由2100-4300磅的範圍有很多的選擇，可以進行加侖里程之檢測，設認為性能不佳者(poor)依抽測之結果有如下之數據：

w = (2100:200:4300)';
poor = [1 2 0 3 8 8 14 17 19 15 17 21]';
total = [48 42 31 34 31 21 23 23 21 16 17 21]';
[w poor total]
ans =
2100           1          48
2300           2          42
2500           0          31
2700           3          34
2900           8          31
3100           8          21
3300          14          23
3500          17          23
3700          19          21
3900          15          16
4100          17          17
4300          21          21
plot(w,poor./total,'ro')

上述資料中，total為總抽測之數量；而 poor為認為性能不佳之數量。性能不佳之比例具有兩界限值，即[0 1]。可以利用這些資料進行適配分析，但一般線性迴歸無法包含兩個界限附近的數值，使用多項式或對數是一個方法。下面之例是利用'binomial'、 'logit'（預設）與使用'probit'進行適配比較：

[bl,dl, s1] = glmfit(w,[poor total],'binomial');
[bp,dp] = glmfit(w,[poor total],'binomial','probit');
[dl dp]

ans =
6.4842    7.5693

就偏異值而言，使用logit模式比probit為小。所以logit模式似乎略勝一疇。由上式的統計資料中，可以得知迴歸係數、標準差、t檢驗及Ｐ值分別如下：

[bl sl.se sl.t sl.p]

ans =
-13.3801    1.3940   -9.5986    0.0000
0.0042    0.0004    9.4474    0.0000

由於Ｐ值均接近於零，故其對應之截距與斜率均在顯著水準。

11.12.2 通化線預測值glmval

利用通化式預測值時，可用指令glmval進行計算，並與特定函數連結，其語法如下：

yfit = glmval(B,X,'link')
[yfit, dlo, dhi] = glmval(B,X,'link',  stats, clev)
[yfit, dlo, dhi] = glmval(B, X, 'link', stats, clev, N, offset, 'const')

這是配合通化線適配指令glmfit之輸出參數B及連結函數'link'，可以利用glmval指令，輸入預測值Ｘ，得到對應之新的觀測值yfit。Ｂ及連結函數'link'應與原glmfit指令使用的相同。yfit值為利用連結函數中之反函數依X*B的方式求得。

在輸入參數中，大體上其定義與glmfit指令之參數相同，Ｂ為利用glmfit求得，clev則為信任水準，其預設值為0.95。其上下限 [yfit-dlo, yfit+dhi]為信任水準之範圍，對應於特定之Ｘ值。此範圍僅適用於適配曲線，不能應用於新的觀察值。 N 為二項式Ｎ參數，配合glmfit之二項式分配曲線使用。offset及'const'的定義則與前面同。

就上述之汽車里程測定為例：

w = [2100:200:4300]';
poor = [1 2 0 3 8 8 14 17 19 15 17 21]';
total = [48 42 31 34 31 21 23 23 21 16 17 21]';
[b2,d2,s2] = glmfit([w w.^2],[poor total],'binomial');
b2'
wnew = (3000:100:4000)';
[yfit,dlo,dhi] = glmval(b2,[wnew wnew.^2],'logit',s2,0.95,30);
errorbar(wnew,yfit,dlo,dhi); 

ans =
  -7.3109    0.0002    0.0000

11.13 迴歸分析regress

11.13.1最小平方法多項式迴歸

多項式簡單迴歸可用線性解法直接求解：

 y=a0 + a1t +a2t2

若利用最小平方法，將模式與觀測點誤差之平方和求取最小值，將可間接求得上項係數。
迴歸分析是實驗工作常需使用之工具。其基本的迴歸為線性，以最小平方法解此線性模式：

y= Xβ+ε
ε~N(0,σ²I)

在此公式中，y為nx1觀察值，Ｘ為nxP之設計矩陣，其列數與y之行數相同，而p則為預測子之變數。β則為px1之向量參數，為待求之解；此外ε為隨機函數。其對應之指令型式如下：

B = regress(y, X)
[B, bint, r, rint, stats] = regress(y, X)
[B, bint, r, rint, stats] = regress(y, X, alpha)

指令 regress(y, X)為以最小平方法進行多重回歸，執行後會傳回線性模式Y = X*B之係數B向量。其中，Ｘ為nxp之設計參數，其列數n必須與觀測值Ｙ之行長度相同，且其第一行必須為1，以包含模式中之常數項。Ｙ則為nx1之觀測量向量。因為Ｆ檢定與Ｐ值之計算必須模式含有常數項。
輸入參數中，alpha為顯著水準，預設值為0.05。輸出值bint則為在此信任水準下，Ｂ之上下區間。r為殘值，其平方則為迴歸平方和與變方和之比。rint則為殘值r之上下區間，可作為偵察離散度之用。stats則為統計參數，含有R平方、Ｆ檢定及Ｐ值等。

在Ｘ矩陣內，行間是線性獨立的，regress指令執行時 ，會先設定Ｂ元素之最大數為零，以先得到基本解，其後再傳回bint，若Ｂ為零值時會傳回零。Ｘ與Ｙ中若含有NaN，則指令會視為遺失值，會自動去除。

例：在MATLAB 工具箱中，存有moore.mat之資料，其中前面五行為設定變數，第六行為回應值。求其回歸係數關係。

>>load moore;
moore
moore =
1.0e+003 *
  1.1250    0.2320    7.1600    0.0859    8.9050    0.0016
  0.9200    0.2680    8.8040    0.0865    7.3880    0.0009
  0.8350    0.2710    8.1080    0.0852    5.3480    0.0007
  1.0000    0.2370    6.3700    0.0838    8.0560    0.0007
  1.1500    0.1920    6.4410    0.0821    6.9600    0.0003
  0.9900    0.2020    5.1540    0.0792    5.6900    0.0004
  0.8400    0.1840    5.8960    0.0812    6.9320    0.0001
  0.6500    0.2000    5.3360    0.0806    5.4000    0.0001
  0.6400    0.1800    5.0410    0.0784    3.1770   -0.0002
  0.5830    0.1650    5.0120    0.0793    4.4610   -0.0002
  0.5700    0.1510    4.8250    0.0787    3.9010         0
  0.5700    0.1710    4.3910    0.0780    5.0020         0
  0.5100    0.2430    4.3200    0.0723    4.6650   -0.0001
  0.5550    0.1470    3.7090    0.0749    4.6420   -0.0002
  0.4600    0.2860    3.9690    0.0744    4.8400   -0.0004
  0.2750    0.1980    3.5580    0.0725    4.4790   -0.0002
  0.5100    0.1960    4.3610    0.0577    4.2000   -0.0002
  0.1650    0.2100    3.3010    0.0718    3.4100   -0.0004
  0.2440    0.3270    2.9640    0.0725    3.3600   -0.0005
  0.0790    0.3340    2.7770    0.0719    2.5990   -0.0000
X = [ones(size(moore,1),1) moore(:,1:5)];y = moore(:,6);
[b,bint,r,rint,stats] = regress(y,X);
[b bint],[r rint],stats

ans = % b        bint
 -2.1561   -4.1154   -0.1969
 -0.0000   -0.0011    0.0011
  0.0013   -0.0014    0.0040
  0.0001   -0.0000    0.0003
  0.0079   -0.0221    0.0379
  0.0001   -0.0000    0.0003
ans = % r      rint
  0.5623    0.2258    0.8988
 -0.1456   -0.5476    0.2565
  0.0885   -0.3262    0.5032
 -0.0479   -0.5515    0.4557
 -0.2307   -0.7043    0.2429
  0.1707   -0.2802    0.6216
 -0.3413   -0.8377    0.1550
 -0.0708   -0.6260    0.4844
 -0.0103   -0.4749    0.4543
 -0.1094   -0.6400    0.4211
  0.1717   -0.3311    0.6745
  0.0504   -0.4907    0.5915
 -0.0399   -0.5938    0.5140
  0.0227   -0.4991    0.5445
 -0.3945   -0.8701    0.0812
  0.0813   -0.4169    0.5795
  0.0730   -0.0879    0.2338
  0.0114   -0.4987    0.5214
 -0.2223   -0.6676    0.2231
  0.3806   -0.0071    0.7682
stats =
  0.8107   11.9886    0.0001    0.0685

由stats之值可知， ，表示此模式已經包括了80%之變化。Ｆ檢定值近於12，而Ｐ值則為0.0001，低於0.05之顯著水率，顯然非常不會有迴歸係數為零之狀況。最後一項有關誤差之估計值則為0.0685，也相當小。

11.13.2 步進迴歸stepwisefit

步進迴歸可以依自進變數之顯著性進行，篩選其指令為stepwisefit，語法如下：

b = stepwisefit(X, y)
[b,se,pval,inmodel,stats,nextstep,history] = stepwisefit(...)
[...] = stepwisefit(X, y, 'Param1',value1,'Param2',value2,...)
stepwise(X, y)
stepwise(X, y, inmodel, penter, premove)

步進迴歸stepwise是一個交談式之迴歸方式而stepwisefit則是直接迴歸曲線，並得到係數值。步進迴歸旨在利用設計之Ｘ參數矩陣預測Ｙ向量值，以此獲得迴歸係數b 。在交談式迴歸中，沒有指定之輸出參數，但等到交談過程中得到滿意的結果後，即可將相關係數輸出。在stepwisefit指令之輸入項中，其相關參數說明如下：

• b為計算得到之係數向量，se為其標準差。
• pval為Ｐ值向量，檢驗是否 b=0？
• inmodel為邏輯向量，其長度與Ｘ之行數相同，旨在說明留在模式之預測子。例如在j處之值為1表丕第j個預測值留在模式內；若為0時則否。
• stats為其他統計資訊。
• nextstep為建議下一步做法，說明那一個預測值要加入或移出。若結果為0時表示沒有意見。
• history為記錄每一步驟資訊結構。

在輸入項中，有'Param1',value1之輸入項，其內容包括如下：

參數名稱 參數值
'inmodel' 邏輯向量，說明那一個預測子在最初迴歸時納入。若全為零，表示沒有預測子。
'penter' 加入新預測子之最大Ｐ值，預設值為0.05。
'premove' 移去預測子之最小Ｐ值，預設值為0.01。
'display' 步驟顯示關關：'on'為顯示資訊；'off'不顯示。
'maxiter' 步驟最高次數（預設值為無最大值）
'keep' 邏輯向量，說明那一個預測子在最初迴歸時要保留。若全為零，表示沒有預測子。
'scale' 'on' 在迴歸前就Ｘ中各行以標準差進行分格；'off'不不格。

例：就工具箱hald.mat之資料進行迴歸分析

load hald
stepwisefit(ingredients, heat, 'penter', .08)

Initial columns included:  none
Step 1, added column 4, p=0.000576232
Step 2, added column 1, p=1.10528e-006
Step 3, added column 2, p=0.0516873
Step 4, removed column 4, p=0.205395
Final columns included:  1 2
ans =
  'Coeff'      'Std.Err.'    'Status'    'P'        
  [ 1.4683]    [  0.1213]    'In'        [2.6922e-007]
  [ 0.6623]    [  0.0459]    'In'        [5.0290e-008]
  [ 0.2500]    [  0.1847]    'Out'       [     0.2089]
  [-0.2365]    [  0.1733]    'Out'       [     0.2054]
ans =
  1.4683
  0.6623
  0.2500
 -0.2365

使用stepwise之交談式分析畫面：

11.13.3 非線性迴歸nlinfit

進行非線性迴歸時可使用nlinfit指令，其語法如下：

beta = nlinfit(X,y,fun,beta0)
[beta,r,J] = nlinfit(X,y,fun,beta0)
[...] = nlinfit(X, y, fun, beta0, options)

非線性迴歸法係利用最小平方法求得非線性函數之係數。其中y為回應向量，Ｘ為設計矩陣，配合所需之獨立變數，其值為一列對應一個Ｙ元素。實際上Ｘ可為任意陣列，只要fun這個函數接受。這個fun函數具有下列之型式：

yhat=myfun(beta, X)

其中，beta為係數向量，而Ｘ為設計矩陣。fun回應一個適配y值之向量yhat。在原來之指令中，beta0則為係數之初值。其正確之係數值beta則出現在等號之左方，其殘值為r，及Jacobian矩陣Ｊ。這些輸出參數都可以作為nlpredci函數指令之輸入，以估計預測值及係數值之相關誤差。
參數中，options則提供輸入的控制參數，以正確執行nlinfit之功能。此options亦可利用statset函數產生。options的欄位包括：

1. MaxIter：最高容許迴圈次數，預設值為100。
2. TolFun：殘餘平方運算終止容許值，預設值為1e-8。
3. TolX：計算beta係數之終止容許值，預設值為1e-8。
4. Display：運算期間之輸出－'off'（預設值）或 'iter'或 'final'

例：在工具箱中，有reaction.mat之資料，這是有關化學反應所到之觀察值，其中牽涉到三種化學反應：氫、n-戊烷、同位戊烷等。其中之函數hougen已經存在，是使用Hougen-Watson模式建立之反應動能量，其回應值為預設之反應速率。首先，下載資料：

load reaction
betafit = nlinfit(reactants,rate,@hougen,beta)
betafit =
   1.2526
   0.062776
   0.040048
   0.11242
 　1.1914

其中之hougen.m可以type hougen，得其內容如下：

function yhat = hougen(beta,x)
%HOUGEN Hougen-Watson model for reaction kinetics.
%   YHAT = HOUGEN(BETA,X) gives the predicted values of the
%   reaction rate, YHAT, as a function of the vector of
%   parameters, BETA, and the matrix of data, X.
%   BETA must have 5 elements and X must have three
%   columns.
%
%   The model form is:
%   y = (b1*x2 - x3/b5)./(1+b2*x1+b3*x2+b4*x3)
%
%   Reference:
%      [1]  Bates, Douglas, and Watts, Donald, "Nonlinear
%      Regression Analysis and Its Applications", Wiley
%      1988 p. 271-272.


%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 2.7.2.1 $  $Date: 2004/01/24 09:34:06 $
%   B.A. Jones 1-06-95.

b1 = beta(1);
b2 = beta(2);
b3 = beta(3);
b4 = beta(4);
b5 = beta(5);

x1 = x(:,1);
x2 = x(:,2);
x3 = x(:,3);


yhat = (b1*x2 - x3/b5)./(1+b2*x1+b3*x2+b4*x3);

另有一個nlintool指令函數則可利用圖形介面執行如下，讀者可試在指令窗下打入執行之。

曲線迴歸工具CFTOOL

對於任何組資料之迴歸，通常可以使用迴曲線迴歸工具cftool，以獲得所要的結果。此工具是採用圖形介面處理方式，故可直接審視資料並將資料與特定曲線作迴歸。利用顯示之圖形，可觀察其變化及與曲線間之適配情形，其中包括使用殘數分析及預測範圍線性等。

這個工具有示範程式，讀者可以自行上網取得。其語法如下：

• cftool
• cftool(xdata,ydata)

其中，xdata、ydata分別為：

xdata	預測值之向量資料
ydata	反應值之向量資料

CFtool(xdata, ydata)是打開此迴歸工具時，同時將一組特定待分析的資料作為輸入。這兩組資料是相對應的，故必須同大小。其中若包括inf與Nan等無效資料時，將會遭剔除。若使用複數資料格式輸入則僅實數部份被採用。

下圖所顯示為迴歸工具執行後，所得到的介面。其所用的資料為示範用的調查資料，置於census.mat。讀者可用load census.mat讀入，然後觀察比較：


在圖形之最上方有四組按鈕，可以選擇其執行功能：

1. Data ：讀取資料
2. Fitting ：回歸分析
3. Exclude ：捨棄某些資料點
4. Plotting：進行繪圖
5. Analysis：進行分析

資料輸入 GUI介面

資料GUI介面可輸入、預覽、取名或取消相關資料，作有效的整理，並去除一些不好的資料。其圖形介面如下：

資料迴歸GUI介面

利用迴歸介面可以產生參數型或非參數型之迴歸方程式，並檢查與比較不同迴歸結果，其間之係數及適配狀況(goodness of fit)之統計分析。在此節中所得之迴歸方程式及相關統計資料都會留下紀錄，並且可以利用程式再現。

下圖之迴歸介面為多項式四次回歸的結果：

隔除 GUI介面

隔除介面提供一個暫時性隔離資料點，使分析之結果更接近滿意的程度。利用隔除介面可以將特定資料點不參與分析，這些點可能獨立的個別點，也可以整區段的點。

下圖所示為頭兩點標示為被隔除點，其移除法則分別給定一個名稱，此處稱為excl。

繪圖 GUI介面

繪圖介面是將控制之資料組合進行繪製，下圖是將調查之資料進行適配後所顯示之結果，其迴歸之程式名為：poly2，如下圖：

分析過程

就本例而言，可以使用外插法作迴歸二次方程式之分析，以預測美國人口自2000人至2050年間每十年之變化。然後再就其結果進行繪圖，其過程如下：

• 在 Analyze at Xi 欄內輸入適當向量。
• 在 Evaluate fit at Xi 欄打勾
• 在 Plot results 及 Plot data set 欄上打勾。
• 按下 Apply 按鈕。

經過數值外插法計算結果如下：

經過外插法計算之結果與原先調查資料合併顯示如下：

儲存分析結果

按下 Save to workspace 按鈕，此時會將外插計算值以結構陣列之方式儲存於工作空間中。

其結果如下所示：

• analysisresults1
  
  analysisresults1 =
        xi: [6x1 double]
      yfit: [6x1 double]

下圖為針對2000至2050年間，以十年為區間之調查資料之預測與分析結果：