awk で最小二乗法

私が入社した十年以上前は Linux はまだまだマイナーで企業の中には SystemV 系の OS が多く、しかも開発環境は gcc も含めて全く用意されていませんでした。 今ではそんな状況を考えることも難しいのですが、そんな中でも「UNIX の三種の神器」である grep, sed, awk は使うことができ、いろいろなことに役立ってくれたのを覚えています。 多分、入社して最初に作った awk スクリプトは最小二乗法で 1 次関数近似を行うものだったと記憶しています。

今は Excel などで簡単に近似関数を求めることができ、しかもグラフまで描かせることができますが、当時は awk + Ngraph というのが極めて普通のグラフ作成スタイルでした。 Excel は一度に読み込める行数が少ないので、数 MB のデータを扱う際には、awk + gnuplot で描かせることもあります。 そんなことはさておき、1 次関数近似の最小二乗法を解いてみます。

一次方程式への近似に式が書かれていますが、このとおりに記述します。 また、データはタブまたは複数のスペースで区切られているものとします。

#! /usr/bin/gawk -f
# least_square.awk

{
    x[NR] = $1;
    y[NR] = $2;
}

END {

    # データの総和
    n = NR;

    # xy の和を求める
    for (i in x) {
        sum_xy += x[i] * y[i];
    }

    # x の和を求める
    for (i in x) {
        sum_x += x[i];
    }

    # y の和を求める
    for (i in y) {
        sum_y += y[i];
    }

    # x ^ 2 の和を求める
    for (i in x) {
        sum_x2 += x[i] ^ 2;
    }

    # y = ax + b に近似する場合の a, b を求める
    a = (n * sum_xy - sum_x * sum_y) / (n * sum_x2 - sum_x ^ 2);
    b = (sum_x2 * sum_y - sum_xy * sum_x) / (n * sum_x2 - sum_x ^ 2);

    print "y = " a " * x + " b;
}

ここで用いるデータは以下のようなものです。

1   1
2   2
3   3

あくまでサンプルですので、少ないデータですが、実際には巨大なデータを扱うこともできます。

$ gawk -f least_square.awk sample.txt
y = 1 * x + 0

ちゃんと解くことができました。

ここで意地悪をしてみますが、x, y ともに awk の rand() 関数を用いるとどうなるでしょうか? 完全に乱数であれば、傾向を持たないため、係数 a は 0 になり、0 から 1 までの乱数ですから切片 b は 0.5 付近に収束すると思われます。 これを試してみましょう。

$ gawk 'BEGIN{for(i=1;i<=100;i++){print rand(), rand()}}' | gawk -f least_square.awk
y = 0.105563 * x + 0.414944
$ gawk 'BEGIN{for(i=1;i<=1000;i++){print rand(), rand()}}' | gawk -f least_square.awk
y = -0.0195874 * x + 0.496436
$ gawk 'BEGIN{for(i=1;i<=10000;i++){print rand(), rand()}}' | gawk -f least_square.awk
y = -0.00758212 * x + 0.501733
$ gawk 'BEGIN{for(i=1;i<=100000;i++){print rand(), rand()}}' | gawk -f least_square.awk
y = 0.00241209 * x + 0.498356

確かに予想は当たっているようです。

awk の得意とするところのひとつにデータ処理があります。 小回りのきく awk ならではの簡単な処理を行うことができます。

tag_nawk.pngtag_nawk.pngtag_nawk.pngtag_nawk.png