実験計画法 １因子-量的因子-一次モデル に関するメモ （Design of experiments & Taguchi method）

業務で必要になる機会が多いので、実験計画法やタグチメゾットについて復習しています。以下は、個人的に備忘としてメモとして残します。

■分散分析と区間推定について

限られた実験から、因子の違いの効果があるかを判定したり、母集団のバラつきを求めれることは多くあります。
そういった際に、分散分析やt分布を使った母集団の区間推定などを利用します。
１因子の質的因子の場合は、分散分析も区間推定もよく解説があるので、１因子の量的因子の一次モデルについて備忘として残します。

■１因子/量的因子/一次モデルについて

母平均μ(x)が、量的因子xについて以下のように表現されるとします。

出力yは一次モデルとして以下のように表現できるとします。

上記を行列表現すると以下のようになります。

得られた実験のデータから未知数βは最小二乗法で以下のように求めることができますね。

■実際のデータでトライ

ここでは、実際の例を使って推定した１次モデルを使って、分散分析と区間推定を実践したいと思います。（分散分析の結果一次モデルの推定が正しいことを検定できないが、一次の効果があるかないかを示すことはできる点に注意してください）

以下のようなデータがとれたとします。（一次モデルをあてはめるのは間違いな気がしますが例題として先に進めます）

因子	入力x	出力y
1	200	4.5
	200	3.5
	200	4.4
	200	2.6
	200	6.5
2	220	4.2
	220	4.2
	220	1.7
	220	3.4
	220	3.7
3	240	5.7
	240	3.4
	240	5.1
	240	3.9
	240	4
4	260	6.1
	260	7.5
	260	10.3
	260	8.9
	260	10.1

■最小二乗法を解いてフィッティング

上記のデータを用いてβをとくと、以下のようになります。

β0 = -10.708

β1 = 0.0691

■要因効果の分散分析による検定（１因子/量的因子/一次モデル）

上記の結果を用いて下記の分散分析（ANOVA）表を作成することを考えます。

要因	S	φ	V	F	p
因子 x（1次）	Sr	φr	Vr=Sr/φr	Vr/Ve
誤差	Se	φe	Ve=Se/φe
	ST	φT

各式は以下のようにかけます。ここでは、μは近似した１次方程式から算出するので注意してください。

分散分析表の意味は以下のように自分は考えています

要因	S	φ	V	F	p
因子 x（1次）	因子の影響を説明する（因子（今回の場合は一次モデル）に対してのバラつき）	自由度 （１次だから１）	自由度で割る	因子の効果が誤差の何倍かを示す	p値
誤差	誤差の影響を説明する（残差のバラつき）	自由度	自由度で割る
	ST	φT

実際に計算すると以下のようになります。

要因	S	φ	V	F	p
因子 x（1次）	47.748	1	47.748	14.113	0.0014
誤差	60.897	18	3.383
	ST	φT

上記より、一次の効果があることはいえますが、一次モデルで十分に説明できることを保障している訳ではないことに注意しましょう。

■区間推定

区間推定は以下の式から導出されます。

（分散分析と同じように、区間推定も１次式に対しての区間推定となるので注意が必要）

95%の信頼区間は以下のようにかけます。

予測区間は以下です。

行列で書いた方が覚えやすいので、以下に行列表現での信頼区間と予測区間を示します。

■scilabで解く

では、上記説明してきたこをscilabで解くと以下のようになります。

//実験データ入力
X = [
1 200 4.5
1 200 3.5
1 200 4.4
1 200 2.6
1 200 6.5
1 220 4.2
1 220 4.2
1 220 1.7
1 220 3.4
1 220 3.7
1 240 5.7
1 240 3.4
1 240 5.1
1 240 3.9
1 240 4
1 260 6.1
1 260 7.5
1 260 10.3
1 260 8.9
1 260 10.1
];

y = X(:,3);
X(:,3)=[];

//グラフ描画用行列の作成
a = [200:260];
b = ones(61,1);
x=[b a'];

//最小二乗法
b=inv(X'*X)*X'*y

//Veの計算
Y = X*b; //一次モデルを計算
Se = y-Y; //残差を計算
Se = Se*Se'; //Seを計算①
Se = diag(Se); //Seを計算②
Se = sum(Se); //Seを計算③
Ve = Se/18; //Veを計算

//t分布の区間推定
for i=1:61
I=[1,i+199];
t(i,1)=I*b+cdft("T",18,0.975,0.025)*sqrt((I*inv(X'*X)*I')*Ve);//上側信頼区間
t(i,2)=I*b-cdft("T",18,0.975,0.025)*sqrt((I*inv(X'*X)*I')*Ve);//下側信頼区間
t(i,3)=I*b+cdft("T",18,0.975,0.025)*sqrt((1+I*inv(X'*X)*I')*Ve);//上側予測区間
t(i,4)=I*b-cdft("T",18,0.975,0.025)*sqrt((1+I*inv(X'*X)*I')*Ve);//下側予測区間
end

//グラフ出力
plot(X(:,2),y,'*r');
plot(x(:,2),x*b,'b');
plot(x(:,2),t(:,1),'g');
plot(x(:,2),t(:,2),'g');
plot(x(:,2),t(:,3),'k');
plot(x(:,2),t(:,4),'k');
legend(["data","linear_model","上側信頼区間","下側信頼区間","上側予測区間","下側予測区間"],2);
a=gca(); //アクティブな軸のオブジェクトを取得
a.data_bounds(:,1)=[200;260];
a.data_bounds(:,2)=[0;10];
xlabel('input x');
ylabel('output y');

１次モデルであてはめるのはナンセンスですが、１次の効果があるとデータからはいえることが分かりました。
また、求めた一次モデルから信頼区間を計算すると上記のようになります。
予測区間がかなり広いのが気になりますが、きちんと求まりました。

■追記　Rで解く

線形回帰推定と、その得られたモデルとの分散分析を実施したいと思います。
Rは、これを実行するのに関数がすでに用意されているので、あてはめるだけです。
（グラフの描画の部分がプログラムとしては長くなってしまいました・・・・）

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

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

#線形回帰分析
d.lm <- lm(y~x,data = d)
summary(d.lm)

#分散分析表の表示
anova(d.lm)

#区間推定の実行
new <- data.frame(x = seq(200 , 260 , 1 ))
d.conf <- predict(d.lm,newdata = new,interval = "confidence",level=0.95)
new <- data.frame(x = seq(200 , 260 , 1 ))
d.pre <- predict(d.lm,newdata = new,interval = "prediction",level=0.95)

#グラフの描画
plot(d$x,d$y,xlim=c(200,260), ylim=c(0,10),xlab="",ylab="")
par(new=T)
abline(d.lm,lwd=1,col="blue")
X = seq(200,260,1)
par(new=T)
plot(X, d.conf[,2] ,type="l", xlim=c(200,260), ylim=c(0,10),col="red",xlab="",ylab="")
par(new=T)
plot(X, d.conf[,3] ,type="l", xlim=c(200,260), ylim=c(0,10),col="red",xlab="",ylab="")
par(new=T)
plot(X, d.pre[,2] ,type="l", xlim=c(200,260), ylim=c(0,10),col="green",xlab="",ylab="")
par(new=T)
plot(X, d.pre[,3] ,type="l", xlim=c(200,260), ylim=c(0,10),col="green",xlab="x",ylab="y")
legend("topleft", legend = c("data","confidence","predoction"),pch=c("○","-","-"), col = c("black","red","green"))

#グラフの描画（推定モデルの確認）
par(mfrow=c(2,2))
plot(d.lm)

結果は以下のような感じです。

Call:
lm(formula = y ~ x, data = d)

Residuals:
   Min     1Q Median     3Q    Max
-2.794 -1.110 -0.294  1.313  3.388

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -10.70800    4.25044  -2.519  0.02143 *
x             0.06910    0.01839   3.757  0.00144 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.839 on 18 degrees of freedom
Multiple R-squared:  0.4395, Adjusted R-squared:  0.4083
F-statistic: 14.11 on 1 and 18 DF,  p-value: 0.001444

Analysis of Variance Table

Response: y

          Df Sum Sq Mean Sq F value   Pr(>F)   

x          1 47.748  47.748  14.113 0.001444 **

Residuals 18 60.897   3.383                    

---

グラフは以下のようになります。

上記のグラフは、残差のプロットとQ-Qグラフなどがあり、推定式が妥当かどうかを判断するのに使うことができます。

こういうグラフがすぐに作れるところは地味ながら嬉しい点です。

二次モデルについてもトライしているので、ご参考に