続きました。
前回vmask法についてよくわからなかったのですが、この本を進めてくださった方から
「vMaskはパッケージがあるよ」
と教えていただけたので早速調べてみました。
vmask法について
vmaskは累積和管理図(CUSUM)で平均値の変化を知るために使われるものです。
とあるVマスクというV字を作って、引っかかったら変化アリというものです。
vマスクとは図中のプロット先端(vertex)から伸びているV字のもののことです。
plot点がoriginだとして、dという距離だけ離した位置にvertexを設置し、originから上下にhだけ離した距離を通る上下の線をVマスクとしています。
vmaskの設計はいくつか種類があり、参考論文によって設計が異なるようです。
以下で簡単な参考文献の違いだけ紹介しておきます。
統計的工程管理(仁科)での参考文献
Barnard(1959)
Johnson and Leone(1962)
パッケージ作者(Abbas Parchami)の参考文献
Montgomery, D.C. (1985)
Goel, A. L. (2011)
Goel, A. L. (1982)
(余談ですがアノニマスの仮面もvmaskというらしいですね)
早速パッケージインストールして使ってみる!
library(vMask)
ちなみにvmaskという名前のパッケージもあるのですがこちらはver=3.5でのみ使えるもののようです。
vMaskには5つの関数が入っておりますが、直観的に使いやすそうなmethod1を紹介します。
一変量のとき
Data = c( 324.925, 324.675, 324.725, 324.350, 325.350, 325.225,
324.125, 324.525, 325.225, 324.600, 324.625, 325.150,
328.325, 327.250, 327.825, 328.500, 326.675, 327.775,
326.875, 328.350 )
vMask.method1( data=Data, d=1, theta=45, s=0 )
管理内である値を緑色にplotし、InControlとして表示させています。
最終ポイントから計算してはずれであるものが赤にplot.
関数の中身が気になりますので覗いてみます。
vMask.method1 <-
function(
data, #Type of inputed data must be a matrix with dim=c(m,n)
mu0 = mean(data), #Target value for process mean, with defult mean(data)
d, #The distance between ...
theta = 14, #The angle (in degrees) between ...
sleep = 1 #Sleep time (in secound)
)
まず引数は
Data：これはnumericでもmatrixでも入力できます。
if( class(data) == "matrix" ){
m = dim(data)[1]#The size of samples
n = dim(data)[2] #The size of each subsample
RowMeans <- rowMeans(data)
}
if( class(data) == "numeric" | class(data) == "integer" ){
m = length(data)#The size of samples
n = 1 #The size of each subsample
RowMeans <- data
}
他のパラメータは、
mu0：平均値
dとthetaを設定
sleepは動かしてみると分かるのですが、plotしていく様がR上で確認できます。
これをsleep="PressEnter"に設定するとひとつずつplotしていく様子が確認できます。
これらのパラメータは実務で使う場合には、値が平均値から変わらないときなどは自分で決めてしまってもいいかもしれません。
検出力を高めたいときはthetaの角度をきつめに設定して外れたものを見つけやすくするなど。
本では
h = d*tan(theta)
tan(theta)=(data.i - data.mean) / 2σ
で設計されているようです。
ちなみにthetaは決めていないけどhは決めているというときには、
method2を使うと都合がよさそうです。
つづき
theta <- deg2rad(theta)
h = d * tan(theta)#The length of decision making
c = cumsum( RowMeans - mu0 )
入力されたパラメータを使って関数内部の引数に格納しています。
あとはplotする内容が書かれています。
色付けやVmaskのアームを書いたり上限加減を書いたり。
私のスクリプトでは全然描けなかったわけです...
関数内部を覗くだけでも勉強になります。
x.win = c(0, m+d+5)
y.win = c(min(c)-h*(m-1+d)/d, max(c)+h*(m-1+d)/d)
BackColor = rep(1,m) #For bg vector
for(i in 1:m) {
plot(x.win, y.win, type="n", #For plotting the window
xlab = "i", ylab = expression(c[i]),
main = list(paste("CUSUM control chart and V-Mask on point (", i, ", c", i, ")")) )
lines(1:m, c, type="b", col=1, cex=1.3 , pch=21, bg=BackColor )
points(i, c[i], type="p", cex=2 ,col=1, pch=21, bg=7)
x.Mask <- c(i+d,
1,
i+d+4,
i+d+4,
1 )
y.Mask <- c(c[i],
c[i]+h*(i-1+d)/d,
c[i]+h*(i-1+d)/d,
c[i]-h*(i-1+d)/d,
c[i]-h*(i-1+d)/d )
polygon( x.Mask, y.Mask,#Plotting the mask
border=2, col="coral1", lwd=1,
density=c(10, 20), angle=c(-45, 45) )
lines( c(i,i+d), c(c[i],c[i]), col="gray60", lty=2 )
lines( c(i,i), c(c[i]+h,c[i]-h), col="gray60", lty=2 )
# abline( h=mean(c), col=2, lty=3 )
for(j in 1:i){
LCL <- c[i] - h*(i-j+d)/d
UCL <- c[i] + h*(i-j+d)/d
# points( c(i,i), c(LCL,UCL), col="red" )
if( !(LCL<=c[j] & c[j]<=UCL) ){
for(jj in 1:i){
LCL <- c[i] - h*(i-jj+d)/d
UCL <- c[i] + h*(i-jj+d)/d
if( !(LCL<=c[jj] & c[jj]<=UCL) )
points(jj, c[jj], type="p", cex=1.6 ,col=1, pch=21, bg="turquoise1")
}
BackColor[i] = 2
break
} else{
BackColor[i] = "green3"
}
}
if(sleep == "PressEnter"){
cat ("Press [enter] in 'R Console' to continue")
line <- readline()
}else { Sys.sleep(sleep) } # Sleeping time for computer
}
return( list(
h = h,
c = c,
OutControl = which(BackColor==2), #Out-Control subsamples
InControl = which(BackColor=="green3") #In-Control subsamples
)
)
関数の紹介になってしまいましたが、やりたいことはこのパッケージを使えばできそうだと確認できましたので良かったと思います。
あと興味のある内容は重回帰解析といろいろ付録についている部分になります。
一旦本記事で区切るつもりですが、いつかつづくかもしれません。
つづく？