Go to Qiita Advent Calendar 2024 Top
LoginSignup
3
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@EjiriAkira(章 江尻)

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

Last updated at Posted at 2018-03-29

#はじめに
基本統計量を計算する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の分析ツールの結果です。

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

3
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?