分散分析（直交表）のRのプログラム その３【2/12更新】

プログラムの続き

# Taguchi SN ratio
# References
#
# Taguchi, Gen-ichi et. al., コンピュータによる情報設計の技術開発,
# "Engeenering R&D of Information Design with Computer Simulation", 2004,
# ISBN:4-452-51115-4
#
# Taguchi, Gen-ichi et. al., 開発・設計段階の品質工学,
# "Quality Engeenering in R&D and Design", 1988,
# ISBN:4-542-51101-4
# 
# Taguchi, Gen-ichi, 品質工学の数理,
# "Mathematics for Quality Engineering", 1999,
# ISBN:4-542-511118-9
#
# Miyagawa, Masami, 品質を獲得する技術　タグチメソッドがもたらしたもの,
# "Technology for Getting Quality What the Taguchi Method Has Brought Us", 2000,
# ISBN:4-8171-0339-6
#
# Eppinger, Stiven, Robust Design: Experiments for Better Products, 2002? ,
# MIT Open Courseware, http://ocw.mit.edu/NR/rdonlyres/Sloan-School-of-Management/15-783JProduct-Design-and-DevelopmentSpring2002/20CA897E-5419-4EC2-B54A-31460AB444F9/0/robust_design.pdf
#
# Funao, Nobuo, The R Tips, 2005, ISBN:4-86167-039-X
#
# RjpWiki,P Johnson tips1,http://www.okada.jp.org/RWiki/index.php?P_Johnson_tips_1
# RjpWiki,P Johnson tips2,http://www.okada.jp.org/RWiki/index.php?P_Johnson_tips_2
# 
# RjpWiki,http://www.okada.jp.org/RWiki/index.php?RjpWiki
#
# R Development Core Team,  R: A language and environment for statistical computing (Ver. 2.21), 2005,
# ISBN 3-900051-07-0, URL http://www.R-project.org.
#

latin.factor.calc<-function(ltable,colno,fun="mean"){

 	fname<-colnames(ltable)

	ycolno<-length(fname)

	nfac<- length(levels(ltable[ ,colno]))

 	rs<-data.frame(fname=rep(fname[colno],nfac),flvls=levels(ltable[ ,colno]),y=as.vector(tapply(ltable[ ,ycolno],ltable[ ,colno],fun)))

 	return(rs)
}

latin.factor.plotone<- function(facdata,ylim=NULL,...){

	fnames<-colnames(facdata)

	matplot(x=facdata[,2],y=facdata[,3],type="b",
		pch=1,xaxt="n",xlab=as.character(facdata[1,1]),ylab=fnames[3],ylim=ylim,...)
 	axis(side=1,at=facdata$flvl)
}

latin.factor.plot<- function(ltable,ncol,fun="mean",mfrow=c(1,1),ylim=NULL,title=NULL,...){

	tmppar<-par(no.readonly=TRUE)

	par(oma=c(0,0,4,0))
	par(mfrow=mfrow)
	dfadd<-NULL

	for(i in 1:ncol){
 		rs<-latin.factor.calc(ltable=ltable,colno=i,fun=fun)

		latin.factor.plotone(rs,ylim=ylim,...)
		dfadd<-rbind(dfadd,rs)

	}
	mtext(side=3,line=1,outer=TRUE,text=title,cex=1.5)
 	par(tmppar)
	return(dfadd)
}

taguchi.param.sn.stdone<- function(N0,noise){


 	nnoise<- ncol(noise)

	k<- nrow(N0)

	L<- rep(0,nnoise)

	for(i in 1:nnoise){
		L[i]<- sum(N0[,1]*noise[,i])

	}
	
	fSb<-  1

	fSnb<- (nnoise-1)

	fSe<-  (nnoise*k-fSb-fSnb)

	fSnbSe<- (fSe+fSnb)

	fSt<-  (nnoise*k)	

	
	St<- sum(noise^2)

	r<- sum(N0^2)

	Sb<- (sum(L)^2)/(nnoise*r)

	Snb<- (sum(L^2)/r)-Sb

	Se<- St -(Sb+Snb)

	SnbSe<- Snb+Se

	Vb<- Sb/fSb

	Vnb<- Snb/fSnb

	Ve<- Se/fSe

	Vn<- (SnbSe)/fSnbSe

	SN<- 10*log10((nnoise*r)/Vn)

	
	Df<- c(fSb,fSnb,fSe,fSnbSe,fSt)

	Sdata<- c(Sb,Snb,Se,SnbSe,St)

	Vdata<- c(Vb,Vnb,Ve,Vn,NA)

	
	atable<- data.frame(Df=Df,SumSq=round(Sdata,digits=4),MeansSq=round(Vdata,digits=4))

	rownames(atable)<- c("b","N*b","e","N+e","T")

	
	res<- list(SN=SN,ANOVA=atable)


	return(res)
}


taguchi.param.sn.std<- function(noise){


 	nfno<- length(levels(noise[,1]))

	nfsig<- length(levels(noise[,2]))

	nfnoise<- length(levels(noise[,3]))

	SN<- rep(0,nfno)

	ANOVA<- list(rep(NA,nfno))


	for (i in 1:nfno){
		Ndata<- matrix(0,ncol=nfnoise-1,nrow=nfsig)

		Nn<- subset(noise,fno==i)

		tmp<- subset(Nn,fnoise==0,noisedata)

		N0<- as.matrix(tmp[[1]])

		for(j in 1:(nfnoise-1)){
			tmp<- as.vector(subset(Nn,fnoise==j,noisedata))

			Ndata[,j]<- tmp[[1]]

		}
		res<- taguchi.param.sn.stdone(N0=N0,noise=Ndata)

		
			SN[i]<- res$SN

			ANOVA[i]<- list(res$ANOVA)

	}

	return(list(SN=SN,ANOVA=ANOVA))

}


taguchi.param.tune<- function(y,m){


	nm<- length(m)

	ny<- length(y)

	if(nm != ny){
		stop("yとmの行数が違います．")
	}
	
	w<- rep(0,nm)

	
	fSt<- nm

	fSb1<- 1

	fSb2<- 1

	fSe<- nm-2

	
	St<- sum(y^2)

	
	L1<- sum(m*y)

	r1<- sum(m^2)

	kr3<- sum(m^3)

	Sb1<- (L1^2)/r1

	b1<- L1/r1

	
	K2<- r1/nm

	K3<- kr3/nm

	K<- K3/K2

	
	for(i in 1:nm){
		w[i]<- (m[i]^2)-K*m[i]

	}
	
	L2<- sum(w*y)

	r2<- sum(w^2)

	Sb2<- (L2^2)/r2

	b2<- L2/r2

	
	Se<- St-Sb1-Sb2

	Ve<- Se/fSe

	
	
	estimate<- c(b1,b2,NA,NA)

	Df<- c(fSb1,fSb2,fSe,fSt)	

	Sumsq<- c(Sb1,Sb2,Se,St)

	meansq<- c(NA,NA,Ve,NA)


	result<- data.frame(

			    Estimate=round(estimate,digits=4),
			    Df=Df,
			    SumSq=round(Sumsq,digits=4),
			    MeansSq=round(meansq,digits=4))
 
 	rownames(result)<- c("b1","b2","e","T")

 
	return(result)
}

taguchi.param.sn.smaller<- function(y){


	nrowy<- nrow(y)

	ncoly<- ncol(y)

	SN<- rep(0,nrowy)

	
	for (i in 1:nrowy){
		St<- sum(y[i,]^2)

		Ve<- St / ncoly

		SN[i]<- -10*log10(Ve)

	}
	return(SN)
}

taguchi.param.sn.larger<- function(y){


	nrowy<- nrow(y)

	ncoly<- ncol(y)

	SN<- rep(0,nrowy)

	
	for (i in 1:nrowy){
		SN[i]<- -10*log10(sum(1/(y[i,]^2))/ncoly)

	}
	return(SN)
}

taguchi.param.sn.nominal<- function(y){


	nrowy<- nrow(y)

	ncoly<- ncol(y)

	SN<- matrix(0,ncol=2,nrow=nrowy)

	colnames(SN)<- c("SN","Sig")

	for(i in 1:nrowy){
		St<- sum(y[i,]^2)

		Sm<- (sum(y[i,])^2)/ncoly

		Ve<- (St-Sm)/(ncoly-1)

		Sig<- (Sm-Ve)/ncoly

		SN[i,1]<- 10*log10((Sig)/Ve)

		SN[i,2]<- 10*log10(Sig)

	}
	return(SN)
	
}

taguchi.param.sn.zero<- function(y){


	nrowy<- nrow(y)

	ncoly<- ncol(y)

	SN<- rep(0,nrowy)

	for (i in 1:nrowy){
		St<- sum(y[i,]^2)

		Ve<- St / ncoly

		SN[i]<- 10*log10(1/Ve)

	}
	return(SN)

}


#
#Wheatstone Bridge Program
#
# References
#
# Taguchi, Gen-ichi et. al., 開発・設計段階の品質工学,
# "Quality Engeenering in R&D and Design", 1988,
# ISBN:4-542-51101-4
# 
# Miyagawa, Masami, 品質を獲得する技術　タグチメソッドがもたらしたもの,
# "Technology for Getting Quality What the Taguchi Method Has Brought Us", 2000,
# ISBN:4-8171-0339-6
#
# Funao, Nobuo, The R Tips, 2005, ISBN:4-86167-039-X
#
# RjpWiki,http://www.okada.jp.org/RWiki/index.php?RjpWiki
#
# R Development Core Team,  R: A language and environment for statistical computing (Ver. 2.21), 2005,
# ISBN 3-900051-07-0, URL http://www.R-project.org.
#

wheatstoneBridge.one<-function(fac){


	f1<-(fac["RB"]*fac["RD"])/fac["RC"]

	f2<-fac["XX"]/(fac["EE"]*fac["RC"]^2)

	f3<-fac["RA"]*(fac["RD"]+fac["RC"])+fac["RD"]*(fac["RB"]+fac["RC"])

	f4<-fac["RB"]*(fac["RC"]+fac["RD"])+fac["RF"]*(fac["RB"]+fac["RC"])

	fac["RRres"]<- f1-f2*f3*f4

	
	return(fac)
	
}

wheatstoneBridge.RB<- function(lvl.result,RR){


	ncol<- ncol(lvl.result)

	for (i in 1:ncol){
		lvl.result["RB",i]<- (RR*lvl.result["RC",i])/lvl.result["RD",i]

	}
	
	return(lvl.result)

}

lvl.make.L18<- function(lvl){

	
	y<- rep(0,18)

	ltable<- latin.table.L18(y)

	ltable<- as.matrix(ltable)

	ltable<- as.numeric(ltable)

	ltable<- matrix(ltable,ncol=9)


	lvl.result<- matrix(rep(0,144),nrow=8,ncol=18)


	for (i in 1:8){
		lvl.result[i,]<- lvl[i,ltable[1:18,i]]

	}
	
	return(lvl.result)

}

wheatstoneBridge.L18<- function(ilvl.result,olvl.result){


	icol<- ncol(ilvl.result)

	ocol<- ncol(olvl.result)


	y.result<- matrix(rep(0,icol*ocol),ncol=ocol,nrow=icol)


	for (i in 1:icol){
		for (j in 1:ocol){

			fac<-ilvl.result[,i]*olvl.result[,j]

			fac<- c(fac,0)

			names(fac)<-c("nF","RA","RC","RD","EE","RF","RB","XX","RRres")

			fac.result<- wheatstoneBridge.one(fac)

			y.result[i,j]<- fac.result["RRres"]

		}
	}

	return(y.result)

}