zjfc.m

MATLAB科学计算与工程分析源代码源程序3

function [result,error,errorDim]=zjfc(s,opt)
% 对正交试验进行方差分析，s是输入矩阵，opt是空列参数向量，给出s中是空白列的列序号
%s=[  1  1   1   1  1 1 1 83.4;
%     1  1   1   2  2 2 2 84;
%     1  2   2   1  1 2 2 87.3;
%     1  2   2   2  2 1 1 84.8;
%     2  1   2   1  2 1 2 87.3;
%     2  1   2   2  1 2 1 88;
%     2  2   1   1  2 2 1 92.3;
%     2  2   1   2  1 1 2 90.4;
%];
%opt=[3,7];
% s的最后一列是各个正交组合的试验测量值，前几列是正交表
[m,n]=size(s);
 p=max(s(:,1));%取水平数
 q=n-1;% 取列数
 sum0=zeros(p,q);
 for i=1:q
    for k=1:m 
         for j=1:p
          if(s(k,i)==j)
              sum0(j,i)=sum0(j,i)+s(k,n);%求和
          end
        end
     end
 end
totalsum=sum(s(:,n));
ss=sum0.*sum0;
levelsum=m/p; %水平重复数
ss=sum(ss./levelsum)-totalsum^2/m; %每一列的S
ssError=sum(ss(opt));
for i=1:q
    f(i)=p-1;    %自由度
end
fError=sum(f(opt));  %误差自由度
ssbar=ss./f;
Errorbar=ssError/fError;
index=find(ssbar<Errorbar);
index1=find(index==opt);
index(index==index(index1))=[];%剔除重复
ssErrorNew=ssError+sum(ss(index));  %并入误差
fErrorNew=fError+sum(f(index));  %新误差自由度
F=(ss./f)/(ssErrorNew./fErrorNew);   %F值
errorDim=[opt,index];

errorDim=sort(errorDim);%误差列的序号
result=[ss',f',ssbar',F'];
error=[ssError,fError;ssErrorNew,fErrorNew];