gawkで基本統計量を計算する

はじめに

基本統計量を計算するshスクリプト(Bash)です。 計算にはgawkを使いました。

基本統計量

統計量	記号	変数	説明	表計算ソフト(Excel)の関数
データ	-	ds[]	データ(標本)	-
データ数	n	n	データ数(標本数)	-
最小値	-	min	最小値	MIN()
最大値	-	max	最大値	MAX()
中央値	-	med	ソートされたデータの中央値 データ数が偶数の場合は中央の２値の平均	MEDIAN
平均値	$ \bar{x} $	ave	平均値(標本平均)	AVERAGE()
データの総和	-	sum	$ \sum_{i=1}^{n} {x_i} $	SUM()
標本分散	$ s^2 $	svar	$ \frac{1}{n}\sum_{i=1}^{n} (x_i-\bar{x})^2 $	VAR.P()
不偏分散	$ u^2 $	uvar	$ \frac{1}{n-1}\sum_{i=1}^{n} (x_i-\bar{x})^2 $	VAR.S()
標本標準偏差	$ s $	ssd	$ \sqrt{\frac{1}{n}\sum_{i=1}^{n} (x_i-\bar{x})^2} $	STDEV.P()
不偏標準偏差	$ u $	usd	$ \sqrt{\frac{1}{n-1}\sum_{i=1}^{n} (x_i-\bar{x})^2} $	STDEV.S()
歪度	-	skew	$ \frac{n}{(n-1)(n-2)}\sum_{i=1}^{n} \frac{(x_i-\bar{x})^3}{u^3} $	SKEW()
尖度	-	kurt	$ \frac{n(n+1)}{(n-1)(n-2)(n-3)}\sum_{i=1}^{n} \frac{(x_i-\bar{x})^4}{u^4}-\frac{3(n-1)^2}{(n-2)(n-3)} $	KURT()

shスクリプトと出力形式

【shスクリプト】

basic_statistics.sh

#!/bin/bash
calc_basic_statistics(){
    gawk '
    function sumf(a, n                 ,sum_p, sum_n, b, d, i){
        # 絶対値が小さいものから加算する。
        sum_p = 0.0;
        sum_n = 0.0;
        asort(a, b);
        for( i = 0; i < n; i++){
            d = b[i+1];
            if( d < 0.0){
                continue;
            }
            else{
                sum_p += d;
            }
        }
        for( i = n-1; i >= 0; i--){
            d = b[i+1];
            if( d > 0.0){
                continue;
            }
            else{
                sum_n += d;
            }
        }
        return  sum_p + sum_n;
    }
    function median(a, n                 ,b, m1, m2){
        if( n == 1){
            return  a[0];
        }
        asort(a, b);
        m1 = int(n/2);
        m2 = int((n-1)/2);
        if( m1 == m2){
            return  b[m1+1];
        }
        else{
            return  (b[m1+1] + b[m2+1])/2;
        }
    }
    BEGIN{
        n = 0;
    }
    {
        # データ
        d = $1;
        ds[n] = d;
        # 最大値、最小値
        if(n == 0){ # 初期値
            max = d;
            min = d;
        }
        else{
            if( d > max){
                max = d;
            }
            if( d < min){
                min = d;
            }
        }
        # データ数
        n += 1;
    }
    END{
        if( n == 0){
            exit 1;
        }
        # データの総和を算出
        sum = sumf(ds, n);

        # 平均値を算出
        ave = sum / n;

        # 偏差
        # 偏差の2乗
        for(i = 0;i < n;i++){
            dev = ds[i] - ave;
            devs[i] = dev;
            dev2s[i] = dev * dev;
        }

        # 偏差の2乗の総和を算出
        sum_dev2 = sumf(dev2s, n);

        # 標本分散と標本標準偏差を算出
        svar = sum_dev2 / n;    # 標本分散
        ssd = sqrt(svar);       # 標本標準偏差

        # 不偏分散と不偏標準偏差を算出
        if( n > 1){
            uvar = sum_dev2 / (n - 1);# 不偏分散
            usd = sqrt(uvar);       # 不偏標準偏差
        }
        else{
            uvar = "@@@";
            usd = "@@@";
        }

        # (偏差/不偏標準偏差)の3乗の総和を算出
        # (偏差/不偏標準偏差)の4乗の総和を算出
        if( usd != "@@@"){
            for(i = 0;i < n;i++){
                dev = devs[i];
                devz3s[i] = (dev / usd)^3;
                devz4s[i] = (dev / usd)^4;
            }
            sumz3 = sumf(devz3s, n);
            sumz4 = sumf(devz4s, n);
        }

        # 歪度
        if( usd != "@@@" && n > 2){
            skew = (n /((n -1 ) * (n - 2))) * sumz3;
        }
        else{
            skew = "@@@";
        }

        # 尖度
        if( usd != "@@@" && n > 3){
            kurt = ((n * (n + 1)) /((n -1 ) * (n - 2) * (n - 3))) * sumz4 - ((3 * (n - 1)^2) / ((n - 2) * (n - 3)));
        }
        else{
            kurt = "@@@";
        }

        # 中央値
        med = median(ds, n);

        printf("%d,", n);
        printf("%G,", min);
        printf("%G,", max);
        printf("%G,", med);
        printf("%G,", ave);
        printf("%G,", sum);
        if( svar != "@@@") {printf("%G,", svar);} else {printf(",");}
        if( uvar != "@@@") {printf("%G,", uvar);} else {printf(",");}
        if( ssd != "@@@") {printf("%G,", ssd);} else {printf(",");}
        if( usd != "@@@") {printf("%G,", usd);} else {printf(",");}
        if( skew != "@@@") {printf("%G,", skew);} else {printf(",");}
        if( kurt != "@@@") {printf("%G,", kurt);} else {printf(",");}

        printf("n=%d\\\\n", n);
        printf("min=%G\\\\n", min);
        printf("max=%G\\\\n", max);
        printf("med=%G\\\\n", med);
        printf("ave=%G\\\\n", ave);
        printf("sum=%G\\\\n", sum);
        if( svar != "@@@") {printf("svar=%G\\\\n", svar);} else {printf("svar=N/A\\\\n");}
        if( uvar != "@@@") {printf("uvar=%G\\\\n", uvar);} else {printf("uvar=N/A\\\\n");}
        if( ssd != "@@@") {printf("ssd=%G\\\\n", ssd);} else {printf("ssd=N/A\\\\n");}
        if( usd != "@@@") {printf("usd=%G\\\\n", usd);} else {printf("usd=N/A\\\\n");}
        if( skew != "@@@") {printf("skew=%G\\\\n", skew);} else {printf("skew=N/A\\\\n");}
        if( kurt != "@@@") {printf("kurt=%G\\\\n", kurt);} else {printf("kurt=N/A\\\\n");}
        printf("\n");

        exit 0;
    }
    ' $@
}

【出力形式】

データ数,最小値,最大値,中央値,平均値,データの総和,標本分散,不偏分散,標本標準偏差,不偏標準偏差,歪度,尖度,テキスト

最後の項目(テキスト)の形式です。
そのままecho -n -eやprintfに渡せば、項目毎に改行されるように、\で\nをエスケープしています。

n=値\\nmin=値\\nmax=値\\nmed=値\\nave=値\\nsum=値\\nsvar=値\\nuvar=値\\nssd=値\\nusd=値\\nskew=値\\nkurt=値\\n

実行例

【データファイル】

q1.txt

0.7
-1.6
-0.2
-1.2
-0.1
3.4
3.7
0.8
0.0
2.0

【その１：データファイルの直接指定】

実行

(. ./basic_statistics.sh; calc_basic_statistics q1.txt)

実行結果

10,-1.6,3.7,0.35,0.75,7.5,2.8805,3.20056,1.6972,1.78901,0.580921,-0.629822,n=10\\nmin=-1.6\\nmax=3.7\\nmed=0.35\\nave=0.75\\nsum=7.5\\nsvar=2.8805\\nuvar=3.20056\\nssd=1.6972\\nusd=1.78901\\nskew=0.580921\\nkurt=-0.629822\\n

出力の整形

実行

(. ./basic_statistics.sh; calc_basic_statistics q1.txt) | (IFS=',' read x x x x x x x x x x x x text;echo -e -n ${text})

実行結果

n=10
min=-1.6
max=3.7
med=0.35
ave=0.75
sum=7.5
svar=2.8805
uvar=3.20056
ssd=1.6972
usd=1.78901
skew=0.580921
kurt=-0.629822

【その２：データファイルをパイプで流し込む】

実行

cat q1.txt | (. ./basic_statistics.sh; calc_basic_statistics)

【その３：shスクリプトからの実行】

実行

#!/bin/bash
. ./basic_statistics.sh
IFS=',' read n min max med ave sum svar uvar ssd usd skew kurt text <<< $(calc_basic_statistics q1.txt)
echo "=============================="
printf "n = %d\n" ${n}
printf "min = %G\n" ${min}
printf "max = %G\n" ${max}
printf "med = %G\n" ${med}
printf "ave = %G\n" ${ave}
printf "sum = %G\n" ${sum}
printf "svar = %G\n" ${svar}
printf "uvar = %G\n" ${uvar}
printf "ssd = %G\n" ${ssd}
printf "usd = %G\n" ${usd}
printf "skew = %G\n" ${skew}
printf "kurt = %G\n" ${kurt}
echo "=============================="
printf "${text}"
echo "=============================="
echo -n -e "${text}"

実行結果

==============================
n = 10
min = -1.6
max = 3.7
med = 0.35
ave = 0.75
sum = 7.5
svar = 2.8805
uvar = 3.20056
ssd = 1.6972
usd = 1.78901
skew = 0.580921
kurt = -0.629822
==============================
n=10
min=-1.6
max=3.7
med=0.35
ave=0.75
sum=7.5
svar=2.8805
uvar=3.20056
ssd=1.6972
usd=1.78901
skew=0.580921
kurt=-0.629822
==============================
n=10
min=-1.6
max=3.7
med=0.35
ave=0.75
sum=7.5
svar=2.8805
uvar=3.20056
ssd=1.6972
usd=1.78901
skew=0.580921
kurt=-0.629822

※以下は検算に使用したExcelの分析ツールの結果です。

環境

ホスト　Windows10 COREi7
VM　　　VirtualBox バージョン 5.2.6 r120293 (Qt5.6.2)
　　　　　CentOS Linux release 7.4.1708 (Core)
GNU bash, バージョン 4.2.46(2)-release (x86_64-redhat-linux-gnu)
GNU Awk 4.2.1, API: 2.0 (GNU MPFR 3.1.1, GNU MP 6.0.0)

おわりに

歪度と尖度の計算には不偏標準偏差を使用しています。表計算ソフト(Excel)に合わせたためですが、本当にこれでよいのか自信がありません。