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

前回は、量的因子（１因子）の一次モデルを解いたので、今回は２次モデルの分散分析と区間推定を実践したいと思います。

■データでトライ

前回と同じデータを用いて、二次モデルの分散分析と区間推定にトライしたいと思います。

今回は、SrとSeに関しては、行列で表現すると以下のように表すことができるのを利用して計算しています。

Jnはやや特殊ですが、要素が1/n（平均値をとるため）の（nｘn）行列のマトリックスになります。

■scilabで解く

上記をscilabで解いた際のソースは以下です。

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

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

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

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

//分散分析表の計算
J = ones(n,n);
J = J*1/n;
Sr = (X*b - J*y)'*(X*b - J*y)
Se = (y-X*b)'*(y-X*b)
Vr = Sr/2
Ve = Se/17
f = Vr/Ve
p = 1 - cdff("PQ",f,2,17)

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

上記の結果区間推定は以下のようなグラフになります。

分散分析表は以下のようにかけます

要因	S	φ	V	F	p
因子 x（1次）	79.2486	2	39.6243	22.914426	0.000015
誤差	39.6243	17	1.7292294
	ST	φT

前回の一次モデルより、Veの値が小さいため、モデルがうまくあてはまっていることが分かります。

■追記 Rで解く

だんだんと自分の中で恒例になってきましたが、Rで解きたいと思います。
線形回帰の際は、lm関数を用いましたが、非線形回帰を行う際は、nls関数を使います。
lm関数から吐き出されるオブジェクトは、下記の点で便利でした。
・plot関数にいれると、細かいグラフができた
・anova関数に単体でいれると元データとの分散分析をしてくれる
・predict関数で区間推定を実施してくれた

しかし、nls関数は、上記のような機能はないようです。

詳細はこちらに詳しいです。

その代わりにscilabで示したように行列演算すれば良いのですが、面倒だったので、分散分析と区間推定は実施していません。（何かスマートな方法があれば教えて欲しいです。）
非線形回帰分析は実施ました。ご参考に。

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

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

#非線形回帰分析
d.nls <- nls(y~a*x^2+b*x+c,data = d,start = c(a=0,b=0,c=0))
summary(d.nls)

#グラフの描画
new <- data.frame(x = seq(200 , 260 , 1 ))
d.predict <- predict(d.nls,newdata = new)
plot(d$x,d$y,xlim=c(200,260), ylim=c(0,10),xlab="",ylab="")
par(new=T)
X = seq(200,260,1)
par(new=T)
plot(X, d.predict,type="l", xlim=c(200,260), ylim=c(0,10),col="red",xlab="",ylab="")
legend("topleft", legend = c("data","estimation"),pch=c("○","-"), col = c("black","red"))

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