C 言語と Mathematica による統計プログラミングの基礎
2009SE046 早川由宏
指導教員：白石高章
し同じ Ij (j = 1, · · · , J) に含まれる添字をもつ母平均ははじめに

1

本研究では, 統計解析を行うために有効な C 言語のプログラムの開発を行った. 多重比較統計量の上側 100α% 点を求めるために二重積分の方程式を解く必要がある. その等しいという帰無仮説を H(I1 , · · · , IJ ) で表す. このとき,

ϕ ⊊ V ⊂ U を満たす任意の V に対して, ある自然数 J と上記のある I1 , · · · , IJ が存在して∧実現のため, 数式処理ソフトウェアの Mathematica を使う手順を示した.

C 言語での計算プログラム

2

Hv = H(I1 , · · · , IJ )v∈V

統計学で必要とされる基礎的なプログラムに加え検定を行うプログラムを C 言語で作成した. 基礎的なプログラム例は以下のとおりである. 多くは文献 [3] のアルゴリズムが成り立つ. T (Ij ) ≡maxi<i′ ,i,i′ ∈Ijk ni1 ∑∑X̄i· − X̄i′ ·, VE ≡(Xij − X̄i· )Tii′ ≡ √n − k i=1 j=1
VE ( n1i + n1′ )i∧
とおき水準 α の帰無仮説
Hv に対する検定を考えるこを参考にした.

v∈V「標準正規分布の密度関数の値を求める」「標本平均, 標本分散, 標本標準偏差, 標本歪度, 標本尖度を求める」「2 次元データを入力し, 各成分の標本平均, 標本分散, 標本相関係数から最良線形単回帰直線, 残差平方和と寄与率を「n 個 の 観 測 値 を 入 力 し, 順 序 統 計 量 を 求 め, 標 本

5, 10, 25, 50, 75, 90, 95 パーセント点および標本範囲, 最大値, 最小値を計算する」「t, χ2 , F, t 分布の上側確率の値を求める」「Mersenne Twister(MT19937) により発生した擬似一様乱数から標準正規分布, 混合正規分布, 異常値をもつ混合正規分布, 対数正規分布, ロジスティック分布, 両側指数分
布, 指数分布に従う乱数を生成する」
MT19937 は 219937 − 1 という周期を持つ擬似乱数生成のソフトウェア (文献 [4]) で上記プログラムはその出力を利用している.
閉検定手順のプログラム

Xij ∼ N (µi , σ ) (j = 1, · · · , ni ; i = 1, · · · , k) とし,すべての Xij は互いに独立とする. k を群の個数, ni は第k∑i 群のサイズで n =ni とする.
i=1

帰無仮説 H(i,i′ ) : µi = µi′ vs.

µ(1) J ≥ 2 のとき l = l1 , · · · , lJ に対して α(M, l) ≡ 1 −(1 − α)l/M で α(M, l) を定義する. 1 ≤ j ≤ J となるある整数 j が存在して ta(lj , m; α(M, lj )) < T (Ij ) ならば帰無∧仮説
Hv を棄却する.

′(2)∧ J = 1 のとき ta(M, m; α) < T (I1 ) ならば帰無仮説
Hv を棄却する.
v∈V(1),(2) の方法で, (i, i′ ) ∈ V ⊂ U を満たす任意の V に∧対して,Hv が棄却されるとき, 多重比較検定として
v∈V

H(i,i′ ) を棄却する. 本研究で作成したプログラムは k = 5以下, 第 i 群の標本数は 100 以下でのみ動作する. 今回検定するデータは n = 75, k = 5 となるので以降はそのデータを解析するための方法を述べる.

Mathematica による計算

3検定を行うために ta(lj , m; α(M, lj )) を求める.

3.12i′とができる.

v∈V求める」

2.1|Tii′ | (j = 1, · · · , J),

A
対立仮説 H(i,i′ ) : µi ̸=分布関数の定義

関数 T A(t|lj ) を∫ ∞T A(t|lj )≡lj
A(ts|lj ) ≡(1 ≤ i < i ≤ k) に対する文献 [1] で提案された閉検定手順のプログラム (文献 [2]p.122 検出力の高い閉検定
手順 I) を C 言語により作成した. そのアルゴリズムを述}lj −1dΦ(x)T A(t|lj ) ≡ 1 − α(M, lj )を満たす解とする. ただし m≡n − k, Φ(x) は N (0, 1) の分布関数とし,g(s)≡v∈Vつかが等しいという仮説となる. I1 , · · · , IJ (Ij ̸= ϕ, j =1, · · · , J) を添え字 1, · · · , k の互いに素な部分集合の組と−∞√Φ(x) − Φ(x − 2·ts)と定義し, ta(lj , m; α(M, lj )) は方程式

べる.

{}U ≡ (i, i′ )|1 ≤ i < i′ ≤ k とおくと帰無仮説のファ{}
ミリー H ≡ Hv |v ∈ U となる. さらに ϕ ⊊ V ⊂ U を∧満たす V に対して,Hv は k 個の母平均に関していく

A(ts|lj )g(s)ds,0∫ ∞{とする.

2mm/2
sm−1 e−ms /2
m/2−1
Γ(m/2)2

3.2

Mathematica での定義

Mathematica で T A(t|lj ) は文献 [5] を参考に「:=」や「^」,「*」などの記号を使って定義をする. 式は数値積分を用いて表現すると
f[x_] = Exp[-x^2/2]/Sqrt[2 Pi]
F[x_] = Integrate[f[z], {z, -Infinity, x}]A[t_] := NIntegrate[f[x] (F[x] - F[x Sqrt[2]*t])^(l - 1), {x, -5, 5}]
G[a_] := Integrate[x^(a - 1) Exp[-x],{x, 0, Infinity}]g[s_] := m^(m/2)/(G[m/2] 2^(m/2 - 1))*s^(m - 1)*Exp[-m*s^2/2]
TA[t_] := l*NIntegrate[A[t*s]*g[s], {s, 0, 5},
Method -> DoubleExponential]
Off[NIntegrate::inumr]h[x_] := TA[x] - 1 + al
al := 1 - (1 - a)^(l/M)となる.

3.5

計算結果

Mathematica で求めた ta(l, m; α(M, l)) の一例を以下に示す.
表 1 α = 0.01 のときの ta(l, 70; α(M, l)) の値
M \l5
432

42
2.976
2.898
*
2.648

33.187
*
3.011

4*
3.229

53.384

C 言語プログラムによる閉検定手順のアルゴリズム
プログラムは main 関数の中に検定をする関数がある. main 関数は変数宣言で始まり Mathematica で得た「Method -> DoubleExponential」は二重指数求積法解をメモリーに格納する. 標本がメモリーに格納されるとを使う指定コマンドで無指定時と比べ計算が高速化する.

統計量 Tii′ などを計算し, ファミリーを 4 次元の配列に格「Off[NIntegrate::inumr]」は警告を非表示にするコマ納する. ファミリーの作成は検定の関数と独立しており kンドで計算を継続させる. これらの式は複雑で方程式を解の値によって場合分けして出力する. 今回使う k = 5 はく専用のコマンド Solve が使えないので二分法を用いる.

る性質を使う. h(x) = 0 を解く際, 区間 [p, q] 内に解がk = 2, 3, 4 と同様のアルゴリズムで作成後足りない帰無仮説を追加する. 検定する関数は得られたファミリーに含まれる i, i′ について |Tii′ | の最大値を求め添字を控える. 最後に Tii′ と ta(lj , m; α(M, lj )) を比較し Hv が棄却されない場合に検定の関数を終了する. この検定の関数を二重の繰り返しで検定を繰り返し行う.

ある時 h(p)h(q) ≤ 0 が常に成り立つのでこの区間を狭く

53.3

二分法のアルゴリズム

この研究では, 方程式の近似解を求める方法として, 二分法を使用する. ただし, 分布関数は連続で単調増加であすれば良い. h((p + q)/2) < 0 ⇒ 区間を [(p + q)/2, q],h((p + q)/2) ≥ 0 ⇒ 区間を [p, (p + q)/2] と新しく定めq − p ≤ 10−5 程度まで繰り返す. Mathematica の繰り返
しコマンド For を用いて最初の区間の決定は
For[p = 1, ! (h[p] < 0 && 0 < h[p + 1]), p++]で区間を縮めるにはFor[q = p + 1, Abs[p - q] > 10^-5,If[h[(p + q)/2] < 0, p = (p + q)/2, q = (p + q)/2]]

である.

3.4

解を求める

Mathematica で１つだけ解を求めるには関数を定義後単純に二分法を行う. 複数の解を求めさせるには For 文をうまく使えば良い. 解 ta(lj , m, α(M, lj )) は l, m, α, M の4 元で表されている. そこで α, m を固定し M, l を動かすことで複数解を得る. コマンドは
a = 0.01; m = 70
For[M = 5, M != 1,For[l = 2, l != 6, If[l == M - 1, l++];For[For[p = 1, !(h[p] < 0 && 0 < h[p + 1]),p++];s = p; t = p + 1, Abs[s - t] > 10^-5,If[h[(s + t)/2] < 0,s = (s + t)/2, t = (s + t)/2]];
Print["M=", M, ",l=", l "のとき t=", N[s, 4]];l++; If[M < l, Break[]]]; M--]とすれば 2≤M ≤5, 2≤l≤5 で必要な解 7 個を得られる.

おわりにC 言語では 2 重積分など複雑な計算は簡単にできない.
C 言語の文法, 数学の近似法などの知識に加えオーバーフローなどコンピュータ特有のエラーにも対処する必要があり, 精度保証は難しいが情報の整理がしやすい. 一方
Mathematica は C 言語のように複雑な手続きはさせにく
いが精度の良い近似値が容易に得られる. それぞれメリット, デメリットを考え適切なソフトを使うことが望ましいことがわかった.

参考文献[1] 白石高章：『多群モデルにおけるすべての平均相違に関する閉検定手順』．計量生物学，2011.
[2] 白石高章：『多群連続モデルにおける多重比較法』．
共立出版，東京，2011.
[3] 白石高章：『FORTRAN による統計プログラミング』．
2011/6/30, 講義ノート.
[4] Mersenne Twister の Web Page
http://www.math.sci.hiroshima-u.ac.jp/~mmat/MT/mt.html[5] Wolfram Mathematica のドキュメントセンター
http://reference.wolfram.com/mathematica/tutorial/
VirtualBookOverview.html