坂元 慶行, 石黒 真木夫, 北川 源四郎:情報量統計学

作成日 : 2017-11-03
最終更新日 :

概要

情報量規準で統計的手法を統一的にまとめようとする本書は,入門解説書であると共に,問題提起書でもある.

感想

日本生まれの統計手法

日本で生まれ育った(赤池)情報量規準(AIC)を用いた統計手法の解説。 基本となる検定問題について新たな考え方が述べられている。統計表がいらない、という 利点もあり、もっとこの方法は知られてよいのではないかと思う。

AIC は特に回帰分析で威力を発揮していて、多変量解析や自己回帰モデルの変数選択では統計処理のパッケージに必ずといっていいほど入っている。 本書では回帰分析のみならず、分割表モデルや分散分析でも AIC が適用可能であることを示していて、本書の価値を高めている。

変数変換

小針:確率・統計入門で誤っていた変数変換の式が、 本書ではどのようになっているか確認したかった。 p.10 における 1・6 変数変換では、正しい式になっている。以下本書から引用する。

`f(x)` を確率変数 `X` の密度関数,`g(x)` を `x` の単調増加関数とする. `Y = g(x)` で新しい確率変数 `Y` を定義すると,`Y` の分布関数は,

`F(y)`` = P(Y le y) = P(X le g^(-1)(y))`
` = int_(-oo)^(g^-1(y)) f(x) dx`

で計算される.ここで `y = g(x)` なる変数変換を行うと,
`(dy)/(dx) = g'(x) = g'(g^-1(y))`
を用いて,
`F(y) = int_-oo^y f(g^-1(y)) / (g'(g^-1(y))) dy`
となる.これは,`Y` の密度関数 `h(y)` が,
`h(y) = f(g^-1(y)) / (g'(g^-1(y))) `
で与えられることを示している.`g` が単調減少関数である場合にも同様に考えて,`g` が単調関数であれば `Y` の密度関数が,
`h(y) = f(g^-1(y)) / abs(g'(g^-1(y))) `
で与えられることが導かれる.

分散分析について

第9章は分散分析である。p.155 では次のように始められている。

(前略)この章では,目的変数が連続的,説明変数が離散的なデータを解析するためのモデルとして,分散分析モデルと名づけるモデルを導入する.

このモデルにおいても,目的変数の動きを説明する必要最小限の説明変数の組を見出すことが重要な課題となる.(中略) 分散分析は,目的変数の標本分散を各説明変数に起因する部分とそれ以外の未知の要因による部分の和に分解し,ある説明変数からの寄与が‘小さ’ければこの変数からの寄与はなかったとみなす方法である. このとき,‘小さい’か否かの判定基準を客観的に設定できないのが重大な欠点であった.

その後、AIC による分散分析の理論と計算式が説明される。「9.1 分散分析モデル B. パラメータの最尤推定量」の p.159 で述べられているのは次のようなことである。 すなわち、特定のパラメータ `bb theta` によって構成されたモデル `"MODEL"(1, 2, cdots, m)` の条件付き対数尤度 `l(bbtheta)` を適切な制約の下で最大化する `bbtheta` は、最尤推定量 `hat bbtheta` であることが述べられる。ここで次の一文がある。

一般の場合は計算が複雑になるのでここでは次の 2 つの仮定を置く.(後略)

この一般の場合の計算については、本書の本文では述べられていない。第Ⅲ編のプログラムを読み解いてその計算を理解しようとしたのだが、残念ならプログラムは動作しなかった。 ではどうすればいいのか。 p.166 で「B. § 9.1B の仮定が成り立たない場合」という項があり、次のようなことが書かれている。

この場合は(9.39)~(9.47)を用いて解を得ることはできないが `"MODEL"(j_1, cdots, j_r)` を仮定したときの `sigma^2` の推定値 `hat sigma^2(j_1, cdots, j_r)` を用いて,

`"AIC"(j_1, cdots, j_r) = n log2 pi + n log hat sigma^2(j_1, cdots, j_r) + n + 2{sum_(j in J) (c_j - 1) + 2}`
が計算できるのは同じである.

そうか、それならば `hat sigma^2(j_1, cdots, j_r)` が計算できればいいわけだ。そこで他書を読みながら無い知恵を絞って考えた結果、やっとたどり着いたのが以下に記すことである。

目的変数 `Y` の値とともに、説明変数として 2 個の特性(因子)`X_1, X_2` が測定されている場合を考える。特性 `X_1` は 2 個の, `X_2` は 3 個のカテゴリー(水準)からなるとする。 `(Y, X_1, X_2)` の測定値 `(y_i, x_(1i), x_(2i)), i = 1, cdots, 9` が与えられているとする。ここで、特性 `X_1` だけを考えると、 どのようにモデルを構成すればいいだろうか。これは次のように考える。まずカテゴリーにかかわらず一定の定数 `theta_0` があるとする。 次に、 2 個のカテゴリーを `[1, 2]` としたときに、`X_1` がカテゴリー 1 に属するときにはある値 `theta_1` が、 `X_1` がカテゴリー 2 に属するときには別の値 `theta_2` が、モデルを構築するために使われると考えればいい。つまり、

`y_i ~= {(theta_0 + theta_1 (x_(1i) = 1)),(theta_0 + theta_2 (x_(1i) = 2)) :}`
である。さらに、`theta_0 + theta_1 = 0` であるから、
`y_i ~= {(theta_0 + theta_1 (x_(1i) = 1)),(theta_0 - theta_1 (x_(1i) = 2)) :}`
となる。

今度は、特性 `X_2` だけを考える。同様に考えると、 3 個のカテゴリーを `[1, 2]` としたときに、`X_2` がカテゴリー 1 に属するときにはある値 `theta_3` が、 `X_2` がカテゴリー 2 に属するときには別の値 `theta_4` が、`X_2` がカテゴリー 3 に属するときにはまた別の値 `theta_5` が、モデルを構築するために使われると考えればいい。つまり、

`y_i ~= {(theta_0 + theta_3 (x_(2i) = 1)),(theta_0 + theta_4 (x_(2i) = 2)), (theta_0 + theta_5 (x_(2i) = 3)) :}`
である。さらに、`theta_3 + theta_4 + theta_5 = 0` であるから、
`y_i ~= {(theta_0 + theta_3 (x_(2i) = 1)),(theta_0 + theta_4 (x_(2i) = 2)),(theta_0 - theta_3 - theta_4(x_(2i) = 3)) :}`
となる。

さて、特性 `X_1` だけを考えたとき、パラメータ `theta_0` とパラメータ `theta_1` はどのように求められるか。これには、次の近似式を用意する。 仮に、`x_(1i)` が、`i` が奇数のとき1、`i` が偶数の時 -1 `であったとすると、次の式が成り立つ。

`[[y_1],[y_2],[vdots],[y_10]] ~= [[1,1],[1,-1],[vdots,vdots],[1,-1]][[theta_0],[theta_1]]`
この `~=` の意味は、次の値 `L` を最小化するという意味である。ここで右辺の行列を `A = (a_(ij))` とおく。この `A` はデザイン行列と呼ばれる。
`L = sum_(i = 1)^10 (y_i - (a_(i1)theta_0 + a_(i2)theta_1))^2`

この `L` を最小化するために `theta_0, theta_1` を求めるために行なう計算として、一つは正規方程式を作って `(theta_0, theta_1)` を求めてから `L` を求めるもので、 `L/n` が `hat sigma^2` を与える。もう一つの方法はデザイン行列 `A` の右に左辺の `y` ベクトルを接合した拡大行列 `X`

`X = [A | y] = [[1,1, |, y_1],[1,-1, |, y_2],[vdots,vdots, |, vdots],[1,-1, |, y_(10)]]`
を作り、この `X` をハウスホルダー変換して得られる3次の上三角行列 `R`
`R = [[S_(11), S_(12), S_(13)],[0, S_(22), S_(23)], [0, 0, S_(33)]]`
から `L = S_(33)^2 / n` として、また `theta_0, theta_1` は次の連立方程式を解いて求められる。
`[[S_(11), S_(12)],[0,S_(22)]][[theta_0],[theta_1]] = [[S_(13)],[S_(23)]]`
詳細は本書の pp.149-153 などを参照のこと。

特性 `X_2` だけを考えた場合のデザイン行列 `A` は次のようになる。ただし、`x_(2i)` は、`i` を 3 で割った時の剰余(ただし剰余が 0 のときは 3)であったとする。

`A = [[1,1,0],[1,0,1],[1,-1,-1],[1,1,0],[vdots,vdots,vdots],[1,1,0]]`

推定すべきパラメータは `sigma^2` のほかは、`theta_0, theta_1,theta_2` である。あとも同様である。

特性 `X_1` と `X_2` を考えた場合のデザイン行列 `A` は次のようになる。ただし、`x_(1i), x_(2i)` は、上記と同様とする。

`A = [[1,1,1,0],[1,-1,0,1],[1,1,-1,-1],[1,-1,1,0],[vdots,vdots,vdots, vdots],[1,-1,1,0]]`

推定すべきパラメータは `sigma^2` のほかは、`theta_0, theta_1,theta_2, theta_3` である。あとも同様である。

(この項 2025-01-27

第Ⅲ編のプログラムについて

p.171 以降の第Ⅲ編では、FORTRAN 77 のプログラムが掲載されている。プログラムがあるのはいいのだが、FORTRAN 77 であるので、稼働させるのは一苦労である。 たとえば、Ⅲ-4 分散分析プログラム VARMOD についてみてみよう。本プログラムを見る限り、下記の点に注意しないといけない。なお、 以下のプログラムは GNU Fortran (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0 で確認している。

CONTINUE の位置

DO ループの末尾には CONTINUE 文を置くのが普通であるが、行数を節約したいがためか、CONTINUE 文をはしょってしまっている。これは紛らわしい。仮に次のようなプログラムがあったとする。


      PROGRAM MAIN
      IMPLICIT NONE
      INTEGER i
      DO 10 i=1,3
      WRITE(*,*)i
   10 WRITE(*,*)'Saluton!'
      END

このプログラムをコンパイルして得られる結果は

					1
	Saluton!
					2
	Saluton!
					3
	Saluton!

であって

					1
					2
					3
	Saluton!

ではない。gfortran は、文番号がある行に continue がないと次の警告を出す。

Warning: Fortran 2018 deleted feature: DO termination statement which is not END DO or CONTINUE with label 10 at (1)

書き換えるには次のようにする。


      PROGRAM MAIN
      IMPLICIT NONE
      INTEGER i
      DO 10 i=1,3
      WRITE(*,*)i
      WRITE(*,*)'Saluton!'
   10 CONTINUE
      END

さらにいいのは、文番号を使わず ENDDO を使う方法だ。


      PROGRAM MAIN
      IMPLICIT NONE
      INTEGER i
      DO 10 i=1,3
      WRITE(*,*)i
      WRITE(*,*)'Saluton!'
      ENDDO

なぞの 1H

p.209 の行番号 00001750 以降のFORMAT 文の最初に 1H1 とか 1H0 のような記号があり、出ていてよくわからないまま写経をしていた。コンパイルすると表示に 1 や 0 が出てくるのでおかしい。 この H はホレリスの H ではないらしい。調べてみると、どうやら行送りの編集記述子だとわかった。 こういう指示はすべて削除することにした。

1H1次のページの最初の行に送る
1H02 行送る
1H_(スペース)1 行送る
1H+行送りしない

参照:アルゴリズムの基礎(FORTRAN演習) ファイル入出力(coastal.nagaokaut.ac.jp)

配列の読み込みと書き出し

p.207 の行番号 00000460 から 00000460 にかけて次のプログラムがある。

      READ( 5,100 )     NFCTR
      READ( 5,100 )     (NLEVEL(I),I=1,NFCTR)

最初はいい。NFCTR という変数にファイルから読み込んだ値を読み込ませればいい。次の行が問題だ。NLEVEL という配列の変数に i = 1, ..., NFCTR まで読み込ませたい、 という意図はわかる。しかし、実際には読み込まれない。次の行の値を読んでしまう。これを解決するには部分配列を使ってこうすればいい。

      READ( 5,100 )     NFCTR
      READ( 5,100 )     NLEVEL(1:NFCTR)

誤植

p.134 の d(1) は 0.01333 という値になっているが、p.135 の表 8.2 における次数 1 の残差分散は 0.01329 という値になっていて、両者が異なる。 私の計算では 0.013312... と続く。

p.164 の `hat theta_1(3)` の計算結果が 79.275 となっているが、正しくは 7.9375 である。

p.212 の天(上部)、ノンブルがある行の中央、Ⅲ-4 分散分析プグロラムとなっているが、正しくは《Ⅲ-4 分散分析プログラム》である。

書誌情報

書名情報量統計学
著者坂元 慶行, 石黒 真木夫, 北川 源四郎
発行日1989 年 6 月10日 (初版5刷)
発行元共立出版
定価
サイズ
ISBN4-320-02171-1
NDC417
その他

まりんきょ学問所統計活用術統計の本 > 坂元 慶行, 石黒 真木夫, 北川 源四郎:情報量統計学


MARUYAMA Satosi