2011年 8月 1日 掲載
「AWK で十分ですよ，わかってくださいよ」
の，第2回である． 今回も基本的な統計量をあつかう．
まずは，サンプルデータ作りである． 今回は 正規乱数列を作る． Box-Muller法で n個の正規乱数列を作るスクリプト normal_random_numbers.awk は，以下のように書ける．
BEGIN{
pi = atan2(0, -1)
while(i < n){
x1 = rand()
x2 = rand()
if(x1 != 0){
print sqrt(-2 * log(x1)) * cos(2 * pi * x2)
i++
if(i >= n) break
print sqrt(-2 * log(x1)) * sin(2 * pi * x2)
i++
}
}
}
このスクリプトを実行するときは， -vオプションで，乱数の個数 n を指定する必要がある． 今回は，100万個の乱数列 nrand.txt を作る．
$ awk -f normal_random_numbers.awk -v n=1000000 > nrand.txt
前回作ったスクリプトで，平均値，最大値，最小値を出してみよう．
$ awk -f average.awk nrand.txt -0.00055316 $ awk -f max.awk nrand.txt 4.91464 $ awk -f min.awk nrand.txt -4.56578
まずは，nrand.txt の分散を計算してみよう． 分散の公式は，
だから，
{
sum_x += $1
sum_x2 += $1 ^ 2
}
END{
if(NR == 0) exit
print sum_x2 / NR - (sum_x / NR) ^ 2
}
と書けば，だいたいの場合は OK だが， このスクリプトでは，計算誤差が蓄積して負の値が出てきたりすることがある． めんどうでも，一度平均値を計算してから，差の 2乗の平均を求めた方が安全である．
分散を計算するスクリプト variance.awk は，以下のように書ける．
{
x[NR] = $1
}
END{
if(NR == 0) exit
for(i in x){
sum_x += x[i]
}
m_x = sum_x / NR
for(i in x){
sum_dx2 += (x[i] - m_x) ^ 2
}
print sum_dx2 / NR
}
nrand.txt の分散を計算してみよう．
$ awk -f variance.awk nrand.txt 1.00036
分散の平方根を計算すれば，標準偏差である．
標準偏差を計算するスクリプト standard_deviation.awk は，
{
x[NR] = $1
}
END{
if(NR == 0) exit
for(i in x){
sum_x += x[i]
}
m_x = sum_x / NR
for(i in x){
sum_dx2 += (x[i] - m_x) ^ 2
}
print sqrt(sum_dx2 / NR)
}
となる．
nrand.txt の標準偏差は，
$ awk -f standard_deviation.awk nrand.txt 1.00018
となる．
nrand.txt が，平均 0, 標準偏差 1 の分布に，ほぼ従っていることが確認できた．
nrand.txt の度数分布を確認しておきたい． 度数分布を計算するスクリプト frequency_distribution.awk は，以下のように書ける．
{
if(binsize <= 0) exit
if($1 < 0){
frequency[int($1 / binsize) - 1] ++
}else{
frequency[int($1 / binsize)] ++
}
}
END{
for(i in frequency){
print i * binsize, frequency[i] | "sort -n"
}
}
実行時に，-vオプションで，階級の幅 binsize を指定する． binsize=0.2 の場合，以下のようになった．
$ awk -f frequency_distribution.awk -v binsize=0.2 nrand.txt -4.6 5 -4.4 7 -4.2 30 -4 42 -3.8 103 -3.6 165 -3.4 387 -3.2 685 -3 1253 -2.8 2103 -2.6 3455 -2.4 5651 -2.2 8977 -2 12990 -1.8 18893 -1.6 26096 -1.4 34499 -1.2 43422 -1 53128 -0.8 62540 -0.6 70494 -0.4 75926 -0.2 79247 0 79214 0.2 76073 0.4 70484 0.6 62322 0.8 53359 1 43806 1.2 33966 1.4 25834 1.6 18859 1.8 13196 2 8900 2.2 5767 2.4 3538 2.6 2126 2.8 1155 3 652 3.2 324 3.4 161 3.6 87 3.8 42 4 21 4.2 9 4.4 4 4.6 2 4.8 1
これだけでも，分布の雰囲気はだいたいわかるが， できればヒストグラムにして見てみたい． そこで，1次元のデータを簡単にプロットするスクリプト plot1.awk を作ってみた．
BEGIN{
width = 50
}
{
x[NR] = $1
}
END{
xmin = xmax = x[1]
for(i in x){
if(x[i] < xmin) xmin = x[i]
if(x[i] > xmax) xmax = x[i]
}
if(xmin == xmax) xmax = xmin + 1
for(i = 1; i <= NR; i++){
ix = int(width * (x[i] - xmin) / (xmax - xmin))
for(j = 0; j < ix; j++) printf " "
printf "*\n"
}
}
width の値は、環境にあわせて適当に変更してほしい。
nrand.txt のヒストグラムは，以下のようにして描ける．
$ awk -f frequency_distribution.awk -v binsize=0.2 nrand.txt | awk '{print $2}' | awk -f plot1.awk
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
きれいな釣鐘カーブになった．
今回の AWKスクリプト：
[ KDO ホームにもどる ]
(C) 2011 KDO