第5章〓统计性数据分析实战

统计工具箱提供了用于描述数据、分析数据以及为数据建模的函数和App。可以使用描述性统计量和绘图进行探索性数据分析，对数据进行概率分布拟合，生成进行蒙特卡罗仿真的随机数，以及执行假设检验。

5.1统计量和统计图

描述性统计分析要对调查总体所有变量的有关数据做统计性描述，主要包括数据的频数分析、数据的集中趋势分析、数据离散程度分析、数据的分布，以及一些基本的统计图形。常见的分析方法包括对比分析法、平均分析法、交叉分析法等。

5.1.1描述性统计量

描述性统计量是指通过生成汇总统计量，包括集中趋势、散度、形状和相关性方面的度量，以数值方式来探索数据。统计工具箱允许计算包含缺失(NaN)值的样本数据的汇总统计量。利用一元图、二元图和多元图，实现数据的可视化。可用选项包括箱线图、直方图和概率图。


1. 数据管理

在MATLAB中，可以使用多种不同的文件格式将数据导入和导出。有效格式包括表格数据、以制表符分隔的文件、Microsoft Excel电子表格以及SAS XPORT文件。统计工具箱还提供了更多数据类型，用于处理分组变量和分类数据。此工具箱还支持MATLAB中许多(但不是全部)可用的数据类型。

统计工具箱中包括表51中的实例数据集。要将数据集加载到MATLAB工作区中，可输入： 

load filename

其中，filename是表中列出的文件之一。

数据集包含单独的数据变量、具有引用的描述变量以及封装数据集及其描述的数据集数组(如果适用)。



表51常用的实例数据集



文件数据集的描述


acetylene.mat具有相关预测变量的化学反应数据
arrhythmia.mat来自 UCI 机器学习存储库的心律失常数据
carbig.mat汽车的测量值，1970—1982
carsmall.matcarbig.mat 的子集。汽车的测量值，1970、1976、1982
census1994.mat来自 UCI 机器学习存储库的成人数据
cereal.mat早餐谷物成分
cities.mat美国大都市地区的生活质量评分
discrim.mat用于判别分析的 cities.mat 版本
examgrades.mat0~100分的考试成绩
fisheriris.matFisher 1936 年的鸢尾花数据
flu.matGoogle 流感趋势估计的美国不同地区的 ILI(流感样疾病)百分比，疾病预防控制中心根据哨点提供商报告对 ILI 百分比进行了加权
gas.mat1993 年马萨诸塞州的汽油价格
hald.mat水泥发热与原料混合
hogg.mat牛奶的不同配送方式中的细菌数量
hospital.mat仿真的医疗数据
humanactivity.mat5种活动的人类活动识别数据： 坐、站、走、跑和跳舞
imports85.mat1985年来自UCI存储库的自动导入数据库
ionosphere.mat来自UCI机器学习存储库的电离层数据集
kmeansdata.mat四维聚类数据
lawdata.mat15所法学院的平均分数和LSAT分数
mileage.mat两家工厂的三种汽车型号的里程数据
moore.mat关于5个预测变量的生化需氧量
morse.mat非编码人员对摩尔斯电码的识别情况
nlpdata.mat从MathWorks文档中提取的自然语言处理数据
ovarianceancer.mat关于 4000 个预测变量的分组观测值
parts.mat36个圆形零件的大小偏差
polydata.mat多项式拟合的样本数据
popcorn.mat爆米花机型和品牌的爆米花产出
reaction.matHougenWatson 模型的反应动力学
sat.dat按性别和测验分列的学术能力测验 (SAT) 平均分(表)
sat2.dat按性别和测验分列的学术能力测验 (SAT) 平均分 (csv)
spectra.mat60份汽油样本的近红外光谱和辛烷值
stockreturns.mat仿真的股票回报







2. 数据类型

统计工具箱还另外提供了两种数据类型。要处理有序和无序的离散非数值数据，可以使用nominal和ordinal数据类型。要将多个变量(包括具有不同数据类型的变量)存储到一个对象中，可以使用dataset数组数据类型。但是，这些数据类型是统计和机器学习工具箱所独有的。要获得更好的跨产品兼容性，可分别使用MATLAB中提供的categorical或 table数据类型。

【例51】使用数据集数组变量及其数据。

(1) 按名称访问变量。

可以通过使用变量(列)名称和点索引来访问变量数据或选择变量子集。加载样本数据集数组，显示hospital中变量的名称。

>> load hospital 

hospital.Properties.VarNames(:)

ans =

7×1 cell 数组

{'LastName'}

{'Sex'  }

{'Age'  }

{'Weight'}

{'Smoker'}

{'BloodPressure'}

{'Trials'}


数据集数组中有7个变量(列)和100个观测值(行)。可以在工作区窗口中双击hospital 以在变量编辑器中查看数据集数组。

(2) 绘制直方图。

绘制变量Weight中数据的直方图，如图51所示。

>> figure

histogram(hospital.Weight)


图51中的直方图显示体重呈双峰分布。

(3) 绘制按类别分组的数据。

绘制按Sex中的值分组(男性和女性)的Weight的箱线图。也就是说，使用变量Sex作为分组变量。

>> figure   

boxplot(hospital.Weight,hospital.Sex)  %效果如图5-2所示



图51Weight中的数据直方图




图52箱线图




图52的箱线图表明性别是体重呈双峰分布的原因。

(4) 选择一个变量子集。

创建一个新数据集数组，其中仅包含变量LastName、Sex和Weight。可以通过名称或列号访问变量。

>> ds1 = hospital(:,{'LastName','Sex','Weight'});

ds2 = hospital(:,［1,2,4］);

数据集数组ds1和ds2是等同的。在对数据集数组进行索引时，使用括号()可保留数据类型； 也就是说，基于数据集数组的子集创建一个数据集数组。还可以使用变量编辑器基于变量和观测值的子集创建一个新数据集数组。


(5) 转换变量数据类型。

将变量Smoker的数据类型从逻辑值转换为名义值，标签为No和Yes。

>> hospital.Smoker = nominal(hospital.Smoker,{'No','Yes'});

class(hospital.Smoker)

ans =

'nominal'

(6) 探查数据。

显示Smoker的前10个元素。

>> hospital.Smoker(1:10)

ans = 

10×1 nominal 数组

Yes 

No 

No 

No 

No 

No 

Yes 

No 

No 

No

如果要更改名义数组中的水平标签，请使用setlabels。

(7) 添加变量。

变量BloodPressure是100×2数组。第一列对应于收缩压，第二列对应于舒张压。将此数组分成两个新变量SysPressure和DiaPressure。

>> hospital.SysPressure = hospital.BloodPressure(:,1);

hospital.DiaPressure = hospital.BloodPressure(:,2);

hospital.Properties.VarNames(:)

ans =

9×1 cell 数组

{'LastName'}

{'Sex'  }

{'Age'  }

{'Weight'}

{'Smoker'}

{'BloodPressure'}

{'Trials'}

{'SysPressure'  }

{'DiaPressure'  }

可见，数据集数组hospital有两个新变量。

(8) 按名称搜索变量。

使用regexp查找hospital中变量名称包含'Pressure'的变量。创建只包含这些变量的新数据集数组。

>> bp = regexp(hospital.Properties.VarNames,'Pressure');

bpIdx = cellfun(@isempty,bp);

bpData = hospital(:,~bpIdx);

bpData.Properties.VarNames(:)

ans =

3×1 cell 数组

{'BloodPressure'}

{'SysPressure'  }

{'DiaPressure'  }

可见，新数据集数组bpData仅包含血压变量。

(9) 删除变量。

从数据集数组hospital中删除变量BloodPressure。

>> hospital.BloodPressure = ［］;

hospital.Properties.VarNames(:)

ans =

8×1 cell 数组

{'LastName'  }

{'Sex'}

{'Age'}

{'Weight' }

{'Smoker' }

{'Trials' }

{'SysPressure'}

{'DiaPressure'}

可见，变量BloodPressure不再在数据集数组中。

5.1.2常用的统计量函数

根据样本数据计算描述性统计量，包括有关集中趋势、散度、形状、相关性和协方差的度量。制作数据的一般报表和交叉表，并计算分组数据的汇总统计量。如果数据中包含缺失(NaN)值，MATLAB算术运算函数将返回NaN。不过，统计工具箱提供的专用函数可以忽略这些缺失值，并返回使用其余值计算的数值。

下面介绍两个较为常用的统计量函数。

1. prctile函数

prctile函数用于计算数据集的百分位数。函数的调用格式如下。

Y=prctile(X,p)： 根据区间［0,100］中的百分比p返回数据向量或数组X中元素的百分位数。

 如果X是向量，则Y是标量或向量，向量长度等于所请求百分位数的个数(length(p))。Y(i)包含第p(i)个百分位数。

 如果X是矩阵，则Y是行向量或矩阵，其中，Y的行数等于所请求百分位数的个数 (length(p))。Y的第i行包含X的每一列的第p(i)个百分位数。

 对于多维数组，prctile在X的第一个非单一维度上进行运算。

Y=prctile(X,p,'all')： 返回X的所有元素的百分位数。

Y=prctile(X,p,dim)： 返回运算维度dim上的百分位数。

Y=prctile(X,p,vecdim)： 基于向量vecdim所指定的维度返回百分位数。例如，如果X 是矩阵，则prctile(X,50,［1 2］)返回X的所有元素的第50个百分位数，因为矩阵的每个元素都包含在由维度1和2定义的数组切片中。

Y=prctile(,'Method',method)： 使用上述任一语法中的输入参数组合，根据method的值，返回精确或近似百分位数。

【例52】计算数组中所有值的百分位数。

>> %创建3×5×2数组X

X = reshape(1:30,［3 5 2］)

X(:,:,1) =

1 4 71013

2 5 81114

3 6 91215

X(:,:,2) =

1619222528

1720232629

1821242730

>>%计算X的元素的第40个和第60个百分位数

Y = prctile(X,［40 60］,'all')

Y =

12.5000

18.5000

2. corr函数

corr函数用于计算线性或秩相关性。函数的调用格式如下。

rho=corr(X)： 返回输入矩阵X中各列之间的两两线性相关系数矩阵。

rho=corr(X,Y)： 返回输入矩阵X和Y中各列之间的两两相关系数矩阵。

［rho,pval］=corr(X,Y)： 还返回pval，它是一个p值矩阵，用于基于非零相关性备择假设来检验无相关性假设。

除了上述语法中的输入参数，［rho,pval］=corr(,Name,Value)还使用一个或多个名称值对组参数指定选项。例如，'Type','Kendall'指定计算Kendall tau相关系数。


【例53】计算两个矩阵之间的相关性，并将其与两个列向量之间的相关性进行比较。

%生成样本数据

>> rng('default')

X = randn(30,4);

Y = randn(30,4);

>> %在矩阵X的第二列和矩阵Y的第四列之间引入相关性

Y(:,4) = Y(:,4)+X(:,2);

%计算 X 和 Y 的列之间的相关性

rho,pval］ = corr(X,Y)

rho =

-0.1686-0.03630.22780.3245

0.30220.0332-0.08660.7653

-0.3632-0.0987-0.0200-0.3693

-0.1365-0.18040.08530.0279

pval =

0.37310.84890.22600.0802

0.10450.86190.64910.0000

0.04850.60390.91660.0446

0.47210.34000.65390.8837

5.1.3统计可视化

在MATLAB中，使用一元图(如箱线图和直方图)研究一元分布，使用二元图(如分组散点图和二元直方图)显示变量之间的关系，使用多元图(如Andrews图和图形符号图)可视化多个变量之间的关系。通过添加记录名称、最小二乘线条和参考曲线来自定义绘图。

下面介绍两个常用的函数。

1. boxplot函数

在统计工具箱中，提供了boxplot函数绘制箱线图。函数的调用格式如下。

boxplot(x)： 创建x中数据的箱线图。如果x是向量，boxplot绘制一个箱子； 如果x是矩阵，boxplot为x的每列绘制一个箱子。

在每个箱子上，中心标记表示中位数，箱子的底边和顶边分别表示第25个和75个百分位数。须线会延伸到不是离群值的最远端数据点，离群值会以'+'符号单独绘制。

boxplot(x,g)： 使用g中包含的一个或多个分组变量创建箱线图。boxplot为具有相同的一个或多个g值的各组x值创建一个单独的箱线。

boxplot(ax,)： 使用坐标区图形对象ax指定的坐标区和任何上述语法创建箱线图。

boxplot(,Name,Value)： 使用由一个或多个Name和Value对组参数指定的附加选项创建箱线图。例如，可以指定箱子样式或顺序。

函数的应用可参考例51。

2. refline函数

在统计工具箱中，提供了refline函数实现将参考线添加到绘图中。函数的调用格式如下。

refline(m,b)： 在当前坐标区中添加一条具有斜率m和截距b的参考线。

refline(coeffs)： 将由向量coeffs的元素定义的线添加到图窗中。

refline(ax,)： 使用上述任一语法中的输入参数，向ax所指定坐标区中的图上添加一条参考线。

hline=refline()： 使用上述任一语法中的输入参数，返回参考线对象hline。在创建参考线后，使用hline修改其属性。


【例54】在均值处添加参考线。

>>%为自变量x和因变量y生成样本数据

x = 1:10;

y = x + randn(1,10);

%创建x和y的散点图，如图5-3所示

scatter(x,y,25,'b','*')

下面实现在散点图上叠加一条最小二乘线，如图5-4所示

>> refline




图53x和y的散点图




图54叠加一条最小二乘线




下面实现在散点图的均值处添加一条参考线，如图55所示。

>> mu = mean(y);

hline = refline(［0 mu］);

hline.Color = 'r';



图55均值处添加一条参考线


5.2概 率 分 布

概率分布，是概率论中的基本概念之一，主要用于表述随机变量取值的概率规律，可用来计算均值和中值等汇总统计量、可视化样本数据、生成随机数等。在MATLAB中可使用概率分布对象、命令行函数或交互式App来处理概率分布。

概率分布可分为离散概率分布和连续概率分布，下面针对这两种分布进行介绍。

5.2.1离散概率分布

离散概率分布是指随机变量只能取有限(或可数无限)数量的值的概率分布。例如，在二项分布中，随机变量X只能取值0或1。统计工具箱提供了几种处理离散概率分布的方法，包括概率分布对象、命令行函数和交互式App。

1. 二项分布

二项分布即重复n次独立的伯努利实验。在每次实验中只有两种可能的结果，而且两种结果发生与否互相对立，并且相互独立，与其他各次实验结果无关，事件发生与否的概率在每一次独立实验中都保持不变，则这一系列实验总称为n重伯努利实验，当实验次数为1时，二项分布服从01分布。

在统计和机器学习工具箱中提供了如下几种处理二项分布的方法。

(1) 通过将概率分布拟合到样本数据(fitdist)或指定参数值(makedist)来创建概率分布对象BinomialDistribution。然后，使用对象函数来评估分布、生成随机数等。

(2) 使用Distribution Fitter应用程序以交互方式处理二项分布。可以从应用程序中导出对象并使用对象函数。

(3) 使用具有指定分布参数的分布特定函数(binocdf、binopdf、binoinv、binostat、binofit、binornd)。特定于分布的函数可以接受多个二项分布的参数。

(4) 使用具有指定分布名称(“二项式”)和参数的通用分布函数(cdf、icdf、pdf、随机)。

二项分布的概率密度函数(pdf)为： 

f(x|N,p)=Nxpx(1－p)N－x,x=0,1,2,…,N

其中，x是成功概率为p的伯努利过程的N次实验中的成功次数。结果是在N次实验中恰好x次成功的概率。对于离散分布，pdf也称为概率质量函数(pmf)。

二项分布的累积分布函数(cdf)为： 

F(x|N,p)=∑xi=0Nipi(1－p)N－i;x=0,1,2,…,N

其中，x是成功概率为p的伯努利过程的N次实验中的成功次数。结果是在N次实验中最多x次成功的概率。


下面对几个常用的二项分布函数进行介绍。

1) fitdist函数

在统计工具箱中，提供了fitdist函数对数据进行概率分布对象拟合。函数的调用格式如下。

pd=fitdist(x,distname)： 通过对列向量x中的数据进行distname指定的分布拟合，创建概率分布对象。

pd=fitdist(x,distname,Name,Value)： 使用一个或多个名称值对组参数指定的附加选项创建概率分布对象。例如，可以为迭代拟合算法指示删失数据或指定控制参数。

［pdca,gn,gl］=fitdist(x,distname,'By',groupvar)： 基于分组变量groupvar对x中的数据进行distname指定的分布拟合，以创建概率分布对象。它返回拟合后的概率分布对象的元胞数组pdca、组标签的元胞数组gn以及分组变量水平的元胞数组gl。

［pdca,gn,gl］=fitdist(x,distname,'By',groupvar,Name,Value)： 使用一个或多个名称值对组参数指定的附加选项返回上述输出参数。

2) pdf函数

在统计工具箱中，pdf函数为概率密度函数。函数的调用格式如下。 

y=pdf('name',x,A)： 返回由'name'和分布参数A指定的单参数分布族的概率密度函数 (pdf)，在x中的值处计算函数值。

y=pdf('name',x,A,B)： 返回由'name'以及分布参数A和B指定的双参数分布族的pdf，在x中的值处计算函数值。

y=pdf('name',x,A,B,C)： 返回由'name'以及分布参数A、B和C指定的三参数分布族的pdf，在x中的值处计算函数值。

y=pdf('name',x,A,B,C,D)： 返回由'name'以及分布参数A、B、C和D指定的四参数分布族的pdf，在x中的值处计算函数值。

y=pdf(pd,x)： 返回概率分布对象pd的pdf，在x中的值处计算函数值。

【例55】对数据进行正态分布拟合。

>> %加载样本数据。创建包含患者体重数据的向量

load hospital

x = hospital.Weight;

%通过对数据进行正态分布拟合来创建正态分布对象

pd = fitdist(x,'Normal')

pd = 



图56分布的pdf图

NormalDistribution

正态 分布


mu = 154  ［148.728, 159.272］

sigma = 26.5714 ［23.3299, 30.8674］


参数估计值旁边的区间是分布参数的95%置信区间。

%绘制分布的pdf，如图5-6所示

>> x_values = 50:1:250;

y = pdf(pd,x_values);

plot(x_values,y,'LineWidth',2)


3) cdf函数

在统计工具箱中，cdf函数为累积分布函数。函数的调用格式如下。

y=cdf('name',x,A)： 基于x中的值计算并返回由'name'和分布参数A指定的单参数分布族的累积分布函数(cdf)值。

y=cdf('name',x,A,B)： 基于x中的值计算并返回由'name'以及分布参数A和B指定的双参数分布族的cdf。

y=cdf('name',x,A,B,C)： 基于x中的值计算并返回由'name'以及分布参数A、B和C指定的三参数分布族的cdf。

y=cdf('name',x,A,B,C,D)： 基于x中的值计算并返回由'name'以及分布参数A、B、C和D指定的四参数分布族的cdf。

y=cdf(pd,x)： 基于x中的值计算并返回概率分布对象pd的cdf。

y=cdf(,'upper')： 使用可更精确计算极值上尾概率的算法返回cdf的补函数。'upper' 可以跟在上述语法中的任何输入参数之后。

【例56】计算正态分布cdf。

%创建均值μ等于0、标准差σ等于1的标准正态分布对象

>> mu = 0;

sigma = 1;

pd = makedist('Normal','mu',mu,'sigma',sigma);

>> %定义输入向量 x 以包含用于计算cdf的值

x = ［-2,-1,0,1,2］;

%基于 x 中的值计算标准正态分布的cdf值

y = cdf(pd,x)

y =

0.02280.15870.50000.84130.9772

y中的每个值对应于输入向量x中的一个值。例如，在值x等于1时，对应的cdf值y等于0.8413。

也可以不创建概率分布对象而直接计算同样的cdf值。使用cdf函数，再使用同样的μ和σ参数值指定一个标准正态分布。

>> y2 = cdf('Normal',x,mu,sigma)

y2 =

0.02280.15870.50000.84130.9772

cdf值与使用概率分布对象计算的值相同。

2. 多项式分布

多项式分布是二项分布的推广。二项分布(也叫伯努利分布)的典型例子是扔硬币，硬币正面朝上概率为p，重复扔n次硬币，k次为正面的概率即为一个二项分布概率。而多项分布就像扔骰子，有6个面对应6个不同的点数。二项分布时事件X只有两种取值，而多项分布的X有多种取值，多项分布的概率公式为： 

f(x|n,p)=n!x1！…xk！px11…pxkk

式中，k是每个实验中可能出现的相互排斥的结果数，n是实验总数。向量x=(x1…xk)是每k个结果的观察数，包含总和为n的非负整数分量。向量p=(p1…pk)是每k个结果的固定概率，包含总和为1的非负标量分量。

在n次实验中结果i的预期观察次数为： 

E{xi}=npi


其中，pi是结果i的概率。方差的结果i是： 

var(xi)=npi(1－pi)

结果i和j的协方差为： 

cov(xi,xj)=－npipj,i≠j

【例57】生成随机数，计算并绘制pdf，以及使用概率分布对象计算多项式分布的描述性统计信息。

(1) 定义分布参数。

创建一个包含每个结果概率的向量p。结果1的概率为1/2，结果2的概率为1/3，结果3的概率为1/6。每个实验的实验次数n为5次，重复次数为8次。

>> p = ［1/2 1/3 1/6］;

n = 5;

reps = 8;

(2) 创建一个多项式概率分布对象。

使用“概率”参数的指定值p创建多项式概率分布对象。

>> pd = makedist('Multinomial','Probabilities',p)

pd = 

MultinomialDistribution

Probabilities:

0.50000.33330.1667

>> rng('default')  % 重复性

r = random(pd)

r =

2

结果表明，本实验结果为2。

(3) 生成一个随机数矩阵。

还可以从多项式分布生成一个随机数矩阵，该矩阵报告了包含多个实验的多个实验结果。生成的矩阵，其中包含n=5次实验和8次重复的实验结果。

>> r = random(pd,reps,n)

r =

3 3 3 2 1

1 1 2 2 1

3 3 3 1 2

2 3 2 2 2

1 1 1 1 1

1 2 3 2 3

2 1 3 1 1

3 1 2 1 1



图57多项分布图

结果矩阵中的每个元素都是一次实验的结果。列对应于每个实验中的5个实验，行对应于8个实验。例如，在第一个实验中(对应于第一行)，5个实验中的一个得出结果1，5个实验中的一个得出结果2，5个实验中的三个得出结果3。

(4) 计算并绘制pdf。

计算分布的pdf，并绘图，如图57所示。

>> x = 1:3;

y = pdf(pd,x);

bar(x,y)

xlabel('结果')

ylabel('概率质量')

title('三项分布')


该图显示了每k个可能结果的概率质量。对于此分布，除1、2或3之外的任何x的pdf值均为0。

(5) 计算描述性统计。

计算分布的均值、中值和标准差。

>> m = mean(pd)   %均值

m =

1.6667

>> med = median(pd)  %中值

med =

1

>> s = std(pd)   %标准差

s =

0.7454

3. 泊松分布

泊松分布适用于涉及计算在给定的时间段、距离、面积等范围内发生随机事件的次数的应用情形。应用泊松分布的例子包括盖革计数器每秒咔嗒的次数、每小时走入商店的人数，以及每1000英尺录像带的瑕疵数。

泊松的概率分布函数为： 

f(x|λ)=λxx!e－λ;x=0,1,2,…,∞

泊松分布是接受非负整数值的单参数离散分布。参数λ既是分布的均值，也是分布的方差。因此，随着泊松随机数的特定样本中的数字变大，数字的变异性也变大。

泊松分布是二项分布的极限情况，其中，N趋向无穷大，p趋向零，而Np= λ。

泊松分布和指数分布是相关的。如果计数的数量遵循泊松分布，则单个计数之间的间隔遵循指数分布。

【例58】计算并绘制参数λ=5的泊松分布的pdf。

>> x = 0:15;

y = poisspdf(x,5);

plot(x,y,'+')

运行程序，效果如图58所示。


4. 离散型均匀分布

均匀分布是一种简单的概率分布，分为离散型均匀分布和连续型均匀分布。此处介绍的是离散型均匀分布。

离散型均匀分布是一个简单的分布，它对从1到N的整数赋予相等的权重。离散型均匀分布的概率公式为： 

y=f(x|N)=1NI(1,…,N)(x)

【例59】绘制离散型均匀分布cdf。

对于所有离散分布，cdf是一个阶跃函数。图59显示了N=10的离散型均匀分布cdf。

>> x = 0:10;

y = unidcdf(x,10);

figure;

stairs(x,y)

h = gca;

h.XLim = ［0 11］;



图58泊松分布的概率分布图




图59离散型均匀分布的cdf图




5.2.2连续分布

连续型随机变量X的分布函数是连续的，它对应的分布为连续分布。常用的连续分布有正态分布、均匀分布、指数分布、伽马分布、贝塔分布等。其中，正态分布是最常用的连续分布，如测量误差、人的身高、年降雨量等都可用正态分布描述。

本节对几个常用的连续分布展开介绍。

1. 正态分布

正态分布，有时也称为高斯分布，是双参数曲线族。使用正态分布建模的通常理由是中心极限定理，该定理(粗略地)指出，随着样本大小趋向无穷大，来自任何具有有限均值和方差的分布的独立样本总和会收敛为正态分布。

统计工具箱提供了以下几种处理正态分布的方法。

(1) 通过对样本数据进行概率分布拟合(fitdist)或通过指定参数值(makedist)来创建概率分布对象NormalDistribution。然后使用对象函数来计算分布、生成随机数等。

(2) 使用Distribution Fitter App以交互方式处理正态分布。

(3) 将分布特定的函数(normcdf、normpdf、norminv、normlike、normstat、normfit、normrnd)与指定的分布参数结合使用。分布特定的函数可以接受多个正态分布的参数。

(4) 将一般分布函数(cdf、icdf、pdf、random)与指定的分布名称('Normal')和参数结合使用。

1) 参数估计

最大似然估计(MLE)是最大化似然函数的参数估计。正态分布的μ和σ2的最大似然估计量分别是： 

x-=∑ni=1xin

和

s2MLE=1n∑ni=1(xi－x-)2


其中，x-是样本x1,x2,…,xn的样本均值。样本均值是参数μ的无偏估计量。但是，s2MLE是参数σ2的有偏估计量，这意味着其预期值不等于参数。

最小方差无偏估计量(MVUE)通常用于估计正态分布的参数。MVUE是参数的所有无偏估计量中方差最小的估计量。正态分布的参数μ和σ2的MVUE分别是样本均值x-和样本方差s2。

s2=1n－1∑ni=1(xi－x-)2

要对数据进行正态分布拟合并求出参数估计值，可使用normfit、fitdist或mle。

(1) 对于未删失数据，normfit和fitdist计算无偏估计值，mle计算最大似然估计值。

(2) 对于删失数据，normfit、fitdist、mle计算最大似然估计值。

与返回参数估计值的normfit和mle不同，fitdist返回拟合的概率分布对象 NormalDistribution。对象属性μ和σ存储参数估计值。

2) 概率密度函数

正态概率密度函数是： 

y=f(x|μ,σ)=1σ2πe－(x－μ)22σ2,x∈R

似然函数是被视为参数函数的pdf。最大似然估计(MLE)是最大化x的固定值的似然函数的参数估计。

3) 累积分布函数

正态累积分布函数表示为： 

p=F(x|μ,σ)=1σ2π∫x－∞－(t－μ)22σ2dt,x∈R

p是参数为μ和σ的正态分布中的一个观测值落入(－∞,x］区间的概率。标准正态累积分布函数Φ(x)在功能上与误差函数erf相关。

Φ(x)=121－erf－x2


其中，

erf(x)=2π∫x0e－t2dt=2Φ2x－1

【例510】拟合正态分布对象。

>>%加载样本数据并创建包含学生考试成绩数据的第一列的向量

load examgrades

x = grades(:,1);

%通过对数据进行正态分布拟合来创建正态分布对象

pd = fitdist(x,'Normal')

pd = 

NormalDistribution

正态 分布

mu = 75.0083［73.4321, 76.5846］

sigma =8.7202［7.7391, 9.98843］

参数估计值旁边的区间是分布参数的95%置信区间。

2. 指数分布

在概率理论和统计学中，指数分布(也称为负指数分布)是描述泊松过程中的事件之间的时间的概率分布，即事件以恒定平均速率连续且独立地发生的过程。指数分布与分布指数族的分类不同，后者是包含指数分布作为其成员之一的大类概率分布，也包括正态分布、二项分布、伽马分布、泊松分布等。

指数函数的一个重要特征是无记忆性(Memoryless Property，又称遗失记忆性)。这表示如果一个随机变量呈指数分布，当s,t≥0时，有P(T>s+t|T>t)=P(T>s)，即，如果T是某一元件的寿命，已知元件使用了t小时，它总共使用至少s+t小时的条件概率，与从开始使用时算起它使用至少s小时的概率相等。

指数分布的概率密度公式为： 

f(x)=λe－λx,x>0
0,x≤0

其中，λ> 0是分布的一个参数，常被称为率参数(rate parameter)，即每单位时间内发生某事件的次数。指数分布的区间是［0,∞)。如果一个随机变量X呈指数分布，则可以写作X~Exponential(λ)。

累积分布函数为： 

F(x;λ)=1－e－λx,x≥0
0,x<0

【例511】将指数分布拟合到数据。

>>%生成100个平均值为700的指数分布随机数样本

x = exprnd(700,100,1); %生成样本

%使用fitdist将指数分布拟合到数据

pd = fitdist(x,'exponential')

pd = 

ExponentialDistribution

指数 分布

mu = 661.084 ［548.486, 812.502］

fitdist返回一个指数分布对象。参数估计值旁边的区间是分布参数的95%置信区间。

>>%使用分布函数估计参数

［muhat,muci］ = expfit(x) %分布特定函数

muhat =

661.0843

muci =

548.4859

812.5023

>> ［muhat2,muci2］ = mle(x,'distribution','exponential') %通用分布函数

muhat2 =

661.0843

muci2 =

548.4859

812.5023

5.3假 设 检 验

统计工具箱中提供参数化假设检验和非参数化假设检验，帮助确定样本数据是否来自具有特定特征的总体。

(1) 分布检验(如AndersonDarling检验和单样本KolmogorovSmirnov检验)可以检验样本数据是否来自具有特定分布的总体。双样本KolmogorovSmirnov检验可以检验两组样本数据是否具有相同的分布。

(2) 位置检验(如z检验和单样本t检验)可以检验样本数据是否来自具有特定均值或中位数的总体。双样本t检验或多重比较检验可以检验两组或多组样本数据是否具有相同的位置值。

(3) 散度检验(如卡方方差检验)可以检验样本数据是否来自具有特定方差的总体。双样本F检验或多样本检验可以比较两个或多个样本数据集的方差。

5.3.1KS检验

KolmogorovSmirnov检验(KS检验)基于累积分布函数，用以检验一个经验分布是否符合某种理论分布或比较两个经验分布是否有显著性差异。

两样本KS检验由于对两样本的经验分布函数的位置和形状参数的差异都敏感而成为比较两样本有用且常规的非参数方法之一。

KS检验的累积分布函数为： 

Fn(x)=1n∑ni=1I－∞,x(Xi)

其中，I－∞,x为指示函数： 

I－∞,x(Xi)=
1,Xi≤x
0,Xi>x

对于一个样本集的累积分布函数Fn(x)和一个假设的理论分布F(x)，KS定义为： 

Dn=supxFn(x)－F(x)

supx是距离的上确界(supremum)，如果Xi服从理论分布F(x)，则当n趋于无穷时，Dn趋于0。

在MATLAB中，提供了kstest函数用来做单个样本的KolmogorovSmirnov检验； 它可以做双侧检验，检验样本是否服从指定的分布； 也可以做单侧检验，检验样本的分布函数是否在指定的分布函数之上或之下。

kstest函数的调用格式如下。

h=kstest(x)： 检验样本x是否服从标准正态分布，原假设是x服从标准正态分布，对立假设是x不服从标准正态分布。当输出h=1时，在显著性水平α=0.05下拒绝原假设； 当h=0时，则在显著性水平α=0.05下接受原假设。

h=kstest(x,CDF)： 检验样本x是否服从由CDF定义的连续分布。这里的CDF可以是包含两列元素的矩阵，也可以是概率分布对象，如ProbDistUnivParam类对象或ProbDistUnivKernel类对象。当CDF是包含两列元素的矩阵时，它的第1列表示随机变量的可能取值，可以是样本x中的值，也可以不是，但是样本x中的所有值必须在CDF的第1列元素的最小值与最大值之间。CDF的第2列是指定分布函数G(x)的取值。如果CDF为空(即［］)，则检验样本x是否服从标准正态分布。

h=kstest(x,CDF,alpha)： 指定检验的显著性水平alpha，默认值为0.05。

h=kstest(x,CDF,alpha,type)： 用type参数指定检验的类型(双侧或单侧)。type参数的可能取值如下。

(1) 当type='unequal'时即为双侧检验，对立假设是总体分布函数不等于指定的分布函数。

(2) 当type='larger'时为单侧检验，对立假设是总体分布函数大于指定的分布函数。

(3) 当type='smaller'时为单侧检验，对立假设是总体分布函数小于指定的分布函数。


其中，后两种情况下算出的检验统计量不用绝对值。

［h,p,ksstat,cv］=kstest(…)： 返回检验的p值、检验统计量的观测值ksstat和临界值cv。

【例512】在20天内，从维尼纶正常生活时的生产报表中看到的维尼纶纤度(纤维的粗细程度的一种度量)的情况，有如下100个数据。

1.36,1.49,1.43,1.41,1.37,1.40,1.32,1.43,1.47,1.39,

1.41,1.36,1.40,1.34,1.42,1.42,1.45,1.35,1.42,1.39,

1.44,1.42,1.39,1.42,1.42,1.30,1.34,1.42,1.37,1.36,

1.37,1.34,1.37,1.37,1.44,1.45,1.32,1.48,1.40,1.45,

1.39,1.46,1.39,1.53,1.36,1.48,1.40,1.39,1.38,1.40,

1.36,1.45,1.50,1.43,1.38,1.43,1.41,1.48,1.39,1.45,

1.37,1.37,1.39,1.45,1.31,1.41,1.44,1.44,1.42,1.42,

1.35,1.36,1.39,1.40,1.38,1.35,1.42,1.43,1.42,1.42,

1.42,1.40,1.41,1.37,1.46,1.36,1.37,1.27,1.37,1.38,

1.42,1.34,1.43,1.42,1.41,1.41,1.44,1.48,1.55,1.37.

试根据这100个样本数据在0.10显著性水平下，用KolmogorovSmirnov检验对维尼纶纤度数据进行正态性检验。(将数据保存到data.txt中。)

解析： 检验的原假设是维尼纶纤度服从正态分布。

其实现的MATLAB代码如下。

>> clear all;

X=load('data.txt');  %加载数据

［mu,sigma］=normfit(X)

x=(X-mu)/sigma;

［h,p,stats,cv］=kstest(x,［］,0.10,0)

运行程序，输出如下。

mu =

1.4038

sigma =

0.0474

h =

logical

0

p =

0.2951

stats =

0.2931

cv =

0.3687

结果表明，接受原假设，即认为维尼纶纤度服从均值为1.4038、标准差为0.0474的正态分布。

此外，MATLAB还提供了kstest2函数用来做两个样本的KolmogorovSmirnov检验，它可以做双侧检验，检验两个样本是否服从相同的分布，也可以做单侧检验，检验一个样本的分布函数是否在另一个样本的分布函数之上或之下。函数的调用格式如下。

h=kstest2(x1,x2)： 检验样本x1与x2是否具有相同的分布，原假设是x1与x2来自相同的连续分布，对立假设是来自于不同的连续分布。当输出h=1时，在显著性水平α=0.05下拒绝原假设； 当h=0时，则在显著性水平α=0.05下接受原假设。这里并不要求x1与x2具有相同的长度。

h=kstest2(x1,x2,alpha,type)： 指定检验的显著性水平alpha，默认值为0.05； 并用参数type指定检验的类型(双侧或单侧)。type参数的可能取值如下。

(1) 当type='unequal'时即为双侧检验，对立假设是两个总体的分布函数不相等。

(2) 当type='larger'时即为单侧检验，对立假设是第1个总体的分布函数大于第2个总体的分布函数。

(3) 当type='smaller'时即为单侧检验，对立假设是第1个总体的分布函数小于第2个总体的分布函数。

［h,p］=kstest2(…)： 返回检验的渐近p值，当p值小于或等于给定的显著性水平alpha时，拒绝原假设。样本容量越大，p值越精确，通常要求：

n1n2n1+n2≥4

其中，n1，n2分别为样本x1和x2的样本容量。

［h,p,ks2stat］=kstest2(…)： 返回检验统计量的观测值ks2stat。

【例513】利用kstest2函数对创建的标准正态随机分布检验是否接受原假设，并绘制其分布曲线图。

>> clear all;

x = -1:1:5;

y = randn(20,1);

［h,p,k］ = kstest2(x,y)

%以下代码用于绘制测试统计图

F1 = cdfplot(x);

hold on

F2 = cdfplot(y);

set(F1,'LineWidth',2,'Color','r')

set(F2,'LineWidth',2)

legend(［F1 F2］,'F1(x)','F2(x)','Location','NW')



图510测试统计图 


运行程序，输出如下，效果如图510所示。

h =

0

p =

0.2387

k =

0.4214

结果表明，由于h=0，所以在默认显著性下接受原假设。


5.3.2t检验

t检验，也称为student t检验(Students t test)，主要用于样本含量较小(例如n<30)、总体标准差σ未知的正态分布资料。

单样本t检验是当总体标准差未知时位置参数的参数化检验。检验统计量的计算公式为： 

t=x-－μs/n


其中，x-是样本均值，μ是假设的总体均值，s是样本标准差，n是样本大小。在原假设下，检验统计量具有Student t分布和n－1个自由度。

在MATLAB中，提供了ttest函数用于实现单样本t检验。函数的调用格式如下。

h=ttest(x)： 使用单样本t检验返回原假设的检验决策，该原假设假定x中的数据来自均值等于零且方差未知的正态分布。备择假设是总体分布的均值不等于零。如果检验在5%的显著性水平上拒绝原假设，则结果h为1，否则为0。

h=ttest(x,y)： 使用配对样本t检验返回针对原假设的检验决策，该原假设假定x－y中的数据来自均值等于零且方差未知的正态分布。

h=ttest(x,y,Name,Value)： 返回配对样本t检验的检验决策，其中使用由一个或多个名称值对组参数指定附加选项。

h=ttest(x,m)： 返回针对原假设的检验决策，该原假设假定x中的数据来自均值为m且方差未知的正态分布。备择假设是均值不为m。

h=ttest(x,m,Name,Value)： 返回单样本t检验的检验决策，其中使用一个或多个名称值对组参数指定附加选项。

［h,p］=ttest()： 还使用上述语法组中的任何输入参数返回检验的p值。

［h,p,ci,stats］=ttest()： 还返回x(对于配对t检验则为x－y)的均值的置信区间ci，以及包含检验统计量信息的结构体stats。

【例514】某种电子元件的寿命X(以小时计)服从正态分布，μ、σ2均未知。现测得16只元件的寿命如下。

160278198200236257270167150250194224137

185167255

问是否有理由认为元件的平均寿命大于220(小时)?

其实现的MATLAB代码如下。

>> clear all;

X=［160278198200236257270167150250194224137185167255］;

［h,p,ci］=ttest(X,220,0.005,1)

运行程序，输出如下。

h =

0

p =

0.8457

ci =

174.4465Inf  %均值225在该置信区间内

结果表明，h=0表示在水平α=0.05下应该接受原假设h0，即认为元件的平均寿命不大于220小时。

5.3.3双样本t检验

根据方差齐与不齐两种情况，应用不同的统计量进行检验。

方差不齐时，检验统计量为： 

t=－x-－y-S2xm+S2yn


式中，x-和y-表示样本1和样本2的均值； S2x和S2y为样本1和样本2的方差； m和n为样本1和样本2的数据个数。

方差齐时，检验统计量为： 

t=－x-－y-Sw1m+1n

式中，Sw为两个样本的标准差，它是样本1和样本2的方差的加权平均值的平方根，为： 

SW=(m－1)S2x+(n－1)S2ym+n+1


在不假设两个数据样本来自具有方差齐性的总体的情况下，原假设下的检验统计量具有近似Student t分布，其自由度的数目由Satterthwaite逼近给出。此检验有时称为Welch t检验。

在MATLAB中，提供了ttest2函数用于实现双样本t检验。函数的调用格式如下。

h=ttest2(x,y)： 使用双样本t检验返回原假设的检验决策，该原假设假定向量x和y中的数据来自均值相等、方差相同但未知的正态分布的独立随机样本。备择假设是x和y中的数据来自均值不相等的总体。如果检验在5%的显著性水平上拒绝原假设，则结果h为1，否则为0。

h=ttest2(x,y,Name,Value)： 返回针对双样本t的检验决策，该检验使用由一个或多个名称值对组参数指定的附加选项。例如，可以更改显著性水平或进行无须假设方差齐性的检验。

［h,p］=ttest2()： 还使用上述语法中的任何输入参数返回检验的p值。

［h,p,ci,stats］=ttest2()： 还返回总体均值差的置信区间ci，以及包含检验统计量信息的结构体stats。

【例515】下面分别给出文学家马克·吐温的8篇小品文以及斯诺德格拉斯的10篇小品文中的3个字母组成的单词的比例。

马克·吐温0.225 0.262 0.217 0.240 0.230 0.229 0.235 0.217

斯诺德格拉斯0.209 0.205 0.196 0.210 0.202 0.207 0.224 0.223 0.220 0.201

设两组数据分别来自正态总体，且两总体方差相等，但参数均未知。两样本相互独立，问两个作家所写的小品文中包含由3个字母组成的单词的比例是否有显著差异？零假设为两个作家对应的比例没有显著差异。

其实现的MATLAB代码如下。

>> clear all;

x=［0.225 0.262 0.217 0.240 0.230 0.229 0.235 0.217］;

y=［0.209 0.205 0.196 0.210 0.202 0.207 0.224 0.223 0.220 0.201］;

［h,signnificance,ci］=ttest2(x,y)

运行程序，输出如下。

h =

1

signnificance =

0.0013

ci =

0.01010.0343

结果表明，h=1，拒绝零假设，认为两个作家所写小品文中包含3个字母组成的单词的比例有显著差异。

5.4方 差 分 析

方差分析(ANOVA)是指确定响应变量的变异是出现在总体组内还是出现在不同总体组之间的过程。统计工具箱提供了单因素/双因素/N因素方差分析(ANOVA)、多元方差分析(MANOVA)、重复测量模型以及协方差分析(ANCOVA)。

5.4.1方差的基本原理

方差分析的基本原理是认为不同处理组的均数间的差别基本来源有以下两个。

(1) 随机误差，如测量误差造成的差异或个体间的差异，称为组内差异，用变量在各组的均值与该组内变量值的偏差平方和的总和表示，记作SSw，组内自由度dfw。

(2) 实验条件，即不同的处理造成的差异，称为组间差异。用变量在各组的均值与总均值之偏差平方和表示，记作SSb，组间自由度dfb。

总偏差平方和SSt=SSw+SSb。

组内SSw、组间SSb除以各自的自由度(组内dfw=n－m，组间dfb=m－1，其中，n为样本总数，m为组数)，得到其均方MSw和MSb，一种情况是处理没有作用，即各组样本均来自同一总体，MSb/MSw≈1。另一种情况是处理确实有作用，组间均方是由于误差与不同处理共同导致的结果，即各样本来自不同总体。那么，MSb>MSw。

MSb/MSw比值构成F分布。用F值与其临界值比较，推断各样本是否来自相同的总体。

5.4.2单因素方差分析

单因素方差分析是指对单因素实验结果进行分析，检验因素对实验结果有无显著性影响的方法。它是用来研究一个控制变量的不同水平是否对观测变量产生了显著影响。这里，由于仅研究单个因素对观测变量的影响，因此称为单因素方差分析。

例如，分析不同施肥量是否给农作物产量带来显著影响，考察地区差异是否影响妇女的生育率，研究学历对工资收入的影响等。这些问题都可以通过单因素方差分析得到答案。

单向方差分析是线性模型的一个简单特例。模型的单向方差分析形式为： 

yij=aj+εij


其中，yij是一个观察值，i代表观察值，j代表预测变量y的不同组(水平)。所有yij都是独立的。aj代表第j组(水平)的总体平均值。εij是独立的正态分布随机误差，具有零均值和常数方差，即εij~N(0,σ2)。

这个模型也称为均值模型。该模型假设y列为常数aj加上误差分量εij。方差分析有助于确定这些常数是否都相同。

方差分析检验了所有组均数相等的假设，以及至少一个组与其他组不同的替代假设。


H0: a1=a2=…=ak
H1: 并非所有组为均组


方差分析基于所有样本总体均为正态分布的假设。众所周知，它对适度违反这一假设具有鲁棒性。可以使用法线图(normclot)直观地检查法线假设。或者，可以使用统计和机器学习工具箱中的一个函数来检查正态性： AndersonDarling检验(adtest)、卡方拟合优度检验(chi2gof)、JarqueBera检验(jbtest)或Lilliefors检验(Lilliefest)。

在MATLAB中，提供了anova1函数实现单因素方差分析。anova1主要是比较多组数据的均值，然后返回这些均值相等的概率，从而判断这一因素是否对实验指标有显著影响。函数的调用格式如下。

p=anova1(X)： 为零假设存在的概率，一般p小于0.05或0.01时，认为结果显著(零假设可疑)。

p=anova1(X,group)： 当X为矩阵时，利用group变量作为X中样本箱线图的标签。

p=anova1(X,group,displayopt)： displayopt为on时，则激活anova1表和箱线图的显示。

［p,table］=anova1(…)： 返回单元数组表中的anoval表。

［p,table,stats］=anova1(…)： 返回stats结构，用于多元比较检验。

【例516】有A、B、C、D、E、F这6个小麦品种产量的比较实验，设置标准品种CK，采用3次重复的对比设计，所得产量结果如表52所示。


表526个小麦品种产量结果分析



品种
各重复的产量/kg

ⅠⅡⅢ


A560582520
B
582
565
525
C
600
600
572
D
525
496
590
E
560
578
615
F
640
662
508
CK
500
510
519


利用anova1实现方差分析，代码如下。

>> clear al;

X=［560 582 600 525 560 640 500;582 565 600 496 578 662 510;520 525 572 590  615 508 519］;

group={'A','B','C','D','E','F','CK'};

［p,table,stats］=anova1(X,group)

运行程序，输出如下，效果如图511及图512所示。

p =

0.1602



图511ANOVA表




图512方差盒形图


5.4.3双因素方差分析

双因素方差分析法是一种统计分析方法，这种分析方法可以用来分析两个因素的不同水平对结果是否有显著影响，以及两因素之间是否存在交互效应。一般运用双因素方差分析法，先对两个因素的不同水平的组合进行设计实验，要求每个组合下所得到的样本的含量都是相同的。

在实际问题的研究中，有时需要考虑两个因素对实验结果的影响。例如饮料销售，除了关心饮料品牌，还想了解销售地区是否影响销售量，如果在不同的地区，销售量存在显著的差异，就需要分析原因。采用不同的销售策略，使该饮料品牌在市场占有率高的地区保持领先地位； 在市场占有率低的地区进一步扩大宣传。如果把饮料的品牌看作影响销售量的因素A，饮料的销售地区则是影响因素B。对因素A和因素B同时进行分析，就属于双因素方差分析的内容，双因素方差分析是对影响因素进行检验： 究竟是一个因素在起作用，还是两个因素都起作用，或是两个因素的影响都不显著。

1. 双因素方差分析法的类型

双因素方差分析有以下两种类型。

(1) 一个是无交互作用的双因素方差分析，它假定因素A和因素B的效应之间是相互独立的，不存在相互关系。

(2) 另一个是有交互作用的双因素方差分析，它假定因素A和因素B的结合会产生出一种新的效应。

例如，如果假定不同地区的消费者对某种颜色有与其他地区消费者不同的特殊偏爱，这就是两个因素结合产生的新效应，属于有交互作用的背景； 否则，就是无交互作用的背景。

2. 双因素方差的模型

双向方差分析是线性模型的特例。模型的双向方差分析形式为： 

yijr=μ+ai+βj+(aβ)ij+εij


式中，

 yijr是对响应变量的观察。i表示行因子A的组i，i=1,2,…,I； j表示列因子B的组j，j=1,2,…,J； r表示复制数，r=1,2,…,R，即总共有N=I×J×R个观测值。

 μ为总体平均值。

 ai是由于行因子B导致的行因子A中的组与总体平均值μ的偏差。ai的值的和为0。

∑Ii=1ai=0

 βj是由于行因子B导致的列因子B中的组与总体平均值μ的偏差。给定列βj中的所有值均相同，且βj的值总和为0。

∑Jj=1βj=0


 (aβ)ij是由于行因子B导致的列因子B中的组与总体平均值μ的偏差。βj的给定列中的所有值都是相同的，并且βj的值总和为 0。

∑Ii=1(aβ)ij=∑Jj=1(aβ)ij=0

 εij是随机干扰。假设它们是独立的、正态分布的，并且具有恒定的方差。

在MATLAB统计工具箱中，提供了anova2函数用于实现无交互作用的双因素方差分析。函数的调用格式如下。

p=anova2(X,reps)： 根据样本观测值矩阵X进行均衡实验的双因素一元方差分析。X的每一列对应因素A的一个水平，每行对应因素B的一个水平，X还应满足方差分析的基本假定。reps表示因素A和B的每一个水平组合下重复实验的次数。

p=anova2(X,reps,displayopt)： 当displayopt为on时，则显示方差分析表和箱线图。

［p,table］=anova2(…)： 返回单元数组表中的ANOVA表。

［p,table,stats］ = anova2(…)： 返回stats结构，用于多元检验。

【例517】执行双向ANOVA以确定汽车型号和工厂对汽车里程等级的影响。

>>%加载并显示示例数据

load mileage

mileage

mileage =

33.300034.500037.4000

33.400034.800036.8000

32.900033.800037.6000

32.600033.400036.6000

32.500033.700037.0000

33.000033.900036.7000

返回的结果中，有三个车型(列)和两个工厂(行)。该数据有六个里程行，因为每个工厂为研究提供了每种型号的三辆汽车(即复制数为三)。第一个工厂的数据在前三行，第二个工厂的数据在后三行，进行双向方差分析。

>> nmbcars = 3;%每个型号的汽车数量，即复制次数

［~,~,stats］ = anova2(mileage,nmbcars);%效果如图5-13所示



图513ANOVA表


执行多重比较以找出三种车型中哪一对有显著差异。

>> c = multcompare(stats)%效果如图5-14所示

注意： 模型中包含一个交互效应项。当模型中包括交互效应时，主效应检验可能很难解释。


c =

1.00002.0000-1.5865-1.0667-0.54690.0004

1.00003.0000-4.5865-4.0667-3.54690.0000

2.00003.0000-3.5198-3.0000-2.48020.0000


在矩阵c中，前两列显示了比较的汽车模型对，最后一列显示了检验的p值。所有p值都很小(0.0004、0和0)，这表明所有车型的平均里程彼此之间存在显著差异。

在图514中，蓝色线为第一款车型平均里程的比较区间，红色条是第二款和第三款车型平均里程的比较区间。第二和第三比较区间均不与第一比较区间重叠，表明第一车型的平均里程与第二和第三车型的平均里程不同。如果单击其他栏之一，可以测试其他车型。如果没有一个比较区间重叠，表明每个车型的平均里程与其他两个有显著差异。



图514多重比较界面


5.4.4多因素方差分析

多因素方差分析，用于研究一个因变量是否受到多个自变量(也称为因素)的影响，它检验多个因素取值水平的不同组合之间、因变量的均值之间是否存在显著的差异。多因素方差分析既可以分析单个因素的作用(主效应)，也可以分析因素之间的交互作用(交互效应)，还可以进行协方差分析，以及各个因素变量与协变量的交互作用。

多因素方差分析是两因素方差分析的一般形式，对三个因素的情况，其模型表达式为： 

yijkl=μ+α.j.+βi..+γ..k+(αβ)ij.+(αγ)i.k+(βγ).jk+(αβγ)ijk+εijkl

式中两个连在一起的标记，如(αβ)ij.，表示两个因素之间的交互作用，参数(αβγ)ijk表示三个因素之间的交互作用。

在MATLAB中，提供了anovan函数用于实现多因素方差分析。函数的调用格式如下。

p=anovan(y,group)： 根据样本观测值向量y进行均衡或非均衡实验的多因素一元方差分析，检验多个因素的主效应是否显著。输入参数group为一个元胞数组，它的每一个元素对应一个因素，是该因素的水平列表，与y等长，用来标记y中每个观测所对应的因素的水平。每个元胞中因素的水平列表可以是一个分类(categorical)数组、数值向量、字符矩阵或单列的字符串元胞数组。输出参数p是检验的p值向量，p中的每个元素对应一个主效应。

p=anovan(y,group,param,val)： 通过指定一个或多个成对出现的参数名与参数值来控制多因素一元方差分析。

［p,table］=anovan(y,group,param,val)： 同时返回元胞数组形式的方差分析表table(包含列标签和行标签)。

［p,table,stats］=anovan(y,group,param,val)： 同时返回一个结构体变量stats，用于进行后续的多重比较。当某因素对实验指标的影响显著时，在后续的分析中，可以调用multcompare函数，把stats作为它的输入，进行多重比较。

［p,table,stats,terms］=anovan(y,group,param,val)： 同时返回方差分析计算中的主效应项和交互效应项矩阵terms。terms的格式与'model'参数的最后一种取值的格式相同。当'model'参数的取值为一个矩阵时，anovan函数返回的terms就是这个矩阵。

【例518】显示了如何对1970—1982年生产的406辆汽车的里程和其他信息的汽车数据进行N向方差分析。

%加载样本数据

>> load carbig

该实例侧重于四个变量，MPG是406辆汽车每加仑行驶的英里数，其他三个变量是因素： cyl4(是否为四缸汽车)、org(汽车起源于欧洲、日本或美国)，以及何时(汽车在该时期的早期、中期或期末)。

%拟合完整模型，要求最多三向交互和类型3平方和

>> varnames = {'Origin';'4Cyl';'MfgDate'};

anovan(MPG,{org cyl4 when},3,3,varnames) %效果如图5-15所示

ans =

0.0000

NaN

0.0000

0.7032

0.0001

0.2072

0.6990



图515N因素方差分析表


请注意，许多术语用#符号标记为没有满秩，其中一个术语的自由度为零，缺少p值。当缺少因子组合且模型具有高阶项时，可能会发生这种情况。在这种情况下，下面的交叉表显示，在这一时期的早期，除了四缸，没有欧洲制造的汽车，如tbl(2,1,1)中的0所示。

>> ［tbl,chi2,p,factorvals］ = crosstab(org,when,cyl4)

tbl(:,:,1) =

827525

0 4 3

3 3 4

tbl(:,:,2) =

122238

232617

122532

chi2 =

207.7689

p =

8.0973e-38

factorvals =

3×3 cell 数组

{'USA'   }{'Early'}{'Other'   }

{'Europe'}{'Mid'  }{'Four'}

{'Japan' }{'Late' }{0×0 double}

由结果可以看出，不可能估计三向相互作用效应，并且在模型中包含三向相互作用项会使拟合奇异。即使使用ANOVA表中有限的可用信息，也可以看到三向交互作用的p值为0.699，因此不显著。

只检查双向互动，如图516所示。



图516约束平方和方差表



>> terms

terms =

100

010

001

110

101

011

现在所有的项都是可估计的。相互作用项4(原点×4Cyl)和相互作用项6 (6Cyl×MfgDate)的p值远大于典型的临界值0.05，表明这些项不重要。可以选择忽略这些项，将它们的影响合并到误差项中。输出terms变量返回一个代码矩阵，每个代码都代表一个术语的位模式。通过从术语中删除条目来从模型中省略术语。

>> ［~,~,stats］ = anovan(MPG,{org cyl4 when},terms,3,varnames)  %得到方差分析表与图5-16一致

stats = 

包含以下字段的 struct:

source: 'anovan'

resid: ［1×406 double］

coeffs: ［30×1 double］

Rtr: ［14×14 double］

rowbasis: ［14×30 double］

dfe: 384

…

现在有一个更简洁的模型，表明这些汽车的里程似乎与所有三个因素有关，并且制造日期的影响取决于汽车的制造地点。对Origin和Cylinder执行多重比较。

>> results = multcompare(stats,'Dimension',［1,2］)  %多重比较表如图5-17所示

results =

1.00002.0000-7.5624-4.0331-0.50380.0144

1.00003.0000-8.4566-4.01740.42180.1024

1.00004.0000-10.3771-8.7025-7.02780.0000

1.00005.0000-14.0054-12.3561-10.70690.0000

1.00006.0000-12.8187-11.1673-9.51580.0000

2.00003.0000-5.46980.01575.50121.0000

2.00004.0000-8.3570-4.6693-0.98170.0042

2.00005.0000-11.9795-8.3230-4.66660.0000

2.00006.0000-10.7977-7.1341-3.47060.0000

3.00004.0000-9.3273-4.6850-0.04280.0464

3.00005.0000-12.9050-8.3387-3.77240.0000

3.00006.0000-11.6841-7.1498-2.61560.0001

4.00005.0000-5.5949-3.6537-1.71250.0000

4.00006.0000-4.4210-2.4648-0.50860.0045

5.00006.0000-0.76521.18893.14300.5092



图517多重比较表