Rでタグチメソッド - how to code something

richbutti.blogspot.com

Latin_Table.csv

A,B,C,D,F,G,H,I,J
L1,L1,L1,L1,L1,L1,L1,L1,L1
L1,L1,L1,L1,L1,L2,L2,L2,L2
L1,L1,L2,L2,L2,L1,L1,L1,L2
L1,L2,L1,L2,L2,L1,L2,L2,L1
L1,L2,L2,L1,L2,L2,L1,L2,L1
L1,L2,L2,L2,L1,L2,L2,L1,L2
L2,L1,L2,L2,L1,L1,L2,L2,L1
L2,L1,L2,L1,L2,L2,L2,L1,L1
L2,L1,L1,L2,L2,L2,L1,L2,L2
L2,L2,L2,L1,L1,L1,L1,L2,L2
L2,L2,L1,L2,L1,L2,L1,L1,L1
L2,L2,L1,L1,L2,L1,L2,L1,L2

Output.csv

y1,y2,y3,y4,y5,y6,y7,y8,y9
15.9,14.5,13.1,14.4,13,13.8,12.9,13.7,12.3
15.7,13.5,11.4,13.4,11.2,13.6,11.1,13.5,11.3
18.3,16.8,15.3,16,14.5,17.5,13.8,16.8,15.3
17.7,16.3,14.9,16.9,15.5,16.4,16.2,17,15.6
16.3,15.6,14.8,15.6,14.8,16.3,14.8,16.3,15.6
18.1,16.6,15.1,16.6,15.1,18.1,15.1,18.1,16.6
18.2,17.4,16.7,16.7,15.9,17.4,15.2,16.7,15.9
16.2,15.4,14.7,14.7,13.9,15.4,13.2,14.7,13.9
17.3,15.1,13,15,12.9,15.2,12.8,15.1,13
18.7,17.2,15.7,17.2,15.7,18.7,15.7,18.7,17.2
17.3,15.9,14.5,16.5,15.1,16,15.8,16.6,15.2
18,15.9,13.7,16.5,14.4,16.7,15,17.4,15.2

L12<-read.table("Latin_Table.csv",header = TRUE,sep = ",")
Output<-read.table("Output.csv",header = TRUE,sep = ",")

Mean <- apply(Output,1,mean)
Mean_Square <- Mean*Mean
Vi <- Mean - Output
Vi <- Vi*Vi
Vi <- apply(Vi,1,sum)
Vi <- Vi/(ncol(Output)-1)
SN <- (Mean_Square-(Vi/ncol(Output)))/Vi
SN <- 10*log10(SN)
SN_Bunsan <- cbind(L12,SN)
SN.lm <- lm(SN~A+B+C+D+F+G+H+I+J,data = SN_Bunsan)
summary(SN.lm)
anova(SN.lm)
tmppar<-par(no.readonly=TRUE)
par(mfrow=c(1,ncol(L12)))
SN_max <- max(SN)
SN_min <- min(SN)
fnames<-colnames(SN_Bunsan)
x_fc <- c(NULL)
y_SN <- c(NULL)
for(i in 1:ncol(L12)){
  fc <- factor(SN_Bunsan[,i])
  x_tmp <- 1:length(levels(fc))
  y_tmp <- tapply(SN_Bunsan$SN,fc,mean)
  plot(x_tmp,y_tmp,type="b",pch=1,xaxp=c(1,length(x_tmp),1),
       ylim=c(SN_min,SN_max),xlab=as.character(fnames[i]),
       ylab="SN_Ratio[dB]",col="red")
}


Sensitivity <- Mean_Square-(Vi/ncol(Output))
Sensitivity <- 10*log10(Sensitivity)
Sens_Bunsan <- cbind(L12,Sensitivity)
Sens_Bunsan
Sens.lm <- lm(Sensitivity~A+B+C+D+F+G+H+I+J,data = Sens_Bunsan)
summary(Sens.lm)
anova(Sens.lm)
par(mfrow=c(1,ncol(L12)))
sens_max <- max(Sensitivity)
sens_min <- min(Sensitivity)
fnames<-colnames(Sens_Bunsan)
for(i in 1:ncol(L12)){
  fc <- factor(Sens_Bunsan[,i])
  x_tmp <- 1:length(levels(fc))
  y_tmp <- tapply(Sens_Bunsan$Sensitivity,fc,mean)
  plot(x_tmp,y_tmp,type="b",pch=1,xaxp=c(1,length(x_tmp),1),
       ylim=c(sens_min,sens_max),xlab=as.character(fnames[i]),
       ylab="Sensitivity[dB]",col="red")
}