統計処理と言えるほどではないが、データ処理して中央値や最頻値を眺めていて思った。
ある程度本格的に統計処理する人はR言語とかOctaveとかPython*1とかを使うと思うし、そうでない人もPerlやRubyなどのスクリプト言語を使うと思う。
シェルスクリプト使いもawk(1)を併用するはずだ。awk(1)は汎用の言語ではないものの、シェルスクリプトで使用する他のコマンドと比較すればスクリプト言語に近いツールだ。
では、awk(1)抜きで頑張るとどうなるだろうか？ seq(1)が吐き出すような改行区切りの数列のテキストを扱うという前提で考えてみた。
ちなみに以下のスクリプトは入力値の検証を全く行っていない。
最小値・最大値
最小値を求めるのは簡単だ。昇順にソートして先頭の値を取り出せばよい。
#!/bin/sh exec sort -n ${@+"$@"} | head -1
降順にソートすれば最大値を得られる。
#!/bin/sh exec sort -n -r ${@+"$@"} | head -1
算術平均
算術平均は、awk(1)を使うならこんな感じ。
#!/usr/bin/awk -f { total += $1 } END { if (NR >= 1) { print total / NR } }
シンプルだ。awk(1)のサンプルに出てきそうな感じである。
これをawk(1)なしで実現できるのか？ ヒント：配列は不要。
#!/bin/sh exec sed ' 1i\ scale = 16\ t = 0\ n = 0 s/^.*$/t += &; n++/ $a\ t / n' ${@+"$@"} | bc -q -s | tail -1
sed(1)でbc(1)用の計算式を組み立てて流し込めばOK。bc(1)には配列的な機能はないが変数はあるので、力づくな式で算術平均を計算することができる。
ちなみに、while readとexpr(1)の組み合わせでも実現できそうだが、データ数が大きくなるほど激遅だろうし、最後の除算で小数部を表示したいならbc(1)を使わなくてはならない。
標準偏差
標準偏差の求め方は、これまた力技になるが算術平均のそれの応用でOK。
#!/bin/sh readonly progname=`basename "$0"` trap 'rm -f "$tmpfile"; exit 1' 1 2 15 tmpfile=`mktemp "$progname.$$.XXXXXXXXXX"` if [ $? -ne 0 ]; then echo "$progname: cannot create temporary file" 1>&2 exit 1 fi cat ${@+"$@"} >"$tmpfile" { sed ' 1i\ scale = 16\ t = 0\ n = 0 s/^.*$/t += &; n++/ $a\ m = t / n' "$tmpfile" sed ' 1i\ t = 0 s/^.*$/t += (& - m) ^ 2/ $a\ sqrt(t / n)' "$tmpfile" } | bc -q -s | tail -1 rm -f "$tmpfile"
算術平均を使う関係でデータを2回舐めることになるので、一時ファイルを使用するのがポイント。算術平均の計算式と標準偏差の計算式を分けて生成し、まとめてbc(1)に流し込んでいる。
ちなみに、式の最後で平方根を求めるのを止めれば分散を得られるはず。
中央値
中央値をawk(1)で求めるなら、こんな感じだろうか？
#!/usr/bin/awk -f { samples[NR] = $1 } END { if (NR >= 1) { mid = int(NR / 2) print (NR%2 == 1) ? samples[mid+1] : (samples[mid] + samples[mid+1]) / 2 } }
これまた、そこそこシンプルである。
これをawk(1)なしで実現するにはどうすればよいのか？
#!/bin/sh readonly progname=`basename "$0"` trap 'rm -f "$tmpfile"; exit 1' 1 2 15 tmpfile=`mktemp "$progname.$$.XXXXXXXXXX"` if [ $? -ne 0 ]; then echo "$progname: cannot create temporary file" 1>&2 exit 1 fi nsamples=`sort -n ${@+"$@"} | tee "$tmpfile" | wc -l` [ "$nsamples" -gt 0 ] || do_exit 0 mid=`expr $nsamples / 2` if [ `expr $nsamples % 2` = 1 ]; then sed -n "`expr $mid + 1`p" "$tmpfile" else sed -n "$mid,`expr $mid + 1`p" "$tmpfile" | paste -d '+' -s | sed 's%^.*$%scale = 16; (&) / 2%' | bc -q -s fi rm -f "$tmpfile"
中央値を求めるにはソートされたデータ一式が必要だが、配列の代わりに一時ファイルを使用する。データ数は一時ファイルの行数と同じになる。
データ数が奇数の場合は、一時ファイルの中央の行を出力すればよい。偶数の場合は中央の値2つの算術平均を求める必要があるので、計算式を組み立ててbc(1)に流し込んでいる。
ちなみにこのシェルスクリプトで初めてpaste(1)を使った。
頻度の一覧
値の頻度の一覧を出力するのは簡単だった。
#!/bin/sh exec sort -n ${@+"$@"} | uniq -c
そう、uniq(1)があったじゃないか。
最頻値
頻度を求めるのは簡単だが、最頻値の場合は一工夫必要だ。
シェルスクリプトで実現するにしても、この例のように部分的にawk(1)を使えると非常に楽だ。
#!/bin/sh exec sort -n ${@+"$@"} | uniq -c | sort -k 1 -n -r | awk 'NR == 1 { f = $1 } $1 != f { exit } { print }' | sort -k 2 -n
uniq(1)で頻度を求めたら、頻度の値の降順でソートする。そうすると最も頻度の大きい値の順になる。ここで、先頭の値の頻度の値（つまり最も大きな頻度の値）と同じ頻度のデータのみを取り出すのにawk(1)を使用すれば、このように1行で解決できる。
このawk(1)の部分を他のコマンドの組み合わせで実現すればよいのだが……。
#!/bin/sh readonly progname=`basename "$0"` trap 'rm -f "$tmpfile"; exit 1' 1 2 15 tmpfile=`mktemp "$progname.$$.XXXXXXXXXX"` if [ $? -ne 0 ]; then echo "$progname: cannot create temporary file" 1>&2 exit 1 fi sort -n ${@+"$@"} | uniq -c | sort -k 1 -n -r | sed 's/^[ \t]*//' >"$tmpfile" freq=`head -1 "$tmpfile" | tr ' ' '\t' | cut -f 1` sed -n " /^$freq/p /^$freq/!q" "$tmpfile" | sort -k 2 -n rm -f "$tmpfile"
こうなった。
おわりに
パーセンタイル順位の計算も実現できそうな気がするが、飽きたのでここで中断する。
*1:『Think Stats』からの連想。
