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

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

#**基本統計量**

| 統計量 | 記号 | 変数 | 説明| 表計算ソフト(Excel)の関数 |
|:--|:--|:--|:--|:--|
| データ  | -  | ds[]  | データ(標本)  |-|
| データ数  | n  | n  | データ数(標本数)  |-|
| 最小値  | -  |  min  | 最小値  |MIN()|
| 最大値  | -  |  max  | 最大値  |MAX()|
| 中央値  | -  |  med |ソートされたデータの中央値<BR>データ数が偶数の場合は中央の２値の平均|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スクリプト】**

~~~bash: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をエスケープしています。
<PRE>n=値\\\\nmin=値\\\\nmax=値\\\\nmed=値\\\\nave=値\\\\nsum=値\\\\nsvar=値\\\\nuvar=値\\\\nssd=値\\\\nusd=値\\\\nskew=値\\\\nkurt=値\\\\n</PRE>


#**実行例**

#####**【データファイル】**

~~~bash:q1.txt
0.7
-1.6
-0.2
-1.2
-0.1
3.4
3.7
0.8
0.0
2.0
~~~


#####**【その１：データファイルの直接指定】**
~~~bash:実行

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

~~~
~~~bash:実行結果

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
~~~

出力の整形

~~~bash:実行
(. ./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})
~~~

~~~bash:実行結果
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
~~~

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


~~~bash:実行
cat q1.txt | (. ./basic_statistics.sh; calc_basic_statistics)
~~~

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


~~~bash:実行
#!/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}"

~~~

~~~bash:実行結果
==============================
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の分析ツールの結果です。

![q1.png](https://qiita-image-store.s3.amazonaws.com/0/234070/f7c639a2-42a8-f03d-bd6f-5f3676c1381a.png)


#**環境**

ホスト　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)に合わせたためですが、本当にこれでよいのか自信がありません。