相関係数行列の行列式に基づいた多変量管理図の研究
M2011MM043 近藤佑亮
指導教員：松田眞一

1

はじめに

管理図は製造工程が安定した管理状態下にあるかどうかを判断するための方法であり，W.A.Shewhart によって考えられた方法である。管理するべき品質特性が複数ある場合，多変量管理図が用いられる。各品質特性の平均や分散の変化を検出するための多変量管理図の研究は多くなされているが，相関関係に着目した多変量管理図の研究はあまりなされていない。安井ら [8] は相関関係の管理に標本相関係数行列の行列式に基づく管理図 (|R| 管理図) を提案している。
本研究では安井ら [8] の |R| 管理図と同様に標本相関係
数行列の行列式に基づく管理図として，ブートストラッ
プ |R| 管理図，|R| 累積和管理図を提案すると共に，その評価を行う。また先行研究で秋月 [1] は相関の変化を検出する管理図として，2 変量のデータに対して移動相関管理図を提案している。本研究では多変量データに対して移動相関係数行列の行列式による管理図を提案し，その評価を行う。

2

れる標本相関係数行列 RI の行列式の分布の (λ1 · · · λp )倍になると主張している。したがって，|R| 管理図における管理限界線は |RI | の分布におけるパーセント点の(λ1 · · · λp ) 倍で求めることができる。
各管理項目における特性値すべてが互いに独立である
とき |Ri | が 1 に近くなり，相関が強くなるほど 0 に近くなることを考慮して打点する統計量を − log |Ri | とする。
管理限界線を求めるのに − log |RI | のパーセント点が必要となるが，− log |RI | の分布を数学的に解くのは困難である。従って乱数を用いて元データと同じ標本の大きさをもつ 100000 群のデータから統計量を計算して，それらの経験分布から − log |RI | のパーセント点を求める。
|R| 管理図の信頼水準 100 × (1 − α)% の管理限界線の計算は以下の手順で行う。
m∑1
1. 母相関係数行列の推定値 P̂Σ = m
Ri を計算し，
i=1

その行列式 |P̂Σ | を計算する。

2. θI = − log |RI | の経験分布を乱数を用いて求める。
3. 上側管理限界線 (UCL)，下側管理限界線 (LCL) はそ
れぞれ，

管理図

管理図は，工程が管理状態であるかどうかを判断するた
めに W.A.Shewhart によって発表された手法である。管理図は一本の中心線 (CL) と上側管理限界線 (UCL)，下側管理限界線 (LCL) からなり，群ごとに統計量を打点していくことで，様々な判定ルール (例えば管理限界線の外に打点された) で工程が統計的管理下に置かれているかを判断する方法である。詳細は Montgomery[5] を参照のこと。

U CL = − log |P̂Σ | + θI,(1−α/2)
LCL = − log |P̂Σ | + θI,(α/2)となる。ここで θI,(1−α/2) および θI,(α/2) はそれぞれ θI の 1 − α/2 パーセント点，α/2 パーセント点である。

4

ブートストラップ法

ブートストラップ法は Efron[3] によって提案された数値計算法である。理論に基づく解析を，計算機を用いた反復
今 m 個の群を持ち標本の大きさが n であり，p 個 計算で置き換えたものであり，母集団に特定の分布形を仮の品質特性を持つデータベクトルを管理することを考 定せずに統計量の標準誤差や信頼区間を求めることがでえる。第 i 群に属する第 j 番目のデータベクトルの第 きる。ここで説明のために観測値 X を用いた工程管理を
k 要素を xijk で表すとき，データベクトルは xij = 考える。X の分布の形は未知であるが，工程を特徴づける(xij1 , · · · , xijk , · · · , xijp ) である。xij は i，j に対して互 パラメータ θ を推定したいとする。θ としては平均，標準
いに独立で同一な多変量正規分布 Np (µ, Σ) に従うとす 偏差などが考えられる。母集団から n 個の標本を抽出し，る。安井ら [8] は相関構造を示す統計量として，標本相関 それらを X1 , · · · , Xn とする。また関心のあるパラメータ
係数行列 Ri の行列式，すなわち，|Ri | を打点する管理図 θ の推定値を θ̂ とする。ブートストラップ法では観測されを提案した。
た標本 X1 , · · · , Xn から無作為復元抽出を行い，B 個の母相関係数行列 PΣ の固有値を λ1 , · · · , λp とするとき， ブートストラップ標本 Xb∗ = (X1∗ , · · · , Xn∗ ), b = 1, · · · , B
母相関係数行列の行列式は単位行列 I を用いて，を作成する。この B 個のブートストラップ標本に対して，それぞれブートストラップ推定値 θ̂∗ を計算する。このpp∏∏
ようにして計算された B 個のブートストラップ推定値|PΣ | =λl =λl |I|∗l=1
l=1
θ̂1∗ , · · · , θ̂Bは推定される θ̂ の標本分布を与える。

3|R| 管理図で 表 さ れ る 。λ1 , · · · , λp に 対 応 す る 固 有 ベ ク ト ル
p1 , · · · , pp からなる直交行列 P = (p1 , · · · , pp ) を用いて，Ri についても直交行列 P で変換すれば，安井ら [8]は行列式 |Ri | の分布は多変量正規分布 Np (0, I) から得ら

5

ブートストラップ |R| 管理図
工程が管理状態にあるとき観測値がモデル

xijk = µk + εijk

に従うと仮定する。ここで µk , k = 1, · · · , p は各品質特 6 移動相関係数行列の行列式による管理図
性 k の母平均である。また εijk は誤差項である。そして，
秋月 [1] は 2 変量のデータに対して移動相関を管理する
Seppala et al.[6] のバランス・サブグループ・ブートスト
移動相関管理図を提案した。本研究ではこれを多変量にラップ法を用いて，ブートストラップ |R| 管理図の管理
拡張したものとして移動相関係数行列の行列式による管限界線を構成する。その手順は以下のようになる。
理図を提案する。その構成法は以下の手順で行う。
1. 元の標本の残差
1. 群の数が m，標本の大きさが 1，特性の数が p であ

e1 = e11 = (e111 , · · · , e112 , · · · , e11p )e2 = e12 = (e121 , · · · , e122 , · · · , e12p )..
.
eN = emn = (emn1 , · · · , emn2 , · · · , emnp )を計算する。ここで N = mn であり，eijk = xijk −x̄i·k である。x̄i·k は群ごとの特性それぞれの標本平均である。

2. B = M m ≥ 10000 となるような整数 M を決める。
(M は計算時間を考慮して決める。)
3. N 個の残差ベクトルを M 回複製する。
e1 , · · · , e1 , e2 , · · · , e2 , · · · , eN , · · · , eN|
| {z } | {z }
{z
}
M 回
M 回
M 回
4. M N 個の残差ベクトルを無作為に並び替えて，それぞれの標本の大きさが n となる B 個の残差ベクトルの群を作る。
e∗11 , · · · , e∗1n , e∗21 , · · · , e∗2n , · · · , e∗B1 , · · · , e∗Bn
|
|
{z
} |
{z
}
{z
}
e∗
1

e∗
2

e∗
B

∗
5. ブートストラップ標本 Xbjk
, b = 1, · · · , B は µk のm∑1¯··k = m
推定量 x̄
x̄i·k を用いて
i=1∗¯··k + ae∗bjk
Xbjk
= x̄
√
となる。ここで a = n/(n − 1), n > 1 である。
∗∗
6. ブートストラップ標本 Xb1, · · · , Xbn
からブートス∗∗
トラップ統計量 Tb = − log |Rb | を計算する。ここで Rb∗ はブートストラップ標本の標本相関係数行列
である。

7. パーセンタイル法により管理限界線を構成する。ブートストラップ |R| 管理図における信頼水準 100 × (1 −α)% の上側管理限界線 (UCL) はブートストラップ
統計量 Tb∗ , b = 1, · · · , B の 1 − α/2 パーセント点∗
Tb,(1−α/2)であり，下側管理限界線 (LCL) は α/2 パー∗
セント点 Tb,(α/2)となる。また中心線 (CL) は元の標本の相関係数行列の行列式を対数変換したものの平均である。
∗
U CL = Tb,(1−α/2)

1 ∑− log |Ri |
m i=1
m

CL =

∗
LCL = Tb,(α/2)

るデータ xi = (xi1 , xi2 , · · · , xip ), i = 1, 2, · · · , m が得られたとする。(時系列データであるならば，前日との差分を取る等を行い，データが過去に依存しな
いようにする)2. あるブートストラップ反復係数 B を与えて，標本ベ
クトル xi , i = 1, 2, · · · , m に対して無作為復元抽出を行い，ブートストラップ標本 x∗b , b = 1, 2, · · · , Bを作成する。
3. ブートストラップ標本 x∗b , b = 1, 2, · · · , B を用いて，ウィンドウ幅 w の移動相関係数行列を計算する。
 ∗(j)∗(j)∗(j)
r1,1 r1,2 · · · r1,p


..
∗(j)


r
.


2,2
Rj∗ = 

..


∗(j)
. rp−1,p 

∗(j)
rp,p
∗(j)

′

ただし，rk,k′ , k = 1, 2, · · · , p, k = 1, 2, · · · , p, j =w, w + 1, · · · , B は，ブートストラップ標本から計算′される第 j 時点での第 k, k 成分間の移動相関であり，(x∗bk , x∗bk′ ), b = j − w + 1, j − w + 2, · · · , j の標
本相関係数である。
4. Rj∗ の行列式 |Rj∗ | を計算し，それを対数変換したものをブートストラップ統計量 Tj∗ = − log |Rj∗ |, j =w, w + 1, · · · , B とする。
5. 移動相関係数行列の行列式による管理図における信頼水準 100 × (1 − α)% の上側管理限界線 (UCL) はブートストラップ統計量 Tj∗ の 1 − α/2 パーセント∗
点 Tj,(1−α/2)であり，下側管理限界線 (LCL) は α/2∗
パーセント点 Tj,(α/2)となる。中心線 (CL) は元のデータの移動相関係数行列の行列式を対数変換したものの平均である。
∗
U CL = Tj,(1−α/2)

CL =

1
m−w+1

m−w+1
∑

− log |Rj |j=1∗
LCL = Tj,(α/2)

実データによる移動相関係数行列の行列式による管理
図の適用例として，YAHOO!ファイナンス [7] の自動車メーカーの株価を用いる。用いたデータはトヨタ，日産，三菱，マツダ，ホンダ，富士，スズキ，ダイハツの 8 社である。使用したデータの範囲は 2008 年 1 月初めから 2011年 12 月末までで，1 日ごとの終値を用いる。管理限界線の作成には 2007 年のデータを用いる。
2 社 (トヨタ，日産) による結果を図 1a，図 1b に，8 社の結果を図 1c，図 1d に示す。全ての図で 2008 年後半に

管理外に飛び出している。これはアメリカでのサブプラ えてよい。ここで新たな統計量 Ti を導入する。その定義
イム問題，リーマン・ショックが原因と考えられる。図 1b は次のようになる。
の 2010 年初めの飛び出しはトヨタのリコール問題と考え{
0, Ci+ = 0
られ，その次の飛び出しは急激な円高が原因と考えられTi =++j, Ci+ ̸= 0, · · · , Ci−j+1
̸= 0, Ci−j=0る。図 1b，図 1c の 2011 年後半の飛び出しは欧州債務危機等から生じた世界同時株安が原因と考えられる。また，
j = 1, 2, · · · , i
図 1d の 2011 年前半の飛び出しは東日本大震災が原因と思われる。
つまり，Ti は累積和統計量 Ci+ が最後に 0 であった所か
ら何点経過したかを表しており，Ti を sprint length と呼ぶ。Distribution-Free 累積和管理図での管理限界線は Tiに応じて決められる。Ti = j, j = 1, 2, · · · , jmax である
Ci+ に対しては管理限界線 hj , j = 1, 2, · · · , jmax が用いられる。Ti が jmax を超える Ci+ に対しては管理限界線h∗ を用いる。ここで jmax は事前に決めておく定数であるが，Chatterjee and Qiu[2] は多くても 50 程で良いとしている。
|R| 累積和管理図の構成法の手順を簡単に説明する。
2009

2010

2011

2008

2009

2010

2011

0.0

0.6

0.5

0.8

1.0

1.0

1.5

2.0

1.2

2008

0

200

400

600

800

1000

0

200

400

day

600

800

1000

1. 目標とする平均的な Ti ，ASL(average sprint length)
を決める。

day

2009

2010

2011

2008

2009

2010

2. ブートストラップ法を用いて Cb , b = 1, 2, · · · , B を計算し，Tb の平均が ASL と十分近くなるように kを反復法により調整する。

2011

3.0

3

3.5

4

4.0

5

4.5

6

5.0

7

5.5

8

6.0

2008

9

(a) 2 社データ (ウィンドウ幅 200) (b) 2 社データ (ウィンドウ幅 50)

0

200

400

600

800

1000

day

0

200

400

600

800

1000

day

(c) 8 社データ (ウィンドウ幅 200) (d) 8 社データ (ウィンドウ幅 50)

図 1: 移動相関係数行列の行列式による管理図|R| 累積和管理図
Upper cusum−1

Cumulative Sum

1

管理図の一種として累積和統計量を管理対象とする累積和管理図がある。累積和管理図は比較的小さな変化に鋭敏であることが知られている。そこで相関係数行列の行列式 |R| に対して累積和統計量をとることが考えられるが，一般的な累積和管理図は正規分布の仮定を置いているため，そのままでは使用することができない。そのため Chatterjee and Qiu[2] の Distribution-Free 累積和
管理図を適用することで，|R| 累積和管理図を構成する。
累積和管理図における統計量は

073. 求められた k を用いて，目標とする ARL(averagesprint length) に十分近くなるように管理限界線{hj , j = 1, 2, · · · , jmax ; h∗ } をブートストラップ
法と反復計算により求める。
|R| 累積和管理図の例を図 2 に示す。図の縦軸の 0 より上に C + ，下に C − が打点されている。また図の曲線が管理限界線を表している。|R| 累積和管理図では C + と C −の k は異なる値をとるため C + と C − が共に 0 より大きい箇所がある。

−2

lower cusum

0

20

40

60

80

100

Subgroup number+
Ci+ = max(0, θi − k + Ci−1), C0+ = 0−
Ci− = max(0, −k − θi + Ci−1), C0− = 0となる。C + は正の方向の変化を，C − は負の方向の変化を検出するための統計量であり，それぞれ上側累積和統
計量，下側累積和統計量と呼ばれる。θi は通常の累積和管理図では標本値や標本平均が用いられるが，|R| 累積和
管理図では相関係数行列の行列式を用いる。以降の説明では C + のみを対象とするが，C − についても同様と考

図 2: |R| 累積和管理図

8

シミュレーション

ブートストラップ |R| 管理図，|R| 累積和管理図が有意水準を守っているかを検証する。特性 p = 2, 6 の 2 種類の多変量正規分布 N (µ, Σ) に従うデータを用いて，ブートストラップ |R| 管理図では有意水準 0.0027，|R| 累積和管理図では有意水準水準 0.0027 に対応する ARL = 370を満たす管理限界線を作成する (有意水準 α と ARL には
α = 1/ARL の関係がある)。また 2 変量データに対しては検出力も評価する。6 変量正規分布のパラメータ µ, Σ
は
 


1
1 0.41 0.64 0.79 0.58 0.35
 2 

1
0.64 0.43 0.19 0.48 
 


 2 

1
0.54 0.42 0.47  µ =  ,Σ = 

 3 

1
0.49 0.37 
 



 1 
1
0.30 

31

表 2: |R| 累積和管理図(a) ARLn2 変量データ

6 変量データ

20
50

299.8(114.4)
296.7(111.9)

313.4(145.7)
303.4(118.9)(b) 検出力 (2 変量データ)

n

ρ2

ρ3

20
50

40.4(31.0)
12.1( 4.9)

31.3(20.4)
9.2( 3.2)を用いる。2 変量正規分布では標準状態の相関係数を
ρ0 = 0.6 とする。また検出力の評価では相関係数をなものであることがわかった。また，シミュレーションにρ1 = 0.4, ρ2 = 0.5, ρ3 = 0.7, ρ4 = 0.8 の 4 種類に変よりブートストラップ |R| 管理図および，|R| 累積和管理化させることで検出力を確かめる。ここで検出力とは変図の管理限界線が第 1 種の過誤を守っていることを示す化が生じた後，平均何点目で検出できるかを表す。
ことができた。|R| 累積和管理図はブートストラップ |R|
ブートストラップ |R| 管理図において検証に用いる群管理図より検出力の点で非常に優れていた。しかし，|R|の数を 100000，管理限界線の作成で用いる群の数 m =累積和管理図は累積和統計量の計算や複数の管理限界線
30，標本の大きさ n = 20, 50，ブートストラップ反復係数があることなど，実際に使用する際には非常に複雑な手
M = 1000 として，上記の設定のもとでシミュレーショ
法であることが欠点であると言える。
ンを行う。また |R| 累積和管理図において，管理限界線の作成で用いる群の数 m = 30，標本の大きさ n = 20, 50， 10 おわりに目標とする ASL を 5 として，シミュレーションを行う。
この研究を通して，各管理図が有用なものであるとわシミュレーション回数は共に計算時間を考量して 100 回 かった。今後の課題としては，シミュレーションで用いとする。
るデータの吟味である。今回多変量データとして多変量シミュレーションの結果を表 1，2 に示す。表の括弧内 正規分布を用いたが，そのパラメータについての議論はの数字は標準偏差を表す。表 1a，2a の結果よりブートス 不十分であるため，しっかりと吟味をする必要がある。
トラップ |R| 管理図，|R| 累積和管理図ともにそれぞれの目標とする有意水準，ARL より少し甘い値となっている。 参考文献
また，どちらの管理図も第 1 種の過誤の点で大きな差は [1] 秋月良康：移動相関管理図の構成とその有用性の研究,ない。一方，検出力は表 1b，表 2b より，ブートストラッ
南山大学数理情報研究科修士論文, 2012.
プ |R| 管理図は n = 50 の時で相関係数の変化が ±0.2 で[2] Chatterjee, S. and Qiu, P. : Distribution-Free Cumuある ρ1 , ρ4 では早く検出できているが，相関係数の変化
lative Sum Control Charts Using Bootstrap-Basedが ±0.1 である ρ2 , ρ3 では検出力が悪い。一方，|R| 累積
Control Limits. The Annals of Applied Statistics, 3,
和管理図は非常に早く検出できている。
349–369, 2009.
表 1: ブートストラップ |R| 管理図[3] Efron, B. : Bootstrap methods: another look at the
jackknife. Annals of Statistics, 7, 1–26, 1979.

(a) 第 1 種の過誤

n2 変量データ

6 変量データ

2050

0.0032(0.0010)
0.0031(0.0009)

0.0032(0.0012)
0.0034(0.0010)[4] Hawkins, D. M. : Cumulative Sum Control Charting: An Underutilized SPC Tool. Quality Engineering, 5, 463-477, 1993.
[5] Montgomery, D. C. : Introduction to Statistical
Quality Control. John Wiley & Sons, 2001.

(b) 検出力 (2 変量データ)

n

ρ1

ρ2

ρ3

ρ4

20
50

59.8(30.1)
8.7( 2.9)

164.0(74.9)
52.1(25.0)

118.1(65.2)
34.4(16.8)

13.9(6.0)
2.6(0.5)[6] Seppala, T. and Moskowitz, H., Plante, R. and Tang,
J. : Statistical Process Control via the Subgroup
Bootstrap. Journal of the Quality Technology, 27,
139–153, 1995.
[7] YAHOO!ファイナンス :http://finance.yahoo.co.jp/ , 2012.

9まとめ

移動相関係数行列の行列式による管理図が一部の株価データに対してのみでしか検証を行っていないが，有用[8] 安井清一，尾島善一，仙石直人 : 相関構造の変化を検出するための多変量管理図. 日本品質管理学会第 80回研究発表会研究発表要旨集, 259–262, 2006.