はじめに

基本統計量を計算するbcコマンドのスクリプトです。

スクリプトと出力形式

【bcスクリプト】

basic_statistics.bc

define set_scale(x, s) {
    auto old_scale
    old_scale = scale
    scale = s
    x /= 1  #scale 設定
    scale = old_scale
    return (x)
}
define round(d, s) {
    auto temp;
    if(scale(d) < s) return (set_scale(d, s))
    temp = (d - set_scale(d, s)) * (10^(s+1))
    if(temp > 5) d += (1/(10^s))
    return (set_scale(d,s))
}
define power(d, i) {
    return(e(i*l(d)))
}

scale = 100
n = read()
for(i = 0; i < n; i++){
    ds[i] = read()
}
sum = 0
for(i = 0; i < n; i++){
    d = ds[i]
    sum += d
    if( i == 0) {
        max = d
        min = d
    } else {
        if(d > max){
            max = d;
        }
        if(d < min){
            min = d;
        }
    }
}
old_scale = scale
if( (n % 2) == 0) {
    scale = 0
    m1 = ds[n/2]
    m2 = ds[(n/2)+1]
    scale = old_scale
    med = (m1 + m2) / 2
} else {
    scale = 0
    med = ds[n/2]
    scale = old_scale
}
ave = sum / n
sum_dev2 = 0
sum_dev3 = 0
sum_dev4 = 0
for(i = 0; i < n; i++){
    dev = ds[i] - ave
    devs[i] = dev
    sum_dev2 += (dev * dev)
    sum_dev3 += (dev * dev * dev)
    sum_dev4 += (dev * dev * dev * dev)
}
svar = sum_dev2 / n
ssd = sqrt(svar)
uvar = sum_dev2 / (n - 1)
usd = sqrt(uvar)

skew = (n /((n -1 ) * (n - 2))) * sum_dev3 / (usd^3)
kurt = ((n * (n + 1)) /((n -1 ) * (n - 2) * (n - 3))) * (sum_dev4 / (usd^4)) - ((3 * (n - 1)^2) / ((n - 2) * (n - 3)))

print "n = "; n
print "min = "; min
print "max = "; max
print "med = "; round(med, 20)
print "ave = "; round(ave, 20)
print "sum = "; round(sum, 20)
print "svar = "; round(svar, 20)
print "uvar = "; round(uvar, 20)
print "ssd = "; round(ssd, 20)
print "usd = "; round(usd, 20)
print "skew = "; round(skew, 20)
print "kurt = "; round(kurt, 20)

【出力形式】

n = データ数
min = 最小値
max = 最大値
med = 中央値
ave = 平均値
sum = データの総和
svar = 標本分散
uvar = 不偏分散
ssd = 標本標準偏差
usd = 不偏標準偏差    
skew = 歪度
kurt = 尖度

実行例

【データファイル】

q11.txt

1
2 2 2
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
4 4
5

データは空白、タブ、改行で区切る。

【実行方法】

実行例

((cat q11.txt | wc -w); sort -g q11.txt) | bc -l -q basic_statistics.bc

データ数、データを空白、タブ、改行で区切って、bcにわたします。
中央値を計算する都合でデータは(数値として)ソートする必要があります。

実行結果

n = 30
min = 1
max = 5
med = 3.00000000000000000000
ave = 2.96666666666666666667
sum = 89.00000000000000000000
svar = .43222222222222222222
uvar = .44712643678160919540
ssd = .65743609744386733227
usd = .66867513545937177714
skew = .03679696392430356743
kurt = 4.50339722463221514906

環境

ホスト　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)
bc 1.06.95

おわりに

「gawkで基本統計量を計算する」の検算用に作りました。
bcを選んだのは以下の理由によります。

多倍長演算(Multiple Precision Arithmetic)をサポートしている
できるだけ標準的なもの（Unix/Linux）

データの読み込みにread()を使っています。bcはスクリプトの読み込みもread()も標準入力を使用しているので、shスクリプトとして作成するのが難しいです。数が少なければ環境変数で渡すのも手ですが、配列に渡すのは骨です。
以上

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