プロシージャ

> restart;

> anova2r:=proc(data,r)
> local i, j, k, a, b, alld, alldr, N, Xbar, ST, SA, SB, SAB, SAxB, Am, Bm, ABm, SE, pT, pA, pB, pAB, pAxB, pE, VA, VB, VAxB, VE, F0A, F0B, F0AxB, L1, L2, L3, L4, L5, L6, L, flga, flgb, flgaxb, pva, pvb, pvaxb, r1, r5, rb1, rb5, rc1, rc5, g1, g2:
> with(plots):with(stats):with(describe):
> a:=nops(data):b:=nops(data[1])/r:
> alld:=[seq(op(data[i]),i=1..a)]:
> N:=nops(alld):
> Xbar:=evalf(mean(alld),5):
> ST:=add((alld[i]-Xbar)^2,i=1..N):
> ABm:=evalf([seq([seq(add(data[k][r*(i-1)+j],j=1..r)/r,i=1..b)],k=1..a)],5):
> alldr:=[seq(op(ABm[i]),i=1..a)]:
> Am:=evalf([seq(mean(ABm[i]),i=1..a)],5):
> SA:=add((Am[i]-Xbar)^2,i=1..a)*b*r:
> Bm:=evalf([seq(mean([seq(ABm[i,j],i=1..a)]),j=1..b)],5):
> SB:=add((Bm[i]-Xbar)^2,i=1..b)*a*r:
> SAB:=add((alldr[i]-Xbar)^2,i=1..a*b)*r:
> SE:=ST-SAB:SAxB:=SAB-SA-SB:
> pT:=N-1: pA:=a-1: pB:=b-1: pAB:=a*b-1: pE:=pT-pAB: pAxB:=pAB-pA-pB:
> VA:=SA/pA: VB:=SB/pB: VAxB:=SAxB/pAxB: VE:=SE/pE:
> F0A:=VA/VE: F0B:=VB/VE: F0AxB:=VAxB/VE:
> statevalf[cdf,fratio[pA,pE]](F0A):pva:=evalf(1-%,5):
> statevalf[cdf,fratio[pB,pE]](F0B):pvb:=evalf(1-%,5):
> statevalf[cdf,fratio[pAxB,pE]](F0AxB):pvaxb:=evalf(1-%,5):
> if pva<0.01 then flga:=`(**)` elif pva<0.05 then flga:=`( *)` else flga:=`( )` fi:
> if pvb<0.01 then flgb:=`(**)` elif pvb<0.05 then flgb:=`( *)` else flgb:=`( )` fi:
> if pvaxb<0.01 then flgaxb:=`(**)` elif pvaxb<0.05 then flgaxb:=`( *)` else flgaxb:=`( )` fi:
> L1:=[` `,phi,S,V,`F0`,`1-P`,` `]:
> L2:=evalf([`A`,pA,SA,VA,F0A,1-pva,flga],5):
> L3:=evalf([`B`,pB,SB,VB,F0B,1-pvb,flgb],5):
> L4:=evalf([`AxB`,pAxB,SAxB,VAxB,F0AxB,1-pvaxb,flgaxb],5):
> L5:=evalf([`E`,pE,SE,VE,`--`,`--`,` `],5):
> L6:=evalf([`T`,pT,ST,`--`,`--`,`--`,` `],5):
> print(`**** ANOVA Table (two-way, repeated data) ****`);
> print(`N=`,N,`a=`,a,`b=`,b,`r=`,r,`X_bar=`,Xbar);
> L:=[L1,L2,L3,L4,L5,L6]:print(convert(L,array));
> # Pooling
> if pvaxb>0.2 then
> pE:=pE+pAxB:SE:=SE+SAxB:VE:=SE/pE:
> F0A:=VA/VE:F0B:=VB/VE:
> statevalf[cdf,fratio[pA,pE]](F0A):pva:=evalf(1-%,5):
> statevalf[cdf,fratio[pB,pE]](F0B):pvb:=evalf(1-%,5):
> if pva<0.01 then flga:=`(**)` elif pva<0.05 then flga:=`( *)` else flga:=`( )` fi:
> if pvb<0.01 then flgb:=`(**)` elif pvb<0.05 then flgb:=`( *)` else flgb:=`( )` fi:
> L1:=[` `,phi,S,V,`F0`,`1-P`,` `]:
> L2:=evalf([`A`,pA,SA,VA,F0A,1-pva,flga],5):
> L3:=evalf([`B`,pB,SB,VB,F0B,1-pvb,flgb],5):
> L5:=evalf([`E`,pE,SE,VE,`--`,`--`,` `],5):
> L6:=evalf([`T`,pT,ST,`--`,`--`,`--`,` `],5):
> print(` `);
> print(`**** ANOVA Table after pooling (two-way, repeated data) ****`);
> L:=[L1,L2,L3,L5,L6]:print(convert(L,array));
> with(plots):
> g1:=plot({seq([seq([k,ABm[j,k]],k=1..b)],j=1..a)},color=black,labels=[`B`,` `]):
> g2:=plot({seq([seq([k,ABm[j,k]],k=1..b)],j=1..a)},color=black,style=point):
> print(display(g1,g2));
> end: