群组函数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