11.1 前言
統計的技巧與資料分析常常形影不離。一般統計使用加法、累加法、平均值，中間值等等，由於處理的對象是矩陣資料，故其基本統計之技巧已經廣為應用，其觀念也會在正常之運作中出現。統計學中比較特殊應用者為機率、亂數、常態分配等，而配合應用者為其相關之圖表。
在MATLAB中，有一個統計學工具箱，內藏各種統計學上需要應用的指令，可以執行上述與統計學有關之內容。這些相關的指令大部份以M-檔案組成，所以可利用type 這個功能檢視其內容。甚至可以更改其檔案名稱與內容，增加自己需要的功能，使其成為新的指令。此外，有些指令尚搭配繪圖介面，因而可以在繪圖模式下，進行資料與圖之適配，形成具體的方程式或實驗式，以供未來研究者使用。
統計工具箱中提供約二百餘個指令檔案，其中對機率分配方面則提供廿餘種機率型態，每種均有其相關的函數，諸如：
- 機率密度函數(pdf)
- 累積分佈函數(cdf)
- 累積分佈函數之反函數
- 亂數產生器
- 均數與變異數
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
中間值亦可利用平均值的指令型式求得，其正式指令名稱為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
要求一矩陣或向量中元素之最大與最小值時，其指令之型式如下：
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
變方值為各樣本品與平均值差之平方和。又稱為變值常態檢定公式，其指令型式為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
標準誤差為各樣本品與平均值間之常態差，其值實際上為上述變方var執行結果之開方值，其計算公式如下：
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
兩個變數相關性可由相關係數求得。其指令型式如下：
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
前面討論到之平均值求法，通常應用於整個陣列之值，若要應用到比較複雜的分組平均問題，則必須使用不同的函數才能達成。此項指令之格式如下：
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
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
在一般大量樣本之情況下，可以利用百分值去確定樣本之合理對應值，由此百分比與對應值之關係可以瞭解資料之外形、位置以及擴散度。其指令格式如下：
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
細分值與百分值之意義類似，但其區間為小數，介於[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
堆疊矩陣之使用，前面也曾述及，其相關語法如下：
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
統計分析之技巧中，通常以觀察值之平均值為參考標，經處理後再進行判定。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.1 前言
統計的技巧與資料分析常常形影不離。一般統計使用加法、累加法、平均值，中間值等等，由於處理的對象是矩陣資料，故其基本統計之技巧已經廣為應用，其觀念也會在正常之運作中出現。統計學中比較特殊應用者為機率、亂數、常態分配等，而配合應用者為其相關之圖表。在MATLAB中，有一個統計學工具箱，內藏各種統計學上需要應用的指令，可以執行上述與統計學有關之內容。這些相關的指令大部份以M-檔案組成，所以可利用type 這個功能檢視其內容。甚至可以更改其檔案名稱與內容，增加自己需要的功能，使其成為新的指令。此外，有些指令尚搭配繪圖介面，因而可以在繪圖模式下，進行資料與圖之適配，形成具體的方程式或實驗式，以供未來研究者使用。
統計工具箱中提供約二百餘個指令檔案，其中對機率分配方面則提供廿餘種機率型態，每種均有其相關的函數，諸如：
- 機率密度函數(pdf)
- 累積分佈函數(cdf)
- 累積分佈函數之反函數
- 亂數產生器
- 均數與變異數
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
中間值亦可利用平均值的指令型式求得，其正式指令名稱為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
要求一矩陣或向量中元素之最大與最小值時，其指令之型式如下：
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
變方值為各樣本品與平均值差之平方和。又稱為變值常態檢定公式，其指令型式為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
標準誤差為各樣本品與平均值間之常態差，其值實際上為上述變方var執行結果之開方值，其計算公式如下：
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
兩個變數相關性可由相關係數求得。其指令型式如下：
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
前面討論到之平均值求法，通常應用於整個陣列之值，若要應用到比較複雜的分組平均問題，則必須使用不同的函數才能達成。此項指令之格式如下：
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
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
在一般大量樣本之情況下，可以利用百分值去確定樣本之合理對應值，由此百分比與對應值之關係可以瞭解資料之外形、位置以及擴散度。其指令格式如下：
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
細分值與百分值之意義類似，但其區間為小數，介於[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
堆疊矩陣之使用，前面也曾述及，其相關語法如下：
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
統計分析之技巧中，通常以觀察值之平均值為參考標，經處理後再進行判定。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則有正負。
最低0.47元/天 解锁文章
7498
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
