11.1 前言
統計的技巧與資料分析常常形影不離。一般統計使用加法、累加法、平均值,中間值等等,由於處理的對象是矩陣資料,故其基本統計之技巧已經廣為應用,其觀念也會在正常之運作中出現。統計學中比較特殊應用者為機率、亂數、常態分配等,而配合應用者為其相關之圖表。在MATLAB中,有一個統計學工具箱,內藏各種統計學上需要應用的指令,可以執行上述與統計學有關之內容。這些相關的指令大部份以M-檔案組成,所以可利用type 這個功能檢視其內容。甚至可以更改其檔案名稱與內容,增加自己需要的功能,使其成為新的指令。此外,有些指令尚搭配繪圖介面,因而可以在繪圖模式下,進行資料與圖之適配,形成具體的方程式或實驗式,以供未來研究者使用。
統計工具箱中提供約二百餘個指令檔案,其中對機率分配方面則提供廿餘種機率型態,每種均有其相關的函數,諸如:
- 機率密度函數(pdf)
- 累積分佈函數(cdf)
- 累積分佈函數之反函數
- 亂數產生器
- 均數與變異數
張貼者: Martin Fon 位於 此文章的連結 Chap11
11.2.1平均值MEAN
統計學上對資料處理常用趨中的處理,求取均值或中間值等,均會取中的特徵。求取一個矩陣或向量之平均值時可用指令MEAN,其格式如下:
M=mean(A,dim)
若A為向量,其結果M為單一值,亦即向量中各元素之平均;若A為矩陣,則結果M為一列向量,其中元素為各行之平均值。dim為方向性參數,其預設值為1,表示結果係行向平均,故M為列向量;若dim=2,則係沿列向平均,結果M為行向量。例如:
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
張貼者: Martin Fon 位於 此文章的連結 Chap11
11.2.2中間值median
中間值亦可利用平均值的指令型式求得,其正式指令名稱為median。但其求得之值若非正好中間值,則會以接近中間值之兩值加以平均,其結果與mean之平均值仍有不同,下面以前述之B矩陣為例,比較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
由上述這些值看來,有些的結果相同,有些則因其分離的情況而有異。中間值及修剪平均值是不管偏異的大小,一般為避免有極端的數值的影響,採用這種平均法。
張貼者: Martin Fon 位於 此文章的連結 Chap11
11.2.3最大值與最小值 MIN & MAX
要求一矩陣或向量中元素之最大與最小值時,其指令之型式如下:
C = max|min(A)
C = max|min(A,B)
C = max|min(A,[],dim)
[C,I] = max(...)
若A為向量,其結果C為單一值,亦即向量中各元素之最大或最小;若A為矩陣,則結果C為一列向量,其中元素為各行之最大或最小。dim為方向性參數,其預設值為1,表示結果係行向取得最大或最小,故C為列向量;若dim=2,則係沿列向操作,結果M為行向量。注意要dim之參數時,需加在第三位置。此外,在輸出項中,I表示最大或最小元素之位置,不過此項功能僅求最大值時適用。例如:
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
在這個指令中,比較特殊的是兩個矩陣A與B之最大或最小比較,其結果C應為與A或B相同的矩陣,但比較A與B中對應元素間之最大與最小。例如:
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
張貼者: Martin Fon 位於 此文章的連結 Chap11
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則為加權值向量,此權值必須為正值,且長度與X相同。
若X為矩陣時,通常預設為行向計算,但可以利用dim=2參數改為以列向為計算基礎,其結果為行向量。var指令會將其元素除以總和,因此權值總和為一。若w值為零,其結果如var(X);若為1則如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
張貼者: Martin Fon 位於 此文章的連結 Chap11
11.2.5標準差STD
標準誤差為各樣本品與平均值間之常態差,其值實際上為上述變方var執行結果之開方值,其計算公式如下:
張貼者: Martin Fon 位於 此文章的連結 Chap11
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為母數,則可增加參數1,即採用 COV(X,1) 或 COV(X,Y,1)指令之型式。
張貼者: Martin Fon 位於 此文章的連結 Chap11
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矩陣,其內有X陣列行間之相關係數。而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)}
輸入項中除XY等資料矩陣外,尚可輸入其他特定變數與常數。這些可以用 'PARAM1',VAL1成對表示,其項目包括:
參數值:
'alpha' 顯著水準,預設值為0.05(即95%信任區間)
'rows' 使用 'all' (預設值)表示使用所有列值;
'complete'表示使用沒有含NaN 值之列;
'pairwise'表示計算R(i,j)時使用不含
NaN值之 i行或 j行。
輸出值中, P表示檢驗無關係假設之P值矩陣。每一個P值代表隨機可以觀察得到之最大值域。若 P(i,j)值很小,例如小於 0.05,則R(i,j) 之關係甚為顯著。
此外,有RLO與RUP代表95%信任水準之下限與上限矩陣,其大小與R相同。
例一
>> 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
張貼者: Martin Fon 位於 此文章的連結 Chap11
11.2.8群組函數grpstats
前面討論到之平均值求法,通常應用於整個陣列之值,若要應用到比較複雜的分組平均問題,則必須使用不同的函數才能達成。此項指令之格式如下:
means = grpstats(X, group)
[means, sem, counts, name] = grpstats(X, group, whichstats)
grpstats(x, group, alpha)
輸入參數中X為求平均值之對象,X可為多行,其平均結果也會多行。group則為與X同列長之陣列,可能由多項分組之向量組成,其內容可為字串列或細胞陣列之文字,如{G1 G2 G3}。若X中之元素同屬分組中之一項,則其平均值會出現在該項下。
在輸出項中,第一項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
張貼者: Martin Fon 位於 此文章的連結 Chap11
11.2.9表列函數tabulate
TABLE = tabulate(x)
這個指令可將一向量X之觀測值製成一表格,其第一行為X向量中之相同數值,第二行為該數值出現之次數,最後一行為該值出現之百分比。若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))
張貼者: Martin Fon 位於 此文章的連結 Chap11
11.2.10百分值prctile
在一般大量樣本之情況下,可以利用百分值去確定樣本之合理對應值,由此百分比與對應值之關係可以瞭解資料之外形、位置以及擴散度。其指令格式如下:
Y = prctile(X, p)
此指令計算X之樣本值中一個大於p%部份之對應值位置,此值並不一定是原有之觀測值,只求其比例位置。輸入參數 p 必須落在[0 100]間,可為常數或向量。若 p = 50% 時,則Y值應對應X之中間值(median)。X之資料可為向量或矩陣,而 p 則可能為一向量或其中之常數。下表說明可能之之四種狀況:
X p Y= prctile(X, p)
==== ======== ===========================================
向量 數值常數 數值等於X中第p級百分值
矩陣 數值常數 列向量含X中每一行內第p級百分值,其中
y(i)=prctile(X(:,j),p)
向量 向量 Y向量與p同,其第i項為X中第p(i)級百分值
y(i)=prctile(X,p(i))
矩陣 向量 Y矩陣中,其第(i,j)項為X中第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之向量相同,與X為列或行向量無關。
%X2若為矩陣,則p先與X之行向量作百段分值。
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)先與X之行向量作百段分值,並Y依序完成行向量之值。
張貼者: Martin Fon 位於 此文章的連結 Chap11
11.2.11細分值quantile
細分值與百分值之意義類似,但其區間為小數,介於[0 1]之間,以配合累積密度函數之使用,其指令格式如下:
Y = quantile(X, p)
Y = quantile(X, p, dim)
其輸出值Y為X觀測值中傳回值,p為數值或累積機率值之向量。當X為向量時,Y之大小與p相同,而Y(i)則是第p(i)之細分對應值。當X為矩陣,則Y之第i列為第p(i)對X之行向量之細分值。
其細分方向亦可利用dim設定,但Y在dim指定的方向長度應與p之長度相同。
細分值常應用於累加機率,故其範圍為[0 1]。若X為一具N元素之向量,則QUANTILE依下列方式運算:
- 1)X值經過排序,並細分為(0.5/N), (1.5/n), ..., ((N-0.5)/N) 細分段
- 2) 線性內插法用於計算 (0.5/N) 與 ((N-0.5)/N)間機率之細分段。
- 3) X中之最大值與最小值作為機率外圍之細段。
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
張貼者: Martin Fon 位於 此文章的連結 Chap11
堆疊矩陣之使用,前面也曾述及,其相關語法如下:
B = repmat(A,m,n)
B = repmat(A,[m n])
B = repmat(A,[m n p...])
這是一個處理大矩陣且內容有重複時使用之。其功能是以A之內容堆疊在一(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。A也可以置放文字串,如:
>>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
張貼者: Martin Fon 位於 此文章的連結 Chap11
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則有正負。