Rで計算する、タグチメゾット（品質工学）を使ったパラメータ設計 に関するメモ② （Design of experiments & Taguchi method）

前回は、静特性のパラメータ設計についてまとめました。

■動特性のパラメータ設計について

今回は、動特性のパラメータ設計について学びます。

入力の水準に応じて、結果が変わる特性を動特性と呼びます。

動特性を行う際の直交表は、内側配置に制御因子を、外側配置に信号因子を用いることが多いようです。また、誤差因子と信号因子を組み合わせた場合もあるようです。

動特性の場合も、静特性の場合の様に、いくつかに大きく区分できます。
ゼロ点比例式、基準点比例式、一次式、などに別れるようです。
以下では、一次式の事例をあげていますが、本質的にどの様に誤差計算を実施しているかが理解できれば、様々なモデルに応用できると思います。

では、動特性の場合の感度とSN比について確認していきたいと思います。

外側配置の信号因子M1,・・・,Mnがで構成されている場合に、その解きの第i番目の制御因子の組み合わせに対する出力をyi1,yi2,・・・,yinとします。

その際に、第i番目の制御因子の水準組み合わせにおいて、応答yと因子Mの間における、以下の関係について考えます。

上記の式から、傾きであるβ1iが大きいほどその出力は敏感であると考えられます。一方で、データは上記の直線の周辺にバラつきます。そのバラつきをあらわすσ_i^2については、小さい方が出力が敏感であることになります。

これらの概念を組み合わせすると、出力の鋭敏さの測度として、以下を用いるとよいことが考えられます。

上記の推定について考えます。

σ_i^2は、第i水準組み合わせでの残差分散V_eiで推定ができます。

では、β_1i^2を推定について考えます。

こちらは、最小二乗法で求めて終わりかと思いきや、βを二乗したものとの不偏性の観点から、以下で推定します。

以上から、鋭敏さの測度の考え方とあわせて、下記でSN比を定義します。

また、感度としては、以下を用いることができます。

S_riは下記で求めることができます。

上記の流れの式は、線形式や入力の二乗和のrなどを用いる式の表現とは異なりますが、計算結果は一致しますのでご安心を。（※モデルが一次モデルなのか、比例モデルなのかは注意してください）

再びですが、計算方法よりも、品質工学の場合は、式が意味するものをしっかり把握しておくのが、計算ミスをなくすのにベターだと思います。

■実際のデータで計算

では、今回も実際のデータで計算したいと思います。

以下のような直交表を考えます。

										1	2	3	4	5	6	7	8	9	10
									M	20	20	40	40	60	60	80	80	100	100
	A	B	C	D	F	G	H	I	J
1	L1	L1	L1	L1	L1	L1	L1	L1	L1	19.9	20	39.9	40.1	60.1	59.9	80	80.2	100	100.5
2	L1	L1	L1	L1	L1	L2	L2	L2	L2	11.9	11.8	24	23.7	35.6	36.1	47.7	47.9	59.9	60
3	L1	L1	L2	L2	L2	L1	L1	L1	L2	27.9	28	55.6	55.9	84.1	83.9	112.5	111.9	140.6	140.7
4	L1	L2	L1	L2	L2	L1	L2	L2	L1	3.9	4	7.9	8	12	11.9	15.9	15.8	20	20
5	L1	L2	L2	L1	L2	L2	L1	L2	L1	19.9	20	40	40	59.3	60.3	79.7	80.1	100.2	100.1
6	L1	L2	L2	L2	L1	L2	L2	L1	L2	11.9	12	24	23.9	36.1	35.8	48	47.9	59.7	59.8
7	L2	L1	L2	L2	L1	L1	L2	L2	L1	28.1	27.9	59.9	56.5	84	83.5	112.4	112.6	140.2	139.2
8	L2	L1	L2	L1	L2	L2	L2	L1	L1	28	27.8	56	56	84.1	83.8	111.5	111.8	140.2	140
9	L2	L1	L1	L2	L2	L2	L1	L2	L2	28	27.9	56.1	55.9	84.1	84	111.7	111.8	139.6	140.1
10	L2	L2	L2	L1	L1	L1	L1	L2	L2	28	27.9	56.1	55.8	83.5	83.9	112.2	112	140.2	139.6
11	L2	L2	L1	L2	L1	L2	L1	L1	L1	20	19.9	40.1	39.8	60	59.8	79.6	80.3	100.3	99.9
12	L2	L2	L1	L1	L2	L1	L2	L1	L2	12	12	23.9	23.9	36	35.8	47.6	47.8	60.5	59.6

制御因子がA~Iで、入力因子として20~100を考えます。今回は、誤差因子がありません。
（誤差因子のない例題はあまり品質工学としては良くない気がします・・・・）
また、各入力因子に対して、２つのデータが取得されています。

上記のように一次式のモデルを考えますが、品質工学ではなるべく避けるべきモデルと言われることが多い様です。（なるべく、エネルギー保存則が成り立つモデルで考えるべき。）

なにはともあれ、Rで解いていきたいと思います。

前回と同様に、Rで処理しやすいように、以下のようにcsvファイルを作成します。

Output.csv

X	Y1	Y2	Y3	Y4	Y5	Y6	Y7	Y8	Y9	Y10	Y11	Y12
20	19.9	11.9	27.9	3.9	19.9	11.9	28.1	28	28	28	20	12
20	20	11.8	28	4	20	12	27.9	27.8	27.9	27.9	19.9	12.9
40	39.9	24	55.6	7.9	40	24	55.9	56	56.1	56.1	40.1	23.9
40	40.1	23.7	55.9	8	40	23.9	56.5	56	55.9	55.8	39.8	23.9
60	60.1	35.6	84.1	12	59.3	36.1	84	84.1	84.1	83.5	60	36
60	59.9	36.1	83.9	11.9	60.3	35.8	83.5	83.8	84	83.9	59.8	35.8
80	80	47.7	112.5	15.9	79.7	48	112.4	111.5	111.7	112.2	79.6	47.6
80	80.2	47.9	111.9	15.8	80.1	47.9	112.6	111.8	111.8	112	80.3	47.8
100	100	59.9	140.6	20	100.2	59.7	140.2	140.2	139.6	140.2	100.3	60.5
100	100.5	60	140.7	20	100.1	59.8	140	140	140.1	139.6	99.9	59.6

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

※Eは誤差と勘違いされるため、あえて抜かしています。

■Rのプログラム

下記が実際のRのコードです。

#作業ディレクトリの設定
setwd("作業ディレクトリを設定")

#データの読み出し
LT<-read.table("Latin_Table.csv",header = TRUE,sep = ",")
str(LT)#読み出しデータ構造の確認
summary(LT)

Output<-read.table("Output.csv",header = TRUE,sep = ",")
str(Output)#読み出しデータ構造の確認
summary(Output)

#線形回帰
colname <-colnames(Output)#headerの名前を抽出
beta <- c(NULL)
y_est <- c(NULL)
for(i in 1:nrow(LT)){
  fm <- paste(colname[1+i],"~",colname[1])#各formulaを文字列から作成
  Output.lm <- lm(fm,data=Output)#各モデルの線形回帰を実施
  beta <- c(beta,Output.lm$coefficients[2])#傾きを抽出しベクトルに格納
  y_est <-cbind(y_est,Output.lm$fitted.values)#推定値をデータフレームに格納
}

#S_riの計算
Mbar<-apply(Output,2,mean)
SUM_Sm <- (Output[,1] - Mbar[1])
SUM_Sm <- SUM_Sm*SUM_Sm
SUM_Sm <- sum(SUM_Sm)/2
S_ri <- beta*beta*SUM_Sm

#V_eiの計算
Output_Y <- Output[,colnames(Output) != colname[1]]#入力因子を削除したデータフレームを用意
S_ei <- y_est - Output_Y
S_ei <- S_ei*S_ei
S_ei <- apply(S_ei,2,sum)
V_ei <- S_ei/8

#感度の計算
Sens <- (S_ri - V_ei)/SUM_Sm
Sens <- 10*log10(Sens)
Sens

#SN比の計算
SN <- (S_ri - V_ei)/SUM_Sm/V_ei
SN <- 10*log10(SN)
SN

#感度の分散分析
Sens_Bunsan <- cbind(LT,Sens)
Sens.lm <- lm(Sens~A+B+C+D+F+G+H+I+J,data = Sens_Bunsan)
anova(Sens.lm)

#感度の要因効果図の作成
par(mfrow=c(1,ncol(LT)))
sens_max <- max(Sens)
sens_min <- -5
fnames<-colnames(Sens_Bunsan)

for(i in 1:ncol(LT)){
  fc <- factor(Sens_Bunsan[,i])
  x_tmp <- 1:length(levels(fc))
  y_tmp <- tapply(Sens_Bunsan$Sens,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")
}

#SN比の分散分析
SN_Bunsan <- cbind(LT,SN)
SN.lm <- lm(SN~A+B+C+D+F+G+H+I+J,data = SN_Bunsan)
anova(SN.lm)

#SN比の要因効果図の作成
tmppar<-par(no.readonly=TRUE)
par(mfrow=c(1,ncol(LT)))
SN_max <- 15
SN_min <- 11
fnames<-colnames(SN_Bunsan)

for(i in 1:ncol(LT)){
  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")
}

結果は以下のようになります。

Analysis of Variance Table

Response: Sens
          Df Sum Sq Mean Sq F value  Pr(>F)
A          1 61.143  61.143 16.7743 0.05476 .
B          1 61.970  61.970 17.0010 0.05409 .
C          1 61.591  61.591 16.8972 0.05440 .
D          1  3.687   3.687  1.0116 0.42043
F          1  3.616   3.616  0.9919 0.42421
G          1  3.568   3.568  0.9788 0.42678
H          1 76.748  76.748 21.0554 0.04436 *
I          1  3.723   3.723  1.0213 0.41860
J          1  1.057   1.057  0.2900 0.64411
Residuals  2  7.290   3.645

Analysis of Variance Table

Response: SN
          Df Sum Sq Mean Sq F value Pr(>F)
A          1  0.759   0.759  0.0342 0.8703
B          1 26.516  26.516  1.1952 0.3884
C          1  4.926   4.926  0.2220 0.6839
D          1  0.705   0.705  0.0318 0.8750
F          1  1.024   1.024  0.0462 0.8498
G          1 13.722  13.722  0.6185 0.5140
H          1 33.273  33.273  1.4998 0.3454
I          1  6.035   6.035  0.2720 0.6540
J          1  4.130   4.130  0.1861 0.7082
Residuals  2 44.371  22.186   

それぞれのグラフは以下になります。
（クリックして拡大してみてください）

SN比は、グラフと分散分析表から分かるように、HとBが大きいのが分かります。
感度は、H,A,B,Cが大きいのが分かります。

これらの因子をもとに、パラメータ設計を実施し、確認試験を実施するのが一般的です。

本ページは下の書籍を参考に自習用にまとめています。



他の参考書籍です



関連ページ
静特性のパラメータ設計について
多因子--応答の最適化について
多因子-2水準-L16直交表の場合
２因子ー質的因子の場合
１因子ー量的因子ー二次モデルの場合
１因子ー量的因子ー一次モデルの場合
１因子ー質的因子の場合