「統計的工程管理」を読んでRで実装【vmask法】　#2(終)


[「統計的工程管理」を読んでRで実装　#1]()の続き。

続きました。
前回vmask法についてよくわからなかったのですが、この本を進めてくださった方から
**「vMaskはパッケージがあるよ」**
と教えていただけたので早速調べてみました。

#vmask法について

vmaskは累積和管理図(CUSUM)で平均値の変化を知るために使われるものです。
とあるVマスクというV字を作って、引っかかったら変化アリというものです。

![vmask.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/01ee8eb8-c1d1-c103-862b-933715d823e3.gif)


[画像参照元](https://v8doc.sas.com/sashtml/qc/chap12/sect16.htm)

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というらしいですね)

![Spot-wholesale-V-movie-theme-V-face-mask-vendetta-mask-V-mask-shape-hacker-V-blush.jpg_640x640.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/05c249c0-ee84-9fc6-197d-f697f0cad498.jpeg)

#早速パッケージインストールして使ってみる!

```R
library(vMask)
```

ちなみにvmaskという名前のパッケージもあるのですがこちらはver=3.5でのみ使えるもののようです。

vMaskには5つの関数が入っておりますが、直観的に使いやすそうなmethod1を紹介します。

**一変量のとき**

```R
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 )
```
![vm1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/b7d3ef9b-bb99-24dc-c259-3df9b36711fc.png)

![vm2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/c9bca479-5348-5c68-ae08-bbf4d3cfa325.png)

管理内である値を緑色にplotし、InControlとして表示させています。
最終ポイントから計算してはずれであるものが赤にplot.

関数の中身が気になりますので覗いてみます。

```R
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でも入力できます。

```R
  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を使うと都合がよさそうです。

つづき

```R
  theta <- deg2rad(theta)
  h = d * tan(theta)#The length of decision making
  c = cumsum( RowMeans - mu0 )
```
入力されたパラメータを使って関数内部の引数に格納しています。


あとはplotする内容が書かれています。
色付けやVmaskのアームを書いたり上限加減を書いたり。
私のスクリプトでは全然描けなかったわけです...
関数内部を覗くだけでも勉強になります。

```R
  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
        )  
   )
```


関数の紹介になってしまいましたが、やりたいことはこのパッケージを使えばできそうだと確認できましたので良かったと思います。

あと興味のある内容は重回帰解析といろいろ付録についている部分になります。
一旦本記事で区切るつもりですが、いつかつづくかもしれません。

つづく？