適当な数の一様分布の平均を取ることで、その分布がガウス分布に従うことを調べました。
ターミナル上で数字遊びをする練習と思ってみていただけたら幸いです。
ガウス分布に従う数を生成
0-1のランダムな数を生成
これは次の３ステップで行なった。
- ランダムな（無限につづく）数字の並びを生成
- 適当な数（欲しい桁数）で折り返す
- 各行の先頭に 0 をつける
以下順番に確認する。
ランダムな（無限につづく）数字の並びを生成
ランダムな数の生成には、/dev/urandomを使う。これ は乱数を返すプログラム的なものであが、ランダムなバイト列であり、文字コードにないものは文字化けする。
そこでtrコマンドを用いて必要なものだけ抽出している。
必要なオプションは次の2つ
- -d: (delete) 指定したものを削除
- -c: (complement) 引数以外を
trの対象にする
Complement は難しい英単語ですが、数学でいう補集合のことです。
# head をつけないと永遠に謎の文字が流れ続けるので注意！ $ cat /dev/random| head -c 50 ???-???3????9??Y^]<?&????35????/?S#@???????a??G? # tr を用いて数字だけを取り出す ## MACの場合は、以下のように LC_ALL=C が必要 $ cat /dev/urandom | LC_ALL=C tr -dc 0-9 | head -c 50 27690019162627559341647523465681647020178540856142
適当な数（欲しい桁数）で折り返す
foldコマンドを使う。
$ cat /dev/urandom | LC_ALL=C tr -dc 0-9 | fold -w 5 73632 89958 55285 ^C
ここまでで、0 - 99999 の乱数が生成できたこととなる
0-1 のランダム変数に変える
各行の先頭に 0 をつけ、変数の範囲を 0~1 に変更する。
ここではsedコマンドを用いている。
cat /dev/urandom | LC_ALL=C tr -dc 0-9 | fold -w 5 | sed 's@^@0.@' > 0.73632 > 0.89958 > 0.55285 > ^C
これで 0-1 の一様分布が完成した。
ガウス分布に持っていく
今回は試しにN＝12個の平均とったもの（ガウス分布に従うと信じているもの）を１万個生成する。これは次の２ステップで行った。
- 平均取りたい数だけ、一様分布を横に並べる
- 各行で平均を取る
12個まず、横に並べる
xargs を使う。これ色々と便利。
$ cat /dev/urandom | LC_ALL=C tr -dc 0-9 | fold -w 5 | sed 's@^@0.@' | xargs -n 12 | head -n 10000
横に並べた12個の平均をとる
awk の中が長くて混乱するかもしれないが、シンプルに考えていくだけなので省略。
$ cat /dev/urandom | LC_ALL=C tr -dc 0-9 | fold -w 5 | sed 's@^@0.@' | xargs -n 12 | head -n 10000 | \
awk '{for(i=0; i<=int(NF) ;i++){{if(i==0){a = 0}else{a += $i}}{if(i == NF){print 2*a/NF}}}}' > gauss
ガウス分布に従う数を１万個生成できた！
これをヒストグラムにしてプロットしてみる
初めは、ソートした後にawkなどを用いてもしこの数が0.1を超えてたら〜〜とかでやろうとしたけど、while文しんどいなーと思い方針を変更
そもそも 0.1 刻みなら小数点第１桁までわかればいいということに気づく。
任意の刻みで丸め込む
まず簡単な0.1刻み
文字列の1-3つ目だけをとってくればいい
cat gauss | awk '{print substr($0,1,3)}'
> 1.2
> 1.1
> 1.0
> 0.9
ここから更に4倍細かくみたければ、次のようにする
cat gauss | awk '{print substr(4*$0,1,3)/4}'
> 1.15
> 1.10
> 0.825
ソートして、同じ数値のものをカウントする
cat gauss | awk '{print substr(6*$0,1,3)/6}' | sort | uniq -c
> 5 0.5
> 3 0.516667
> 5 0.533333
> 8 0.55
> 8 0.56666
数値とカウントを入れ替える
これは完全に、のちのプロットのしやすさのため
一列目と二列目を入れ替えるには次のように awk を使う
cat gauss | awk '{print substr(6*$0,1,3)/6}' | sort | uniq -c | awk '{print $2,$1}'
> 0.5 5
> 0.516667 3
> 0.533333 5
> 0.55 8
> 0.566667 8
最後にファイルに出力しておく
cat gauss | awk '{print substr(6*$0,1,3)/6}' | sort | uniq -c | awk '{print $2,$1}' > gauss_hist
平均をとる数（N）を変えて、色々と実験
下の図において、gauss は N=12 であり、gauss_n は 「N=12*n」と N を変えて実験したもの。
Nを増やすほど綺麗な正規分布に近づいていることがわかる。
また、標準偏差も小さくなっている。
これは自然科学において「実験を繰り返し統計をためることが重要である」ことと対応してる。
おまけ
0-1 の一様平均が正しいのかの確認
前半で、
これで 0-1 の一様分布が完成した。
と自信満々に言ったが、果たしてこれは正しいのか確認しておく。
以下は、0.05 刻みで 0-1 を分割し、それぞれの分布を調べるためのものである。
$ cat /dev/urandom | LC_ALL=C tr -dc 0-9 | fold -w 5 | sed 's@^@0.@' | head -n 100000 | \ awk '{a[int($1/0.05)]++}END{for(i=0; i<20; i++){print 0.05*i"~"0.05*(i+1),a[i]}}' | tee tmp 0~0.05 4934 0.05~0.1 4958 0.1~0.15 5104 0.15~0.2 5007 0.2~0.25 5012 0.25~0.3 5090 0.3~0.35 4954 0.35~0.4 4784 0.4~0.45 5117 0.45~0.5 5050 0.5~0.55 4929 0.55~0.6 4994 0.6~0.65 4976 0.65~0.7 5034 0.7~0.75 4994 0.75~0.8 5083 0.8~0.85 4989 0.85~0.9 4957 0.9~0.95 5104 0.95~1 4930
ぱっと見正しいように見えるが、これは本当に一様分布と言っていいのか？
それを調べるための統計的手法が、カイ二乗検定である。
カイ二乗検定
カイ二乗検定とは、「観測された事象の相対頻度が、ある頻度分布に従う」という帰無仮説を検定するためのものであり、カイ二乗分布が使われる。カイ二乗分布に必要なパラメータとしては、カテゴリ数 k と有意水準 α である。
具体的な計算
観測された値から計算されるカイ二乗値
を計算し、自由度 k-1 のカイ二乗分布 χ2(k-1, α) と比較しその大小で求める。
今回は、上で 20 個にカテゴライズした数値を用い（k=20）、有意水準は 0.01 として計算してみる。つまり
# カイ二乗の値を求める。期待度数はどのブロックでも 5000。 $ cat tmp | awk '{sum += ($2-5000)**2/5000}END{print "chi_sq_0",sum}' chi_sq_0 24.322
この値は、χ2(19, 0.01) ~ 36.19 よりも小さく、有意水準 1% で棄却されない。
よって、偏りがあるとは言えない。
