ある平均値の正規分布のサンプル

ときどきの雑記帖 濫觴編 にインスパイヤされて、ある平均値を持つ正規分布のサンプルを作成してみます。 正規分布をするサンプルが欲しい時は多いのですが、なかなかサンプルを作ろうとすると大変だったりします。

正規乱数 というのを用いると、うまく作れそうです。 ここでは同じく、ボックス=ミューラー法 (Box-Muller transform) で求めてみます。

#! /usr/local/bin/nawk -f
# normal_distributuon.awk
# 正規分布をする 0 から 100 までの整数の乱数を N 個発生させる
# usage: nawk -f normal_distribution.awk

BEGIN {
    srand();

    sample_num = sample_num ? sample_num : 100000;
    average    = average    ? average    : 50;
    stdev      = stdev      ? stdev      : 10;
    max        = max        ? max        : 100;

    for (i = 1; i <= sample_num; i++) {
        score = (i % 2) ? \
                int(stdev * sqrt(-2 * log(rand())) * cos(2 * pi() * rand()) + average) : \
                int(stdev * sqrt(-2 * log(rand())) * sin(2 * pi() * rand()) + average) ;

        if (score >= 0 && score <= max) {
            print score;
        }
    }
}

# pi():     PI を計算する
#   in:     なし
#   out:    PI
function pi() {
    return 4 * atan2(1, 1);
}

さて、実行してみましょう。

$ gawk -f normal_distribution.awk
40
50
58
37
45
42
69
69
38
42
<snip>

と表示されますが、本当に正規分布しているかどうかこれだけでは分かりません。 そこで、以下のような集計プログラムを用意しました。

#! /usr/local/bin/nawk -f
# histgram.awk
# ヒストグラムを作成する
# usage: nawk -f histgram.awk foo.txt

BEGIN {
    max = -1 * 2 ^ 31;
    min =      2 ^ 31;

    hist_div  = hist_div  ? hist_div  : 10;
    bar_width = bar_width ? bar_width : 50;
}

{
    line[NR] = $0;

    # 合計の計算
    sum += line[NR];

    # 平均の計算
    ave = sum / NR;

    # 最大値の計算
    if (max < line[NR]) {
        max = line[NR];
    }

    # 最小値の計算
    if (min > line[NR]) {
        min = line[NR];
    }
}

END {

    for (i = 1; i <= NR; i++) {
        sum2 += (line[i] - ave) ^ 2;
    }
    stdev  = sqrt(sum2 / (NR - 1));

    print "-------------------------------------------------";
    print "Number of Data: " NR;
    print "Summation     : " sum;
    print "Average       : " ave;
    print "Maximum       : " max;
    print "Minimum       : " min;
    print "Std. Deviation: " stdev;
    print "-------------------------------------------------";

    for (i = 1; i <= NR; i++) {
        for (j = 1; j <= hist_div; j++) {
            if (line[i] >= min + (max - min) / hist_div * (j - 1) &&
                line[i] <  min + (max - min) / hist_div * j &&
                j != hist_div) {
                num[j]++;
            }
            if (line[i] >= min + (max - min) / hist_div * (j - 1) &&
                line[i] <  min + (max - min) / hist_div * j &&
                j == hist_div) {
                num[j]++;
            }
        }
    }

    for (i in num) {
        if (max_num < num[i]) {
            max_num = num[i];
        }
    }

    for (i = 1; i <= hist_div; i++) {
        printf("%8s-%8s %s (%s)\n",
                min + (max - min) / hist_div * (i - 1),
                min + (max - min) / hist_div * i,
                repeat_str("#", num[i] / max_num * bar_width),
                num[i]);
    }
}

# repeat_str():     文字列を num 回連接する関数
#   in:     文字列 str
#           回数 num
#   out:    文字列 str を num 回連接させた文字列
function repeat_str(str, num,    i, result) {
    for (i = 1; i <= num; i++) {
        result = sprintf("%s%s", result, str);
    }

    return result
}

パイプで接続してみましょう。

$ gawk -f normal_distribution.awk | gawk -f histgram.awk
-------------------------------------------------
Number of Data: 100000
Summation     : 4950854
Average       : 49.5085
Maximum       : 95
Minimum       : 6
Std. Deviation: 9.98371
-------------------------------------------------
       6-    14.9  (32)
    14.9-    23.8  (444)
    23.8-    32.7 ##### (3895)
    32.7-    41.6 ######################### (16714)
    41.6-    50.5 ################################################## (32967)
    50.5-    59.4 ############################################# (30092)
    59.4-    68.3 ################### (13011)
    68.3-    77.2 ### (2595)
    77.2-    86.1  (237)
    86.1-      95  (12)

上から、個数、合計、平均、最大、最小、標準偏差を示しています。 確かにヒストグラムからも、正規分布のようですが、分割数を 30 にしてみます。

$ gawk -f normal_distribution.awk | gawk -v hist_div=30 -f histgram.awk
-------------------------------------------------
Number of Data: 99999
Summation     : 4951475
Average       : 49.5152
Maximum       : 97
Minimum       : 7
Std. Deviation: 10.0141
-------------------------------------------------
       7-      10  (7)
      10-      13  (9)
      13-      16  (18)
      16-      19  (50)
      19-      22  (183)
      22-      25 # (379)
      25-      28 ### (771)
      28-      31 ###### (1454)
      31-      34 ########## (2562)
      34-      37 ################# (4222)
      37-      40 ########################## (6188)
      40-      43 ################################## (8296)
      43-      46 ########################################### (10256)
      46-      49 ################################################ (11606)
      49-      52 ################################################## (11871)
      52-      55 ############################################### (11174)
      55-      58 ######################################### (9801)
      58-      61 ################################ (7620)
      61-      64 ###################### (5422)
      64-      67 ############### (3597)
      67-      70 ######### (2232)
      70-      73 #### (1161)
      73-      76 ## (640)
      76-      79 # (265)
      79-      82  (139)
      82-      85  (49)
      85-      88  (16)
      88-      91  (6)
      91-      94  (3)
      94-      97  (1)

乱数で発生させているので、これだけきれいになることも少ないのですが、きれいな正規分布になっているようです。

tag_nawk.png tag_nawk.png tag_nawk.png tag_nawk.png