多群指数モデルにおけるすべての平均相違の多重比較法
2009SE056 廣瀬由幸 2009SE178 森下史章
指導教員：白石高章=はじめに

1y1y n−1 e− 22n Γ(n)= (χ22n の密度関数)

３つ以上の群で, 平均に関する場合に行う多重比較法について考察する.

X1 , X2 , · · · , Xn1 は互いに独立で同一の密度関数 fX (x) =(−x/µ1 )を も つ 指 数 分 布 EX(1/µ1 ) に 従 い,ferroni 法（ボンフェローニ）, Holm 法（ホルム）, Dunnett (1/µ1 )e
Y,Y,···,
Yは互いに独立で同一の密度関数
fY (x) =n2法（ダネット), Tukey-Kramer 法（テューキー・クレーマー） 1 2(−x/µ2 )(1/µ2 )eをもつ指数分布 EX(1/µ2 ) に従うとする.などがある.
また,X,Y(i=1, · · · , n1 , j = 1, · · · , n2 ) も互いに独立.
ij
この中で今回の解析では, 正確な検定をするためボン
2n12
フェローニ法を用いて, 多重比較検定と同時信頼区間を調 命題 1 χ2n1 は自由度 2n1 のカイ 2 乗分布を表し, F2n2 はべる. また, 地震の発生回数はポアソン分布に従うが, 地震 自由度 (2n1 , 2n2 ) の F 分布を表すものとする. このとき,多重比較法にはさまざまなものがある. たとえば, Bon-の発生間隔は, 指数分布に従うことがわかっている. この研T(µ1 ,µ2 ) ≡究では, 2010 年 11 月から 2012 年 6 月までに発生した地震の M5 以上の地震が起きてから次の M5 以上の地震が起µ2 X̄n12n1〜 F2n2µ1 Ȳn2(1)が成り立つ. ただし, X̄n1 と Ȳn2 はそれぞれ, Xi と Yj のきるまでの時間のデータを使って, 解析を行う.

標本平均とする.

2分布論

2.1

指数分布

証明 TX , TY は互いに独立とする. このとき, 参考文献 [1]の定義 5.5 と参考文献 [3] の定理 3.20 より,

ま ず は, 補 題 1 と し て, 知 ら れ て い る 結 果 を 述 べ る .

TX 〜 χ2m , TY 〜 χ2n ⇒補題 1 Z1 , Z2 , · · · , Zn は互いに独立で同一の密度関数

fZ (x) = (1/µ)e(−x/µ) I(0,∞) (x) をもつ指数分布 EX(1/µ)に従うとする. このとき,( )2nT ≡Z̄n 〜 χ22nµn∑が成り立つ. ただし, Z̄n ≡ (1/n)

Zi とする.

証明 Z1 , Z2 , · · · , Zn 〜 EX(1/µ) = GA(1, 1/µ)
ただし, GA(1, 1/µ) は形状母数 1, 尺度母数 1/µ のガンマ分布である.
参考文献 [3] の系 3.16 より,)
1W ≡Zi 〜 GA n,µi=1(n∑( )nfW (x) =T1 ≡1µΓ(n)

2n1となり, (2) で,m ← 2n1 , TX ← T1 , n ← 2n2 , TY ← T2 のように適用
すると,
T1µ2 X̄n12n12n1=〜 F2n2T2µ1 Ȳn22n2となる.

□命題 1 より

(2n1F2n(1 − α/2)2<µ2 X̄n1µ1 Ȳn2)<2n1F2n(α/2)
2= 1−αとなり, 確率の中を書き換えることにより, µ2 /µ1 の信頼xn−1 − µxe

係数 1 − α の信頼区間
2n1Ȳn2 F2n(1 − α/2)2となる.

T = 2W/µ であるので, a = 2/µ, b = 0 として参考文献 [3]の系 2.34 を適用すると,()
1y−bfT (y) = fW|a|a( )n ( )n−11yy1µ=( )e−µ( 2µ )22Γ(n) µµ(2)X̄n1 〜 χ22n1µ12n2T2 ≡Ȳn 〜 χ22n2µ2 2Pとなり, W の密度関数は,〜 Fnm

補題 1 より,□i=1TXmTYnX̄n1<µ2µ1<2n1Ȳn2 F2n(α/2)
2X̄n1を得る.
帰無仮説 H0 : µ1 = µ2 vs. 対立仮説 H1 : µ1 ̸= µ2を考える.
命題 1 より, 帰無仮説 H0 の下で,T ≡X̄n12n1〜 F2n2Ȳn2(3)であるので,( 2n1)
2n1P F2n(1 − α/2) < T < F2n(α/2) = 1 − α22P (帰無仮説 H0i が正しいにもi=1(5)

この多重比較について最大タイプ

2n12n1T > F2n(α/2) または T < F2n(1 − α/2)
22FWER を α 以下にするためには, (5) 式の右辺に現れているそれぞれの項につのときに H0 を棄却し, µ1 ̸= µ2 と判定する.
2n1
ここで, F2n(α/2) は自由度 (2n1 , 2n2 ) の F 分布の上側
2100α% 点を表す.
F 分布

いて,P (帰無仮説 H0i が正しいにもかかわらず棄却される)≤ α/m(i = 1, 2, …, m)が成立すれば,命題 1 の (1) と (3) より, 次の補題 2 を得る.

m∑

補題 2 Xij 〜 EX(1/µi ) (j = 1, · · · , ni , i = 1, · · · , k)

i=1

P (帰無仮説 H0i が正しいにもかかわらず棄却される)≤ m ×（α/m）= αとし, {Xij | j = 1, · · · , ni , i = 1, · · · , k} は互いに独立ととなる. したがって,する. このとき,µ µ̂ii′Tii′ (µ) =µi µ̂i′

P (帰無仮説のうち少なくとも 1 つの H0i が正しいにもかかわらず棄却される確率)≤ α2ni〜 F2ni′となり, 帰無仮説 H(i,i′ ) : µi = µi′ の下で,µ̂iTii′ =µ̂i′

3.12ni〜 F2ni′が成り立つ. ただし, µ̂i ≡ (1/ni )ni∑

Xij (i = 1, · · · , k) とj=1する.

□補題( 2 よりµi′ µ̂i2ni(α/2)< F2ni′µi µ̂i′2ni(1 − α/2) <F2ni′Pm∑

かかわらず棄却される)となる. したがって,2.2かわらず棄却される) ≤)すべての平均相違の場合のボンフェローニ法A′{ 帰無仮説 H(i,i′ ) vs. 対立仮説 H(i,i′ ) |1 ≤ i < i ≤ k}
に対する多重比較検定を補題 2 を用いてボンフェローニ法で行う方法を論述する.
手順 1. 推測の対象となるファミリーを明示し, ファミリーに含まれる帰無仮説の個数 k C2 = {k(k − 1)}/2 個を求める.

= 1−αF = {H(i,i′ ) |1 ≤ i < i′ ≤ k}となる, 確率の中を書き換えることにより, µi′ /µi の信頼手順 2. 有意水準 α を定める. 有意水準は α = 0.05 もしく

係数 1 − α の同時信頼区間は 0.01 と定めることが多い.

2niµ̂i′ F2n(1 − α/2)i′µ̂i<µi′µi<2niµ̂i′ F2n(α/2)
i′

手順 3. ファミリーに含まれるそれぞれの帰無仮説に対して検定統計量 Tii′ を選ぶ.

µ̂iを得る.
したがって,手順 4. データを取り, 検定統計量 Tii′ を計算する.

A
帰無仮説 H(i,i′ ) : µi = µi′ vs. 対立仮説 H(i,i′ ) : µi ̸= µi′に対する水準 α の検定は,Tii′>

手順 5. 各検定統計量 Tii′ について有意水準 α/{k(k−1)/2}2ni(α/2) または Tii′
F2ni′

<2niF2n(1 − α/2)i′に対する検定は,2niTii′ > F2n(α/{k(k − 1)}) またはi′2niTii′ < F2ni′ (1 − α/{k(k − 1)})のときに H(i,i′ ) を棄却し, µi ̸= µi′ と判定する.

3ボンフェローニ法のときに H(i,i′ ) を棄却する.

参考文献 [2] を基にボンフェローニ法を説明する. ボン

H0 : µ1 = · · · = µk の下で,フェローニ法とは, m 個の事象 Di (i = 1, 2, …, m) に対

ボンフェローニの方程式より,して,P(m∪i=1Di ) ≤m∑

P (Di )

(4)

i=1

∪

P0 (∑Aii′ ) ≤1≤i<i′ ≤k∑=(4) の不等式に基づく多重比較法である.
(4) は参考文献 [3] の命題 2.5 で証明されている.

P0 (Aii′ )1≤i<i′ ≤k1≤i<i′ ≤k{{= k C2 ×

推測の対象となるファミリーをF = {H01 , H02 , …, H0m }}αk(k−1)2α}k(k−1)
2=αとする. Di を帰無仮説 H0i が正しいにもかかわらず棄却される事象と考える. これを (4) 式に当てはめると,P (帰無仮説のうち少なくとも 1 つの H0i が正しいにもかとなる. ただし,
2niAii′ ≡ {Tii′ > F2n(α/{k(k − 1)}) またはi′2ni(1 − α/{k(k − 1)})}
Tii′ < F2ni′

表 5. 2011 年 9 月とする. 同様に, 信頼係数 1 − α の同時信頼区間は,2niµ̂i′ F2n(1 − α/{2k(k − 1)})i′<µ̂iµ i′2niµ̂i′ F2n(α/{2k(k − 1)})i′平均分散最小
最大<µi(1 ≤ i < i′ ≤ k)µ̂iで与えられる．

42011 年 9 月

日

時間

1.666
39 : 59
3.963
95 : 070.009
00 : 146.765 162 : 22
標本数 18データ解析

5地震データM5以上
以上

表 1〜表 5 は, 参考文献 [4] のデータをもとに 2010 年 1118月から 2012 年 6 月までの M 5 以上で日本近海で発生した地震の, ある地震が発生してから次の地震が発生するまで161412の時間差を求め, 基本統計量を計算した表である. 東北地
方太平洋沖地震が発生した 2011 年 3 月が 203 回と最も発生回数が多く, 地震が発生してから, 次の地震が発生するまでの時間差の平均は 3 時間 17 分で, 最も短いことがわかった. 標本数の最大値は, 2011 年 3 月の 203 回が最大で, 他

1086M5以上
420の月は, ほとんどの月が 10 回から 20 回程度となっていて,多い月でも 40 回程度である. このことからも, やはり東北地方太平洋沖地震の発生した 2011 年 3 月は異常であった.
表 1. 2010 年 11 月

2010 年 11 月
平均
分散
最小最大日

時間

7.59.598
4.737
11.397
標本数 5180 : 00
230 : 22113 : 42273 : 32

図 1 M 5 以上の地震の発生間隔の平均
図 1 と, 発生回数の表より, 特徴があった 2011 年 3 月

16 日〜31 日, 2011 年 9 月, また東北地方太平洋沖地震の直前の月である 2011 年 2 月 1 日〜15 日, 2010 年 11 月で,3.1 節の手順に添って多重比較による解析を行う.
表 6. データ
群

表 2. 2011 年 2 月 1 日〜15 日

2011 年 2 月
平均
分散
最小最大日

2.54.022
0.611
5.048
標本数 6時間60 : 00
96 : 31
14 : 40
121 : 10

第1群
第2群
第3群
第4群

年月
2010 年 11 月
2011 年 2 月 1 日〜15 日
2011 年 3 月 16 日〜31 日
2011 年 9 月

標本平均

標本数 ni

6.000
2.500
0.280
1.666

5657
18

手順 1. ファミリーの個数は, 4 C2 = 6 個

F = {H(1,2) , H(1,3) , H(1,4) , H(2,3) , H(2,4) , H(3,4) }表 3. 2011 年 3 月 16 日〜31 日

2011 年 3 月
平均
分散
最小最大日

0.280
0.1180.004
1.516
標本数 57

時間06 : 4302 : 5000 : 06
36 : 24

表 4. 2011 年 4 月

2011 年 4 月
平均分散
最小
最大日0.6382
0.54030.00073.493
標本数 47

時間

15 : 19
12 : 58
00 : 01
83 : 51

手順 2. 有意水準を α = 0.05 とする.
手順 3. 検定統計量をTii′ =µ̂iµ̂i′

2ni〜 F2ni′とする.
手順 4. データを取り, 検定統計量 Tii′ を計算する.

(1)T12 = 6.000
2.500 = 2.400(1) と同様に計算すると, (2) T13 = 21.429, (3) T14 =3.601 , (4) T23 = 8.928, (5) T24 = 1.501, (6) T34 = 0.168となる.
手順 5. µ2 /µ1 の信頼係数 1 − α の同時信頼区間は,0.070 <µ2µ1< 2.212である. 検定統計量 T12 に対する検定は,2n1T12 > F2n(α/12) = 5.308 または22n1(1 − α/12) = 0.169
T12 < F2n2のときに H(1,2) を棄却する.

3. 有意水準 α と 1 − α を計算する.

れない. H(1,3) , H(1,4) , H(2,3) , H(2,4) , H(3,4) についても同

4. kentei 関数により, 同時信頼区間, 棄却域を計算出力
し, 第 i 群と第 j 群が帰無仮説 H0 : µi = µj の下で棄却されるか, 採択されるかを出力する.

様に計算を行っていくと, 以下の結果を得る.

5. kekka 関数により, 検定の結果を表にして出力する.

よって, 0.169 < T12 = 2.400 < 5.308 となるので, 棄却さ

表 7. データ解析結果
H(i,i′ )H(1,2)
H(1,3)
H(1,4)
H(2,3)
H(2,4)
H(3,4)

5.1

検定の結果棄却されない
棄却される
棄却される
棄却される
棄却されない棄却される

信頼区間の結果
大小関係はわからないµ3 < µ1µ4 < µ1µ3 < µ2
大小関係はわからないµ4 > µ3注) このプログラムでは C++ 言語の特徴であるオブジェ
クト指向は用いていないが, 情報システム数理学科の実習の際, Windows に C++ 言語のコンパイラをインストールしたため C 言語ではなく, C++ 言語を用いてプログラムの作成をした. なお, 本論で述べた変更方法により, C コンパイラでも実行可能形式のプログラムを作成できる.

6.1

データ解析結果

実行結果

以下はプログラムの実行結果の一部である.

検定の結果より, H(1,2) と H(2,4) が棄却されなかった.

誤差 0.000010 以下の自由度 ( 8 , 12 ) の F 分布の上側

このことから, 1 群と 3 群, 1 群と 4 群, 2 群と 3 群, 3 群と100 α パーセント点は 5.583336,
1 − 100 α パーセント点は 0.134972µ2 /µ1 の信頼係数 1-αの同時信頼区間は,0.044991 < µ2 /µ1 < 1.861112
検定統計量 T(1,2) に対する検定は,T(1,2) > F(8,12;0.004167) = 5.583336
または T(1,2) < F(8,12;0.995833) = 0.134972のときに H(1,2) を棄却する.
0.134972 < T(1,2) = 3.000000 < 5.583336 となるので,棄却されない.のように出力される.

4 群の間で違いがあることがわかった.
信頼区間の結果より,µ3 < µ4 < µ1 , µ3 < µ2⇐⇒
地震の平均発生時間間隔に関して,
2011 年 3 月 16 日〜31 日 <2011 年 9 月 <2010 年 11 月,
2011 年 3 月 16 日〜31 日 <2011 年 2 月 1 日〜15 日となる.

6

C++ 言語によるボンフェローニ法での検定
プログラム
以下のプログラムは C++ げんごによるボンフェローニおわりに7法での検定を行うプログラムである. ページ数の都合上, C

この研究では, 指数分布を用いているため, 地震の発生時言語によるボンフェローニ法での検定プログラムは, main間間隔でボンフェローニ法を用いて検定を行った. 発生時プログラムだけとし, その説明を行う.

間間隔が短いということは, 地震の発生回数は多いという

double alpha1,alpha2;
nyuuryoku(); /*群の個数と有意水準αを入力する関数*/suuchi(); /*標本平均と標本数を入力する関数*/
printf(”/n”);
toukeiryou(avr); /*検定統計量を計算する関数*/
printf(”/n”);
alpha1=alpha/(k*(k-1)); /*有意水準 α/k(k-1) を計算する式*/
alpha2=(double)1.0-alpha1; /*有意水準 1-α/k(k-1) を計算する式*/
kentei(alpha1,alpha2,k); /*信頼区間, 検定統計量に対する検定を多重比較する関数*/
kekka(avr,k); /*検定の結果を出力する関数*/
return 0;
プログラムの流れは,ことになる. よって, 2011 年の 3 月が最も地震の発生回数

1. nyuuryoku 関数と suuchi 関数によって, 扱うデータの群の個数, 有意水準, 標本平均, 標本数 (ni ) を入力する.
2. toukeiryou 関数で検定統計量を計算する.

が多いことが分かった. また, 2011 年の 9 月の方が 2010年の 11 月よりも多く地震が発生していることが分かった.
さらに, C++ というプログラム言語を用いて, データ解析を行うプログラムを作成した. このプログラムは地震のデータだけでなく指数分布に従う他のデータにも適用することができる.

参考文献
[1]

白旗慎吾:『統計学』.
ミネルヴァ書房, 京都, 2008.

[2]

永田靖・吉田道弘:『統計的多重比較法の基礎』.
株式会社

[3]

サイエンティスト社, 東京, 1997.

白石高章：『統計科学の基礎 データと確率の結びつきがよくわかる数理』.
日本評論社, 東京, 2012.

[4]『日本気象協会』 .

http://www.jwa.or.jp/