11/28/2006

11.2.8 群組函數grpstats

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