%a program for 單向分組資料的協方差分析，k組資料，每組有n對觀察值（正確）
clear;clc;
x0=xlsread('danxiangxiefangchax');%P366
y0=xlsread('danxiangxiefangchay');

% x0=xlsread('X3');
% y0=xlsread('Y31');
% y0=xlsread('Y32');
% y0=xlsread('Y33');

% x0=xlsread('X413');
% y0=xlsread('Y413');

x=x0';
y=y0';
k=size(x,2);
n=size(x,1);
x1=x(:);
y1=y(:);
z=[x1,y1];
m=size(z,1);
%SSTx的自由度dfTx
dfT=k*n-1;
%SStx的自由度dftx
dft=k-1;
spT=cov(z)*dfT;
ti=zeros(k,2);
z1=zeros(n,2,k);
for i=1:k
    z1(:,:,i)=z(n*(i-1)+1:n*i,:);
end
for s=1:k
    for i=1:n
        for j=1:2
            ti(s,j)=ti(s,j)+z1(i,j,s);
        end
    end
end
spt=dft*cov(ti)/n;
spe=spT-spt;
%各樣本內（處理內）平方和SSeix(i)
zi=ti/n;
SSei=zeros(k,2);
for s=1:k
    for i=1:n
        for j=1:2
            SSei(s,j)=SSei(s,j)+(z1(i,j,s)-zi(s,j))^2;
        end
    end
    dfei(s)=n-1;%分別具有自由度dfei(i)
end
%SSex的自由度dfex
dfe=k*(n-1);
MSt=spt/dft;
MSe=spe/dfe;
%%%%cut
%將SPe分解為SPei(i),i=1:k
SPei1=zeros(k,1);SPei2=zeros(k,1);
for i=1:k
    for j=1:n
        SPei1(i)=SPei1(i)+x0(i,j)*y0(i,j);
    end
    SPei2(i)=SPei2(i)+ti(i,1)*ti(i,2)/n;
    SPei(i,1)=SPei1(i)-SPei2(i);
end
%%%%%%%%%
Q=0;
for i=1:k
    %     bei(i)=SPei(i)/SSei(i,1)
    %     ai(i)=zi(i,2)-b(i)*zi(i,1)
    Qi(i)=SSei(i,2)-SPei(i)^2/SSei(i,1);
    Q=Q+Qi(i);
end
dfQ=k*(n-2); %P371
% %%%%%%%
be=spe(1,2)/spe(1,1);
for i=1:k
    ae(i)=zi(i,2)-be*zi(i,1);
end
Qe=spe(2,2)-spe(1,2)^2/spe(1,1);
dfQe=dfe-1;
% %%%%%%%%
% bt=SPt/SStx;
% at=ybar-bt*xbar;
% Qt=SSty-SPt^2/SStx;
% dfQt=dft-1;
% %%%%%%%
% bT=SPT/SSTx;
% aT=ybar-bT*xbar;
QT=spT(2,2)-spT(1,2)^2/spT(1,1);
dfQT=dfT-1;
% %以上為八個線性方程
%%%%%%%%
sym('[各品種的y依x用共同的回歸系數be的方程如下：]')
traits=['a';'b'];
fangchenge=[];
for i=1:k
    fangchenge=[fangchenge;ae(i),be];
end
len=size(fangchenge,2);
fangchenge=vpa(fangchenge,5);
line1=vpa(ones(1,len),2);
fangchenge=[line1;fangchenge];
fangchenge(1,1)=traits(1,1);
fangchenge(1,2)=traits(2,1);
fangchenge
zb=zeros(1,2);
for i=1:k
    for j=1:2
        zb(j)=zb(j)+zi(i,j);
    end
end
zb=zb/k;
z2=z1;
for s=1:k
    for i=1:n
        z2(i,2,s)=z1(i,2,s)-be*(z1(i,1,s)-zb(1,1));%第1列是x值，第2列是y的矯正值
    end
end
%處理平均數yibar的回歸矯正
for i=1:k
    yibarjiaozheng(i)=zi(i,2)-be*(zi(i,1)-zb(1,1));
end
y2=yibarjiaozheng; %平均數yibar的矯正值
%對觀察值y進行矯正
ti2=zeros(k,2);
for s=1:k
    for i=1:n
        for j=1:2
            ti2(s,j)=ti2(s,j)+z2(i,j,s);
        end
    end
end
for j=1:n
for i=1:k
       z3(((1+(i-1)*n):i*n),:)=z2(:,:,i);%矯正的觀察值，對應于z(m,2)
   end
end
spT2=cov(z3)*dfT;
spt2=dft*cov(ti2)/n;
spe2=spT2-spt2;
% %各樣本內（處理內）平方和SSeix(i)
% zi=ti2/n;
% SSei2=zeros(k,2);
% for s=1:k
%     for i=1:n
%         for j=1:2
%             SSei2(s,j)=SSei2(s,j)+(z1(i,j,s)-zi(s,j))^2;
%         end
%     end
% end
be2=spe2(1,2)/spe2(1,1);
% for i=1:k
%     ae2(i)=zi(i,2)-be2*zi(i,1);
% end
% traits=['a';'b'];
% fangchenge2=[];
% for i=1:k
%     fangchenge2=[fangchenge2;ae2(i),be2];
% end
% len=size(fangchenge2,2);
% fangchenge2=vpa(fangchenge2,5);
% line1=vpa(ones(1,len),2);
% fangchenge2=[line1;fangchenge2];
% fangchenge2(1,1)=traits(1,1);
% fangchenge2(1,2)=traits(2,1);
% fangchenge2
%%%%%%%%%
MSt=spt/dft;
MSe=spe/dfe;
f2=MSt./MSe;
a=1-fcdf(f2,dft,dfe);
if a(1,2)>=0.05
    ans=sym('[y與x互作效應不顯著]')
else 
    ans=sym('[y與x互作效應顯著]')
end
%%%%%%%%%%%%%%
% ask=sym('[是否k個處理的總體平均數yibar相等？]')
if a(2,2)>=0.05
    ans=sym('[k個處理的總體平均數yibar沒有顯著差異]')
    ae2=zb(1,2)-be2*zb(1,1);
elseif a(2,2)<0.05
    ans=sym('[k個處理的總體平均數yibar有顯著差異：]')
    sy12bar=sqrt(2*MSe/n);
    for j=1:2
        ci=0;
        for s1=1:k
            for s2=1:k
                if s1<s2
                    ci=ci+1;
                    t1(ci,j)=(zi(s1,j)-zi(s2,j))/sy12bar(j,j);
                    ii(ci,j)=s1;jj(ci,j)=s2;
                end
            end
        end
    end
    b=2*(1-tcdf(abs(t1),dfe));
    c1=0;
    for i=1:ci
        for j=1:2
            if b(i,j)>=0.05
                c1=c1+1;xz(i,j)=0;
            elseif b(i,j)>=0.01 
                c1=c1+1;xz(i,j)=1;
            elseif b(i,j)<0.01 
                c1=c1+1;xz(i,j)=2;
            end  
        end
    end
    isx=[];
%     traits=['xih';'xjl';'xxz';'yih';'yjl';'yxz'];
    traits=['yih';'yjl';'yxz';'y12'];
    len=size(traits,1);
    line1=vpa(ones(1,len),6);
    isx=[line1;isx];
    for i=1:len
        isx(1,i)=traits(i,:);
    end
    for i=1:ci
        y12(i)=y2(ii(i,2))-y2(jj(i,2));
    end
    y12=vpa(y12,6);
    for i=1:ci
%         isx(i+1,:)=[ii(i,1),jj(i,1),xz(i,1),ii(i,2),jj(i,2),xz(i,2)]; 
isx(i+1,:)=[ii(i,2),jj(i,2),xz(i,2),y12(i)];
    end
    isx
end
%%%%%%%%%%%%%%%%%
sym('[方差協方差分析：]')
sa(1)=sym('[t]');sa(2)=sym('[e]');sa(3)=sym('[T]');
df(1)=dft;df(2)=dfe;df(3)=dfT;
ss(1)=spt(2,2);ss(2)=spe(2,2);ss(3)=spT(2,2);
ms(1)=ss(1)/df(1);ms(2)=ss(2)/df(2);ms(3)=0;
f(1)=f2(2,2);f(2)=0;f(3)=0;
sp(1)=spt(1,2);sp(2)=spe(1,2);sp(3)=spT(1,2);
mp(1)=sp(1)/df(1);mp(2)=sp(2)/df(2);mp(3)=0;
fp(1)=f2(1,2);fp(2)=0;fp(3)=0;
fangcha=[];
for i=1:3
    fangcha=[fangcha;sa(i),df(i),ss(i),ms(i),f(i),sp(i),mp(i),fp(i)];
end
len=size(fangcha,2);
fangcha=vpa(fangcha,6);
traits=['SOV';'DF ';'SS ';'MS ';'F  ';'SP ';'MP ';'Fp '];
line1=vpa(ones(1,len),6);
fangcha=[line1;fangcha];
for i=1:len
fangcha(1,i)=traits(i,:);
end
ANOVA=fangcha