統計計算

作成日 : 2022-04-07
最終更新日 :

統計計算については、統計活用法というページで、はるか昔にまとめたことがあった。新たにわかったことを記す。

カイ二乗分布

カイ二乗分布の応用例に検定がある。 たとえば、クロス表（分割表）から複数項目の独立性を検定するのに使われる。カイ二乗分布の表は統計の教科書の付録のページにあって、 たとえばこんな表になっている。

カイ二乗分布の表
α
`n` 0.99 0.98 0.95 0.90 0.80 0.70 0.50 0.30 0.20 0.10 0.05 0.02 0.01 0.0050.001
1 0.00 0.00 0.00 0.02 0.06 0.15 0.45 1.07 1.64 2.71 3.84 5.41 6.63 7.8810.83
2 0.02 0.04 0.10 0.21 0.45 0.71 1.39 2.41 3.22 4.61 5.99 7.82 9.21 10.6013.82
3 0.11 0.18 0.35 0.58 1.01 1.42 2.37 3.66 4.64 6.25 7.81 9.84 11.34 12.8416.27
4 0.30 0.43 0.71 1.06 1.65 2.19 3.36 4.88 5.99 7.78 9.49 11.67 13.28 14.8618.47
5 0.55 0.75 1.15 1.61 2.34 3.00 4.35 6.06 7.29 9.24 11.07 13.39 15.09 16.7520.52
6 0.87 1.13 1.64 2.20 3.07 3.83 5.35 7.23 8.56 10.64 12.59 15.03 16.81 18.5522.46
7 1.24 1.56 2.17 2.83 3.82 4.67 6.35 8.38 9.80 12.02 14.07 16.62 18.48 20.2824.32
8 1.65 2.03 2.73 3.49 4.59 5.53 7.34 9.52 11.03 13.36 15.51 18.17 20.09 21.9626.12
9 2.09 2.53 3.33 4.17 5.38 6.39 8.34 10.66 12.24 14.68 16.92 19.68 21.67 23.5927.88
10 2.56 3.06 3.94 4.87 6.18 7.27 9.34 11.78 13.44 15.99 18.31 21.16 23.21 25.1929.59
11 3.05 3.61 4.57 5.58 6.99 8.15 10.34 12.90 14.63 17.28 19.68 22.62 24.72 26.7631.26
12 3.57 4.18 5.23 6.30 7.81 9.03 11.34 14.01 15.81 18.55 21.03 24.05 26.22 28.3032.91
13 4.11 4.77 5.89 7.04 8.63 9.93 12.34 15.12 16.98 19.81 22.36 25.47 27.69 29.8234.53
14 4.66 5.37 6.57 7.79 9.47 10.82 13.34 16.22 18.15 21.06 23.68 26.87 29.14 31.3236.12
15 5.23 5.98 7.26 8.55 10.31 11.72 14.34 17.32 19.31 22.31 25.00 28.26 30.58 32.8037.70
16 5.81 6.61 7.96 9.31 11.15 12.62 15.34 18.42 20.47 23.54 26.30 29.63 32.00 34.2739.25
17 6.41 7.26 8.67 10.09 12.00 13.53 16.34 19.51 21.61 24.77 27.59 31.00 33.41 35.7240.79
18 7.01 7.91 9.39 10.86 12.86 14.44 17.34 20.60 22.76 25.99 28.87 32.35 34.81 37.1642.31
19 7.63 8.57 10.12 11.65 13.72 15.35 18.34 21.69 23.90 27.20 30.14 33.69 36.19 38.5843.82
20 8.26 9.24 10.85 12.44 14.58 16.27 19.34 22.77 25.04 28.41 31.41 35.02 37.57 40.0045.31
21 8.90 9.91 11.59 13.24 15.44 17.18 20.34 23.86 26.17 29.62 32.67 36.34 38.93 41.4046.80
22 9.54 10.60 12.34 14.04 16.31 18.10 21.34 24.94 27.30 30.81 33.92 37.66 40.29 42.8048.27
23 10.20 11.29 13.09 14.85 17.19 19.02 22.34 26.02 28.43 32.01 35.17 38.97 41.64 44.1849.73
24 10.86 11.99 13.85 15.66 18.06 19.94 23.34 27.10 29.55 33.20 36.42 40.27 42.98 45.5651.18
25 11.52 12.70 14.61 16.47 18.94 20.87 24.34 28.17 30.68 34.38 37.65 41.57 44.31 46.9352.62
26 12.20 13.41 15.38 17.29 19.82 21.79 25.34 29.25 31.79 35.56 38.89 42.86 45.64 48.2954.05
27 12.88 14.13 16.15 18.11 20.70 22.72 26.34 30.32 32.91 36.74 40.11 44.14 46.96 49.6455.48
28 13.56 14.85 16.93 18.94 21.59 23.65 27.34 31.39 34.03 37.92 41.34 45.42 48.28 50.9956.89
29 14.26 15.57 17.71 19.77 22.48 24.58 28.34 32.46 35.14 39.09 42.56 46.69 49.59 52.3458.30
30 14.95 16.31 18.49 20.60 23.36 25.51 29.34 33.53 36.25 40.26 43.77 47.96 50.89 53.6759.70
40 22.16 23.84 26.51 29.05 32.34 34.87 39.34 44.16 47.27 51.81 55.76 60.44 63.69 66.7773.40
50 29.71 31.66 34.76 37.69 41.45 44.31 49.33 54.72 58.16 63.17 67.50 72.61 76.15 79.4986.66
60 37.48 39.70 43.19 46.46 50.64 53.81 59.33 65.23 68.97 74.40 79.08 84.58 88.38 91.9599.61
70 45.44 47.89 51.74 55.33 59.90 63.35 69.33 75.69 79.71 85.53 90.53 96.39 100.43 104.22112.32
80 53.54 56.21 60.39 64.28 69.21 72.92 79.33 86.12 90.41 96.58 101.88 108.07 112.33 116.32124.84
90 61.75 64.63 69.13 73.29 78.56 82.51 89.33 96.52 101.05 107.57 113.15 119.65 124.12 128.30137.21
100 70.06 73.14 77.93 82.36 87.95 92.13 99.33 106.91 111.67 118.50 124.34 131.14 135.81 140.17149.45

出典は青木先生の下記ページ。
http://aoki2.si.gunma-u.ac.jp/CGI-BIN/txxp.html
ただし`alpha=0.005` の値は薩摩先生の「確率・統計」p.214 より転載（ただし n = 21 から n= 29 を除く）。 `alpha=0.005, n = 21 ～ 29 ` の場合は、青木先生の下記ページにより計算した値を使った。
http://aoki2.si.gunma-u.ac.jp/CGI-BIN/pxx.html

自由度 `n` のカイ２乗分布（`chi^2` 分布）は次の式で与えられる。ここで `Gamma(x)` はガンマ関数である。

`f_n(x) = 1/(2^(n//2)Gamma(n//2)) x^(n//2-1) e^(-x//2)`

カイ2乗分布のグラフは次のようになる。

`n` を固定したとき、塗りつぶしの部分の面積が `alpha` になるときの `x = t` が上記の「カイ二乗分布の表」である。 上記のグラフの例では、`n=1, x = 3.84` であるので、塗りつぶしの部分の面積は表を見て、`alpha = 0.05` だとわかる。 言い換えれば、上の表の数値の値を `t` とすると、`f_n(x)` と `n, t, alpha` の間には次の関係式がある。

`int_t^oo f_n(x)dx = alpha`

`t` と `n` を与えて `alpha` を求める計算を JavaScript で実装し、計算ができるようにした。 `t = 3.84, n = 1` の場合は、`alpha = 0.05004352 cdots` と計算されるはずだ。
`t` : `n` : `alpha` :

`t` 分布

統計的検定として重要度が高いのは、カイ二乗分布による検定の他、`t` 分布による検定がある。 たとえば、2 つの群の平均値に差があるかどうかの検定に使われる。 平均値の差の検定やWelch の t 検定を参照してほしい。 `t` 分布の表も統計の教科書の付録のページにあることが多い。表の一部は次の通り。

`t` 分布の表

出典は稲垣宣生・山根芳知・吉田光雄:統計学入門「付表 2 ティー分布表」。 なお、青木先生の下記ページも参照した。
http://aoki2.si.gunma-u.ac.jp/CGI-BIN/ttxp.html

自由度 `nu` の `t` 分布は次の式で与えられる。ここで `Gamma(x)` はガンマ関数である。

`f_nu(x) = (Gamma((nu+1)//2)) /(sqrt(pi nu)Gamma(nu//2)) (1+x^2/nu)^(-(nu+1)//2)`

通常 `nu` は自然数であるが、Welch の `t` 検定などでは `nu` を正の実数に拡張して使う。 正の実数に対して t 分布を計算する方法として、奥村晴彦「C 言語による標準アルゴリズム事典」に掲載されたアルゴリズムを使った。

`t` 分布のグラフは次のようになる。ν = 1 : 橙、ν = 2 : 赤、ν = 3: 緑、ν = 4: 青である。

`nu` を固定したとき、塗りつぶしの部分の面積が片側 α になるときの `x=t` が上記の`t` 分布の表である。 上記のグラフの例では、`nu=1,x=3.078` であるので、塗りつぶしの部分の面積は表を見て、α=0.10 だとわかる。すなわち、 `t` 分布の表（片側）の数値を `t` とすると、`f_n(x)` と `n, t, alpha` の間には次の関係式がある。

`int_t^oo f_nu(x)dx = alpha`

`t` と `nu` を与えて `alpha` を求める計算を JavaScript で実装し、計算ができるようにした。 `t = 3.078, nu = 1` の場合は、`alpha = 0.0999903 cdots` と計算されるはずだ。
`t` : `nu` : `alpha` :

なお、佐藤郁郎氏による超幾何関数を用いた確率分布の計算(ikuro-kotaro.sakura.ne.jp) も参照のこと。

`F` 分布

統計では、`F` 分布も有用である。`F` 分布を使った検定に `F` 検定がある。`F` 検定は分散の違いの検定や分散分析に使われる。 `F` 分布の表は文献 [2] の巻末などにもあるが、以下の４表は青木先生のページから引用した。

上側確率 `alpha = 0.05` の表

青木先生のＦ分布の上側確率(aoki2.si.gunma-u.ac.jp)から引用した。

`F` 分布の表(`int_alpha^oo G(x; phi, varphi)dx = 0.05`となる `alpha`)

上側確率 `alpha = 0.025` の表

青木先生のＦ分布の上側確率(aoki2.si.gunma-u.ac.jp)から引用した。

`F` 分布の表(`int_alpha^oo G(x; phi, varphi)dx = 0.025`となる `alpha`)

上側確率 `alpha = 0.01` の表

青木先生のＦ分布の上側確率(aoki2.si.gunma-u.ac.jp)から引用した。

`F` 分布の表(`int_alpha^oo G(x; phi, varphi)dx = 0.01`となる `alpha`)

上側確率 `alpha = 0.005` の表

青木先生のＦ分布の上側確率(aoki2.si.gunma-u.ac.jp)から引用した。

`F` 分布の表(`int_alpha^oo G(x; phi, varphi)dx = 0.005`となる `alpha`)

`F` 分布の解説

自由度が `n_1` のカイ２乗分布に従う変数 `x_1` と、自由度が `n_2` のカイ２乗分布に従う変数 `x_2` に関して、次の

`y = (x_1//n_1)/(x_2//n_2)`
が従う確率分布を自由度 `(n_1, n_2)` の `F` 分布という。`F` 分布をとる確率密度関数 `G(x; n_1, n_2)` は、`x le 0` のとき 0 で、`x gt 0` のとき
`G(x; n_1, n_2) = (n_1^(n_1//2) n_2^(n_2//2) x^((n_1-2)//2))/(B(n_1/2, n_2/2)(n_1 x + n_2)^((n_1+n_2)//2))`

の形をとる。ここで `B(*,*)` はベータ関数である。

`alpha = int_t^oo G(x; n_1, n_2)dx`

上記の式で、`t` と `n_1`, `n_2` を与えて `alpha` を求める計算は検定で頻出する。この計算を JavaScript で実装した。 `t = 3.1, n_1 = 3, n_2 = 20` の場合は、`alpha = 0.04992453304567901` と計算されるはずだ。
`t` : `n_1` : `n_2` : `alpha` :

回帰分析

回帰分析は統計の代表的な手法だ。回帰分析のパラメータを得る方法には、 正規方程式と呼ばれる連立一次方程式を計算する方法、QR 分解を使う方法、特異値分解を使う方法がある。 正規方程式の解法について、以下述べる。

修正コレスキー分解

正規方程式を解くための方法として、コレスキー分解を使う方法がある。 ただし、オリジナルのコレスキー分解では対角成分を求めるために平方根を求める計算が必要である。 そこで修正コレスキー分解という、平方根を使わない方法を紹介する。 これは、対称正値行列 `A` を次のような行列の積に分解する方法である。 ここで `L` は下三角行列、`D` は対角行列、そして `{::}^tL` は下三角行列 `L` の転置である。

`A = LD{::}^tL`
既知の `A` と `b` から未知の `x` を求めるには、次の連立一次方程式
`Ax = b`
をとけばよい。それには、
`L(D{::}^tLx) = b`
であることから、`D{::}^tLx = y` とおいて、`Ly = b` となる `y` を求めたのち、`D{::}^tLx = y` から `x` を求める。 具体的には、
`[[l_(11),0,cdots ,0],[l_(21),l_(22),cdots,0],[vdots,vdots,cdots,0], [l_(n1),l_(n2),cdots,l_(n n)]][[y_1],[y_2],[vdots],[y_n]] = [[b_1],[b_2],[vdots],[b_n]]`
から `y` が、
`[[y_1],[y_2],[vdots],[y_n]] = [[b_1 // l_(11)],[(b_2 - l_(21)b_1) // l_(22)],[vdots], [(b_n - l_(n1)b_1 - cdots - l_(n,n-1)b_(n-1)) // l_(n n)]]`
のように、`y_1, y_2, cdots, y_n` の順で求められる（前進消去、elimino antaŭen）。最後に `x` は、
`[[l_(11),cdots,l_(n-1,1) ,l_(n1)], [vdots,ddots,vdots,vdots], [0,cdots,l_(n-1,n-1),l_(n,n-1)], [0,cdots,0,l_(n n)]] [[x_1],[vdots],[x_(n-1)],[x_n]] = [[y_1 // d_1],[vdots],[y_(n-1) // d_(n-1)],[y_n // d_n]]`
であることから、
`[[x_1],[vdots],[x_(n-1)],[x_n]] = [ [(y_1//d_1 - l_(n1)x_n - cdots - l_(21)x_2) // l_(11)], [vdots],[(y_(n-1)//d_(n-1) - l_(n,n-1)x_n) // l_(n-1,n-1)], [(y_n // d_n )// l_(n n)]]`
のように、`x_n, x_(n-1), cdots, x_1` の順で求められる（後退代入、rea anstataŭo）。

`L` の対角成分 `l_(ii)` と `D` の対角成分 `d_(ii)` の決め方は一意ではない。 よく見られるのは `l_(ii) = 1` とする方法だが、 `l_(ii)d_(ii) = 1` とする方法もある。後者は前者よりプログラムの行数が短くなる。 以下の例はすべて後者の分解を作う。 実際の分解のアルゴリズムについては、 連立1次方程式：三角分解 (pbcglba.jp) を参照のこと。

例

問題1. 次の対称行列を修正コレスキー分解せよ。ただし、 分解後の対角行列の対角成分と下三角行列の対角成分の積が 1 に等しくなるようにすること。

(1) `[[1,3],[3,11]]` 。(2) `[[1,-2,3],[-2,5,-4],[3,-4,14]]` 。(3) `[[1,1,0,0],[1,2,1,0],[0,1,2,1],[0,0,1,2]]`

解答1.

(1) `[[1,0],[3,2]] [[1,0],[0,1/2]] [[1,3],[0,2]] ` 。 (2) `[[1,0,0],[-2,1,0],[3,2,1]] [[1,0,0],[0,1,0],[0,0,1]] [[1,-2,3],[0,1,2],[0,0,1]]` 。 (3) `[[1,0,0,0],[1,1,0,0],[0,1,1,0],[0,0,1,1]] [[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]] [[1,1,0,0],[0,1,1,0],[0,0,1,1],[0,0,0,1]] `。

これらの解答で、対角行列の成分がほとんど 1 となっているが、これは計算が楽なように問題を作ったからだろう。 もとは吉田耕作・加藤敏夫の「応用数学」からとった。

問題2. 次の連立一次方程式を解け。

(1) `[[1,3],[3,11]]x = [[26],[7]]`。 (2) `[[1,-2,3],[-2,5,-4],[3,-4,14]]x = [[12],[-19],[48]]` 。 (3) `[[1,1,0,0],[1,2,1,0],[0,1,2,1],[0,0,1,2]] x = [[5],[8],[11],[14]]`

解答2.

(1) `x=[[5],[7]]` 。 (2) `x=[[8],[1],[2]]` 。 (3) `x=[[4],[1],[2],[6]]`。

これらの解は、解答1. の修正コレスキー分解を使って求めている。検算のために、 `Ly = b` となる `y` 、すなわち `D{::}^tLx = y` となる `y` の値も掲げる。

(1) `y=[[26],[7]]` 。 (2) `y=[[12],[5],[2]]` 。 (3) `y=[[5],[3],[8],[6]]`。

問題3. 次の連立一次方程式を解け。

`[[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20]]x = [[34],[68],[96],[125]]`。

解答3.

`x=[[1],[2],[3],[4]]` 。

この例は下記サイトから拝借した。
http://fornext1119.web.fc2.com/9903/0808.html
なお、上記サイトの説明した下三角行列と本ページでは下三角行列と対角行列の定義が異なるが、 同じ解を得る。

問題4. 次の連立一次方程式を解け。

`[[11, 5.5, 3.85],[5.5, 3.85, 3.025],[3.85, 3.025,2.5333]]x = [[1.430],[1.244],[1.11298]]`。

解答4.

`x=[[0.0354895],[-0.492051],[0.9729603]]` 。

この例は「情報量統計学」pp.133-134 の例から拝借した。
同書によれば解が `x = {::}^t[0.03582, -0.49218, 0.97237]` となっているが、 同書の計算結果は正規方程式の前の原データから得ているのかもしれない。

`A` `x` `=` `b`

`x`参考 `y`参考 分解後 `A`

QR 分解

QR 分解は行列の固有値や最小二乗法などで使われる。ここでは最小二乗法のみに焦点を当てる。

例題 5

次の行列 `A` を QR 分解せよ。（Wikipedia より）

`A = [[12,-51,4],[6,167,-68],[-4,24,-41]]`

解答 5

下記は Scilab によって求めた値で一例である。

`Q = [[-0.8571429, 0.3942857, 0.3314286],[-0.4285714, -0.9028571, -0.0342857],[0.2857143, -0.1714286, 0.9428571]]`

`R = [[-14, -21, 14],[0, -175, 70],[0, 0, -35]]`

Wikipedia と比べると、`Q` は第１列と第２列の要素すべてで符号が反対である。`R` は第３行第３列の要素のみ符号が反対である。

下記の[問題5-1]ボタンをクリックすると、上記の `A` がテキストエリアにコピーされる。[QR 分解計算]ボタンをクリックすると R のみ求める（Qは求めない）。 本ページの計算で求めた `Q` の値は、第３行第３列のみ Scilab と逆である。

`R`

数式とグラフ

数式の表現には ASCIIMath を用いている。

グラフは 関数グラフ描画スクリプト・SVGGraph.js (defghi1977.html.xdomain.jp) を使っている。

参考文献

  1. 久米均:統計解析への出発
  2. 小針晛宏:確率・統計入門

まりんきょ学問所 ＞ JavaScript 手習い ＞ 統計計算

MARUYAMA Satosi