分析方法
01方差分析
定义：用数理统计分析试验结果、鉴别各因素对结果影响程度的方法称为方差分析（Analysis Of Variance），记作 ANOVA。
我们关心的试验结果称为：指标
试验中需要考察、可以控制的条件称为：因素
因素所处的状态称为：水平
根据因素数目的不同可以划分为单因素方差分析和双因素方差分析。
01.1单因素方差分析
只考虑一个因素A对所关心的指标的影响， A取几个水平，在每个水平上作若干个试验，试验过程中除 A 外其它影响指标的因素都保持不变，我们的任务是从试验结果推断，因素 A 对指标有无显著影响，即当 A 取不同水平时指标有无显著差别。
单因素方差分析表：
Matlab函数：Matlab 统计工具箱中单因素方差分析的命令是：anoval
若各组数据个数相等，称为均衡数据。若各组数据个数不等，称非均衡数据。
例子(e01)：为考察 5 名工人的劳动生产率是否相同，记录了每人 4 天的产量，并算出其平均值，如表。你能从这些数据推断出他们的生产率有无显著差别吗？
x=[256 254 250 248 236
242 330 277 280 252
280 290 230 305 220
298 295 302 289 252];
p=anova1(x);
alpha=0.05;
if p>alpha
fprintf('with p=%.4f,we can accept H0\n',p);
else
fprintf('with p=%.4f,we can not accept H0\n',p);
end结果：
>> e01
with p=0.1109,we can accept H0
非均衡数据用法：p=anova1(x,group)。x为向量，从第1组到第r组数据依次排列；group为与x同长度的向量，标志x中数据的组别（在与x第i组数据相对应的位置处输入整数i=1,2,3,..,r）。
例子(e02)：用4种工艺生产灯泡，从各种工艺制成的灯泡中各抽出了若干个测量其寿命，结果如下表，试推断这几种工艺制成的灯泡寿命是否有显著差异。
x=[1620 1580 1460 1500
1670 1600 1540 1550
1700 1640 1620 1610
1750 1720 1680 1800];
x=[x(1:4),x(16),x(5:8),x(9:11),x(12:15)];
g=[ones(1,5),2*ones(1,4),3*ones(1,3),4*ones(1,4)];
p=anova1(x,g);
alpha=0.05;
if p>alpha
fprintf('with p=%.4f>alpha(0.05),we can accept H0\n',p);
else
fprintf('with p=%.4f<=alpha(0.05),we can not accept H0\n',p);
end结果：
>> e02
with p=0.0331<=alpha(0.05),we can not accept H0
01.2双因素方差分析
如果要考虑两个因素A,B对指标的影响,A,B各划分几个水平，对每一个水平组合作若干次试验，对所得数据进行方差分析，检验两因素是否分别对指标有显著影响，或者还要进一步检验两因素是否对指标有显著的交互影响。
如果根据经验或某种分析能够事先判定两因素之间没有交互影响，每组试验就不必重复，过程大为简化。无交互双因素方差分析表：
关于交互效应的双因素方差分析表：
Matlab求解双因素方差分析：p=anova2(x,reps)
其中 x 不同列的数据表示单一因素的变化情况，不同行中的数据表示另一因素的变化情况。如果每种行—列对（“单元”）有不止一个的观测值，则用参数 reps 来表明每个“单元”多个观测值的不同标号，即 reps 给出重复试验的次数 t 。下面的矩阵中，列因素有 3 种水平，行因素有两种水平，但每组水平有两组样本，相应地用下标来标识：
例子(e03)：一种火箭使用了四种燃料、三种推进器，进行射程试验，对于每种燃料与每种推进器的组合作一次试验，得到试验数据如下表。
问各种燃料之间及各种推进器之间有无显著差异？
解：记燃料为因素A，它有4个水平，水平效应为𝑎𝑖,i=1,2,3,4。推进器为因素B,它有 3 个水平，水平效应为𝛽𝑗 ， j=1,2,3。我们在显著性水平 a=0.05下检验:
x=[58.2 56.2 65.3
49.1 54.1 51.6
60.1 70.9 39.2
75.8 58.2 48.7];
[p,t,st]=anova2(x);
alpha=0.05;
if sum(p>alpha)>1
fprintf('with p=%.4f,%.4f,we can accept H0\n',p);
else
fprintf('with p=%.4f,%.4f,we can not accept H0\n',p);
end结果：
>> e03
with p=0.4491,0.7387,we can accept H0
例子(e04)：一火箭使用了4种燃料,3种推进器作射程试验，每种燃料与每种推进器的组合各发射火箭2次，得到如下表结果。试在水平0.05 下，检验不同燃料（因素A）、不同推进器（因素B）下的射程是否有显著差异？交互作用是否显著？
x0=[58.2,52.6 56.2,41.2 65.3,60.8
49.1,42.8 54.1,50.5 51.6,48.4
60.1,58.3 70.9,73.2 39.2,40.7
75.8,71.5 58.2,51.0 48.7,41.4];
x1=x0(:,1:2:5);
x2=x0(:,2:2:6);
x=zeros(8,3);
for i=1:4
x(2*i-1,:)=x1(i,:);
x(2*i,:)=x2(i,:);
end
[p,t,st]=anova2(x,2);
alpha=0.05;
if sum(p>alpha)>2
fprintf('with p=%.4f,%.4f,%.4f,we can accept H0\n',p);
else
fprintf('with p=%.4f,%.4f,%.4f,we can not accept H0\n',p);
end结果：
>> e04
with p=0.0035,0.0260,0.0001,we can not accept H0
02回归分析
一元线性回归：
多元线性回归：
求解：最小二乘法
Matlab求解：b=regress(Y,X)或：[b,bint,r,rint,stats]=regress(Y,X,alpha)
式中：X,Y为给定数据；b为回归系数估计值beta；bint为回归系数估计值的置信区间；r，rint为残差及其置信区间；stats是用于检验回归模型的统计量。残差去及其置信区间可以用rcoplot(r,rint)绘图。
例子(e05)：合金的强度 y 与其中的碳含量 x 有比较密切的关系，今从生产中收集了一批数据如下表。试先拟合一个函数 ) (x y ，再用回归分析对它进行检验。
%% 绘制散点图观察可能的曲线关系
x=0.1:0.01:0.18;
y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0];
plot(x,y,'r*')
%% 通过观察散点图，可假设函数关系为线性关系：y=ax+b
x1=(0.1:0.01:0.18)';
y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0]';
x=[ones(9,1),x1];
[b,bint,r,rint,stats]=regress(y,x);
b,bint,stats,figure,rcoplot(r,rint)
%% 结果分析
% 由p=0.0012<0.05知道在0.05的水平下拒绝原假设（H0：a=0），说明上述假设成立
% 观察残差图发现除了第8组数据外残差置信区间均包含零点，因此可认为第8组数据是异常的
% 剔除后再进行计算
x1=[0.1:0.01:0.16,0.18]';
y=[42,41.5,45.0,45.5,45.0,47.5,49.0,50.0]';
x=[ones(8,1),x1];
[b,bint,r,rint,stats]=regress(y,x);
b,bint,stats,figure,rcoplot(r,rint)结果：
>> e05
b =
27.4722
137.5000
bint =
18.6851 36.2594
75.7755 199.2245
stats =
0.7985 27.7469 0.0012 4.0883
b =
30.7820
109.3985
bint =
26.2805 35.2834
76.9014 141.8955
stats =
0.9188 67.8534 0.0002 0.8797
练习(p01)：某厂生产的一种电器的销售量 y 与竞争对手的价格x1和本厂的价格x2有关。下表是该商品在10个城市的销售记录。试根据这些数据建立y与x1,x2的关系式，对得到的模型和系数进行检验。若某市本厂产品售价 160（元），竞争对手售价 170（元），预测商品在该市的销售量。
x1=[120 140 190 130 155 175 125 145 180 150]';
x2=[100 110 90 150 210 150 250 270 300 250]';
y=[102 100 120 77 46 93 26 69 65 85]';
%% 绘制散点图观察
figure,plot(y,x1,'r*');
hold on,plot(y,x2,'b+');hold off;
%% 由散点图发现x2与y有较明显的线性关系，但是x1与y就不太明显
%首先使用线性假设进行尝试：y=ax1+bx2+c
x=[ones(10,1),x1,x2];
[b,bint,r,rint,stats]=regress(y,x);
b,bint,stats
%% 结果分析
% p=0.0247，在0.05的假设水平下可用。
结果：
>> p01
b =
66.5176
0.4139
-0.2698
bint =
-32.5060 165.5411
-0.2018 1.0296
-0.4611 -0.0785
stats =
0.6527 6.5786 0.0247 351.0445
多项式回归分析
如果从数据的散点图上发现 y 与 x 呈较明显的二次（或高次）函数关系，或者用线性模型的效果不太好，就可以选用多项式回归
一元多项式回归：Matlab中使用polyfit实现一元多项式回归
例子(e06)：将 17 至 29 岁的运动员每两岁一组分为7组，每组两人测量其旋转定向能力，以考察年龄对这种运动能力的影响。现得到一组数据如表。试建立二者之间的关系。
%% 一元多项式回归
x0=17:2:29;x0=[x0,x0];
y0=[20.48 25.13 26.15 30.0 26.1 20.3 19.35 24.35 28.11 26.3 31.4 26.92 25.7 21.3];
%% 绘制曲线图观察
plot(x0,y0,'r*');
%% 有曲线发现是较明显的拱形，可假设为二次函数关系式
[p,s]=polyfit(x0,y0,2);
p %输出结果是计算得到的拟合参数
[y,delta]=polyconf(p,x0,s)% 获取y的拟合值和预测值的置信区间半径
%% 使用polytool函数来可视化拟合曲线及其拟合半径，并可以获取预测值
polytool(x0,y0,2); 
多元二项式回归：Matlab中提供rstool工具箱来作相关图
用法： rstool(x,y,model,alpha)
式中：x，y是输入数据；alpha是显著性检验水平；model可选以下值：
练习(p01)：改进解法，使用rstool工具箱求解。某厂生产的一种电器的销售量 y 与竞争对手的价格x1和本厂的价格x2有关。下表是该商品在10个城市的销售记录。试根据这些数据建立y与x1,x2的关系式，对得到的模型和系数进行检验。若某市本厂产品售价 160（元），竞争对手售价 170（元），预测商品在该市的销售量。
%% 使用多元二项式求解
x1=[120 140 190 130 155 175 125 145 180 150]';
x2=[100 110 90 150 210 150 250 270 300 250]';
y=[102 100 120 77 46 93 26 69 65 85]';
x=[x1 x2];
%% 可以通过工具箱调整结果并导出数据
%beta(回归系数)，rmse（剩余标准差），residuals(残差)
rstool(x,y,'purequadratic') % 选择model为purequadratic结果：
03层次分析
层次分析法（Analytic Hierarchy Process，简称 AHP）是对一些较为复杂、较为模糊的问题作出决策的简易方法，它特别适用于那些难于完全定量分析的问题
例子(e08)：挑选合适的工作。经双方恳谈，已有三个单位表示愿意录用某毕业生。该生根据已有信息建立了一个层次结构模型，如下图所示：
(e08)构造判断矩阵：层次结构反映了因素之间的关系，但准则层中的各准则在目标衡量中所占的比重并不一定相同，在决策者的心目中，它们各占有一定的比例。根据这些比重，构建出判断矩阵。
(e08)准则层判断矩阵
(e08)方案层判断矩阵
(e08)层次总排序
(e08)根据层次总排序权值，该生最满意的工作为工作 1
clc,clear
fid=fopen('data.txt','r');
n1=6;n2=3;
a=[];
for i=1:n1
tmp=str2num(fgetl(fid));
a=[a;tmp]; %读准则层判断矩阵
end
for i=1:n1
str1=char(['b',int2str(i),'=[];']);
str2=char(['b',int2str(i),'=[b',int2str(i),';tmp];']);
eval(str1);
for j=1:n2
tmp=str2num(fgetl(fid));
eval(str2); %读方案层的判断矩阵
end
end
ri=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45]; %一致性指标
[x,y]=eig(a);
lamda=max(diag(y));
num=find(diag(y)==lamda);
w0=x(:,num)/sum(x(:,num));
cr0=(lamda-n1)/(n1-1)/ri(n1)
for i=1:n1
[x,y]=eig(eval(char(['b',int2str(i)])));
lamda=max(diag(y));
num=find(diag(y)==lamda);
w1(:,i)=x(:,num)/sum(x(:,num));
cr1(i)=(lamda-n2)/(n2-1)/ri(n2);
end
cr1, ts=w1*w0, cr=cr1*w0data.txt
1 1 1 4 1 1/2
1 1 2 4 1 1/2
1 1/2 1 5 3 1/2
1/4 1/4 1/5 1 1/3 1/3
1 1 1/3 3 1 1
2 2 2 3 3 1
1 1/4 1/2
4 1 3
2 1/3 1
1 1/4 1/5
4 1 1/2
5 2 1
1 3 1/3
1/3 1 1/7
3 7 1
1 1/3 5
3 1 7
1/5 1/7 1
1 1 7
1 1 7
1/7 1/7 1
1 7 9
1/7 1 1
1/9 1 104判别分析
在生产、科学研究和日常生活中，经常会遇到对某一研究对象属于哪种情况作出判断。例如要根据这两天天气情况判断明天是否会下雨；医生要根据病人的体温、白血球数目及其它症状判断此病人是否会患某种疾病等等。Matlab 的统计工具箱提供了判别函数 classify。
调用格式为：[CLASS,ERR] = CLASSIFY(SAMPLE,TRAINING,GROUP, TYPE)
其中：SAMPLE 为未知待分类的样本矩阵；
TRAINING 为已知分类的样本矩阵，它们有相同的列数 m ，设待分类的样本点的个数，即 SAMPLE 的行数为 s ，已知样本点的个数，即 TRAINING 的行数为 t ，则 GROUP 为 t 维列向量。
TYPE 为分类方法，缺省值为'linear'，即线性分类，TYPE 还可取值'quadratic'，'mahalanobis'（mahalanobis 距离）。
返回值 CLASS 为 s 维列向量，给出了 SAMPLE 中样本的分类
ERR 给出了分类误判率的估计值
例子(e09)：已知 8 个乳房肿瘤病灶组织的样本，其中前 3 个为良性肿瘤，后 5 个为恶性肿瘤。数据为细胞核显微图像的10个量化特征：细胞核直径，质地，周长，面积，光滑度。根据已知样本对未知的三个样本进行分类。已知样本的数据为：
| 13.54 | 14.36 | 87.46 | 566.3 | 0.09779 | 13.54 | 14.36 | 87.46 |
| 13.08 | 15.71 | 85.63 | 520 | 0.1075 | 13.08 | 15.71 | 85.63 |
| 9.504 | 12.44 | 60.34 | 273.9 | 0.1024 | 9.504 | 12.44 | 60.34 |
| 17.99 | 10.38 | 122.8 | 1001 | 0.1184 | 17.99 | 10.38 | 122.8 |
| 20.57 | 17.77 | 132.9 | 1326 | 0.08474 | 20.57 | 17.77 | 132.9 |
待分类数据为：
| 16.6 | 28.08 | 108.3 | 858.1 | 0.08455 |
| 20.6 | 29.33 | 140.1 | 1265 | 0.1178 |
| 7.76 | 24.54 | 47.92 | 181 | 0.05263 |
%% 层次分析法解答
a=[13.54,14.36,87.46,566.3,0.09779
13.08,15.71,85.63,520,0.1075
9.504,12.44,60.34,273.9,0.1024
17.99,10.38,122.8,1001,0.1184
20.57,17.77,132.9,1326,0.08474
19.69,21.25,130,1203,0.1096
11.42,20.38,77.58,386.1,0.1425
20.29,14.34,135.1,1297,0.1003];
x=[16.6,28.08,108.3,858.1,0.08455
20.6,29.33,140.1,1265,0.1178
7.76,24.54,47.92,181,0.05263];
g=[ones(3,1);2*ones(5,1)];
[class,err]=classify(x,a,g)结果：
>> e09
class =
1
1
1
err =
0
795
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
