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

個人的に目的だった、タグチメゾットを使ったパラメータ設計についてまとめたいと思います。

誤差因子を取り上げることで、ロバストな設計を目指す点が、今まで復習してきた実験計画法とは大きく異る点です。

今回は静特性のパラメータ設計についてまとめたいと思います。

■静特性のパラメータ設計

静特性のパラメータ設計では、目標値が特定の値である望目特性、値が負にならず小さいほどよい望小特性、値が大きいほどよい望大特性、値が0であることを望む０望目特性などに大別されます。

通常パラメータ設計を行う場合は、以下の２つのステップでパラメータ設計を実施していきます。

ステップ１：誤差因子の水準変動に対してロバストな制御因子の水準を見出す
ステップ２：ステップ１でのロバストネスとは無関係に応答の平均を調整しうる因子を見出し、それにより平均を目標値に近づける。

ばらつきを小さくするのに、タグチメゾットでは、以下のSN比を大きくする制御因子の水準を探し出します。

また、SN比には影響を与えないが、感度μ^2に影響を与える因子を見出し平均値の調整を実施していきます。

上記のSN比や感度は、望目特性、望小特性、望大特性、０望目特性において以下のように異なる定義がされる点に注意が必要です。

１）望目特性

２）望小特性

３）望大特性

４）望０特性

■実際のデータで計算

実例として、望目特性のパラメータ設計を実践したいと思います。

									S	1	2	3	3	1	2	2	3	1
									R	1	2	3	2	3	1	3	1	2
									Q	1	2	3	1	2	3	1	2	3
									P	1	1	1	2	2	2	3	3	3
	A	B	C	D	F	G	H	I	J
1	1	1	1	1	1	1	1	1	1	15.9	14.5	13.1	14.4	13	13.8	12.9	13.7	12.3
2	1	1	1	1	1	2	2	2	2	15.7	13.5	11.4	13.4	11.2	13.6	11.1	13.5	11.3
3	1	1	2	2	2	1	1	1	2	18.3	16.8	15.3	16	14.5	17.5	13.8	16.8	15.3
4	1	2	1	2	2	1	2	2	1	17.7	16.3	14.9	16.9	15.5	16.4	16.2	17	15.6
5	1	2	2	1	2	2	1	2	1	16.3	15.6	14.8	15.6	14.8	16.3	14.8	16.3	15.6
6	1	2	2	2	1	2	2	1	2	18.1	16.6	15.1	16.6	15.1	18.1	15.1	18.1	16.6
7	2	1	2	2	1	1	2	2	1	18.2	17.4	16.7	16.7	15.9	17.4	15.2	16.7	15.9
8	2	1	2	1	2	2	2	1	1	16.2	15.4	14.7	14.7	13.9	15.4	13.2	14.7	13.9
9	2	1	1	2	2	2	1	2	2	17.3	15.1	13	15	12.9	15.2	12.8	15.1	13
10	2	2	2	1	1	1	1	2	2	18.7	17.2	15.7	17.2	15.7	18.7	15.7	18.7	17.2
11	2	2	1	2	1	2	1	1	1	17.3	15.9	14.5	16.5	15.1	16	15.8	16.6	15.2
12	2	2	1	1	2	1	2	1	2	18	15.9	13.7	16.5	14.4	16.7	15	17.4	15.2

上記のようなA〜Jの制御因子をもち、S〜Pのような誤差因子をもつ場合を考えてみます。
それぞれの交差する部分が、出力値となります。

Rで解きやすいように以下のように２つに分けたcsvファイルを用意して、データのSN比と感度を求め、それぞれの分散分析を実施したいと思います。

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

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("作業ディレクトリを指定")

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

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

#望目特性　SN比の計算
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比の分散分析
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)

#SN比の要因効果図の作成
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.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")
}

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

> anova(SN.lm)
Analysis of Variance Table

Response: SN
          Df Sum Sq Mean Sq  F value   Pr(>F)
A          1  0.099   0.099   0.3403 0.618651
B          1 17.234  17.234  59.5159 0.016390 *
C          1  9.842   9.842  33.9876 0.028185 *
D          1  0.832   0.832   2.8737 0.232122
F          1  0.840   0.840   2.9009 0.230642
G          1  0.017   0.017   0.0586 0.831283
H          1  0.213   0.213   0.7353 0.481522
I          1  0.662   0.662   2.2859 0.269693
J          1 55.707  55.707 192.3830 0.005158 **
Residuals  2  0.579   0.290                  

> anova(Sens.lm)
Analysis of Variance Table

Response: Sensitivity
          Df  Sum Sq Mean Sq   F value    Pr(>F)  
A          1 0.41370 0.41370  581.3289 0.0017158 **
B          1 2.36922 2.36922 3329.1803 0.0003002 ***
C          1 1.73940 1.73940 2444.1792 0.0004089 ***
D          1 1.07440 1.07440 1509.7284 0.0006617 ***
F          1 0.00486 0.00486    6.8237 0.1206041  
G          1 0.96297 0.96297 1353.1482 0.0007382 ***
H          1 0.00026 0.00026    0.3584 0.6101908  
I          1 0.00140 0.00140    1.9734 0.2952651  
J          1 0.00204 0.00204    2.8676 0.2324580  
Residuals  2 0.00142 0.00071       

 SN比の要因効果図

感度の要因効果図

要因効果図や分散分析表から以下のことが分かります。

最もSN比に効果が大きいのは、因子Jであり、次に因子B,Cである。これらの因子がロバスト性に重要な因子であることが分かります。
一方で、感度に影響が大きいのは、因子B,C,D,G,Aとなります。ただし、因子B,CはSN比を大きくするために用いるので、調整因子としてD,G,Aを使うのが望ましいということになります。

J,B,Cの因子でSN比を最大にし、D,G,Aの因子で欲しい特性値に近づけることで、パラメータ設計を実施することができます。

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



他の参考書籍です



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