Rで実験計画法L8の確認

L8の実験関数
v1:第一水準
v2:第二水準
f:関数

L8<-function(x=1.0,fun=f,v1,v2){
	
	mat<-cbind(v1,v2)
	e1<-c(1,1,1,1,1,1,1)
	e2<-c(1,1,1,2,2,2,2)
	e3<-c(1,2,2,1,1,2,2)
	e4<-c(1,2,2,2,2,1,1)
	e5<-c(2,1,2,1,2,1,2)
	e6<-c(2,1,2,2,1,2,1)
	e7<-c(2,2,1,1,2,2,1)
	e8<-c(2,2,1,2,1,1,2)
	e<-cbind(e1,e2,e3,e4,e5,e6,e7,e8)

	out<-NULL
	for(i in 1:length(e[1,])){
		q1<-which(e[,i]==1)
		q2<-which(e[,i]==2)
		v<-NULL
		v[q1]<-mat[q1,1]
		v[q2]<-mat[q2,2]

		out[i]<-f(v)
	}
	out

f<-function(v){
#	5*v[1]+4*v[2]+10*v[4]+rnorm(1,mean=0,sd=1);

	5*v[1]+4*v[2]+1*5*sin(v[4])+rnorm(1,mean=0,sd=2);
}


L8<-function(x=1.0,fun=f,v1,v2){

	out<-NULL
	out[1]<-f(v1)
	v=c(v1[1],v1[2],v1[3],v2[4],v2[5],v2[6],v2[7])	
	out[2]<-f(v)

	v=c(v1[1],v2[2],v2[3],v1[4],v1[5],v2[6],v2[7])	
	out[3]<-f(v)

	v=c(v1[1],v2[2],v2[3],v2[4],v2[5],v1[6],v1[7])	
	out[4]<-f(v)


	v=c(v2[1],v1[2],v2[3],v1[4],v2[5],v1[6],v2[7])	
	out[5]<-f(v)


	v=c(v2[1],v1[2],v2[3],v2[4],v1[5],v2[6],v2[7])	
	out[6]<-f(v)


	v=c(v2[1],v2[2],v1[3],v1[4],v2[5],v2[6],v1[7])	
	out[7]<-f(v)


	v=c(v2[1],v2[2],v1[3],v2[4],v1[5],v1[6],v2[7])	
	out[8]<-f(v)
	out
}



v1=rep(1,8)
v2=rep(3,8)
q<-L8(fun=f,v1=v1,v2=v2)
q

a<-c(1,1,1,1,2,2,2,2)
b<-c(1,1,2,2,1,1,2,2)
c<-c(1,2,1,2,1,2,1,2)

sa<-mean(sum((q[which(a==1)]-mean(q)))^2+(sum(q[which(a==2)]-mean(q)))^2)/4
sb<-mean(sum((q[which(b==1)]-mean(q)))^2+(sum(q[which(b==2)]-mean(q)))^2)/4
sc<-mean(sum((q[which(c==1)]-mean(q)))^2+(sum(q[which(c==2)]-mean(q)))^2)/4

#cf<-sum(q)^2/length(q)
#sa1<-mean(sum(q[which(a==1)])^2+(sum(q[which(a==2)]))^2)/4-cf
#sb1<-mean(sum(q[which(b==1)])^2+(sum(q[which(b==2)]))^2)/4-cf
#sc1<-mean(sum(q[which(c==1)])^2+(sum(q[which(c==2)]))^2)/4-cf

st=sum((q-mean(q))^2)
se=st-sa-sb-sc
dfa<-1
dfb<-1
dfc<-1
dcf<-1
dfe=8-dfa-dfb-dfc-dcf

pfa<-1-pf(sa/dfa/se*dfe,dfa,dfe)
pfb<-1-pf(sb/dfb/se*dfe,dfa,dfe)
pfc<-1-pf(sc/dfc/se*dfe,dfa,dfe)

#a[which(a==1)]=v1[1];
#a[which(a==2)]=v2[1];

#b[which(a==1)]=v1[2];
#b[which(a==2)]=v2[2];

#c[which(a==1)]=v1[4];
#c[which(a==2)]=v2[4];


d<-data.frame(A=a,B=b,C=c,y=q)
ans<-aov(d$y~d$A+d$B+d$C)
ans2<-anova(ans)
print(ans2)