統計計算

作成日 : 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` 分布の表
α
ν 0.20 0.10 0.05 0.02 0.01 0.002 0.001 両側
0.10 0.05 0.025 0.01 0.005 0.001 0.0005 片側
1 3.078 6.314 12.706 31.821 63.657 318.309 636.619
2 1.886 2.920 4.303 6.965 9.925 22.327 31.599
3 1.638 2.353 3.182 4.541 5.841 10.215 12.924
4 1.533 2.132 2.776 3.747 4.604 7.173 8.610

出典は青木先生の下記ページ。
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) も参照のこと。

回帰分析

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

修正コレスキー分解

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

`A = LD{::}^tL`
既知の `A` と `b` から未知の `x` を求めるには、次の連立一次方程式
`Ax = b`
をとけばよい。それには、
`L(D{::}^tLx) = b`
であることから、`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_nl_(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` とする方法もある。後者は前者よりプログラムが短くなる。 以下の例はすべて後者の分解を作う。 実際の分解のアルゴリズムについては下記を参照のこと。
http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?Numerical%20Calculation
ここの「連立1次方程式：三角分解」にある。

例

問題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]` となっているが、 同書の計算結果は正規方程式の前の原データから得ているのかもしれない。

数式とグラフ

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

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

参考文献

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

MARUYAMA Satosi