古い記事
ランダムジャンプ
新しい記事
スタージェスの公式(Sturges' formula)とは、ヒストグラムを書くときに、観測値の数 n に対して、柱の数(ビンの数、階級数)k の最適値、とは言わないまでも、そこそこ有用な値を求める手段の一つ。

histgram

公式と計算例


スタージェスの公式:


Perl による計算例。
% perl -MPOSIX -le '$n=shift; $k=1+log($n)/log(2); print ceil$k' 10
5
% perl -MPOSIX -le '$n=shift; $k=1+log($n)/log(2); print ceil$k' 1000
11
% perl -MPOSIX -le '$n=shift; $k=1+log($n)/log(2); print ceil$k' 100000
18
% perl -MPOSIX -le '$n=shift; $k=1+log($n)/log(2); print ceil$k' 10000000
25

つまり、観測値が1000個ある場合は、11本の柱(ビン、階級)によるヒストグラムを作ればよいのである。

なお、柱の幅hは以下の式で求める。


ヒストグラム出力スクリプト


以下、スタージェスの公式を利用した、テキスト形式のヒストグラム出力スクリプト(Perl)。入力は一行一観測値。階級数はスタージェスで自動決定。

■ソース(hist.pl):
#!/usr/bin/perl
use strict;
use warnings;
use POSIX 'ceil';
use List::Util qw(max min sum);

my $bar_width = 60;

### データ読み込み
my @vs;
while (<>) {
    next if not /^(-?[\d\.e-]+)/;
    push @vs, $1;
}

### 値補正:平均値が1以下だったら最小値を1以上にすべく
### すべての値を10^X倍する。
if (sum(@vs)/@vs <= 1) {
    my $l = abs(int(log(min(grep {$_} @vs))/log(10)))+1;
    @vs = map {$_*(10**$l)} @vs;
}

my ($max, $min) = (max(@vs), min(@vs));
my $k = ceil(1 + log(@vs)/log(2)); # 階級値
my $h = ceil(($max - $min) / $k); # 階級幅

### 各ビンの頻度
my @freq;
foreach my $t (sort {$a <=> $b} @vs) {
    my $a = ceil(($t-$min)/$h) || 1;
    $freq[$a]++;
}

### ビンの値の範囲(表示用)
my @kgr = ($min);
for (my $i = 1; $i <= $k; $i++) {
    $kgr[$i] = $kgr[$i-1] + $h;
}

my $max_f = max(map {$_||0} @freq);
my $keta_f = length($max_f);
my $keta_m = length(int($kgr[-1]));
for (my $i = 1; $i <= $k; $i++) {
    my $f = $freq[$i]||0;
    my $bar = "*" x ($max_f <= $bar_width ? $f : $f/($max_f/$bar_width));
    printf "%02d %${keta_f}d ( %${keta_m}d-%-${keta_m}d ) %s\n",
    $i, $f, $kgr[$i-1], $kgr[$i], $bar;
}
※平均値が1以下のときに、最小値が1以上10以下になるように全体の値を補正(10のN乗をかけている)している。

■実行例:
% cat hist.test
1
2
3
4
5
20
% head -5 hist.test | ./hist.pl 
01 2 ( 1-2 ) **
02 1 ( 2-3 ) *
03 1 ( 3-4 ) *
04 1 ( 4-5 ) *
% head -6 hist.test | ./hist.pl 
01 5 (  1-6  ) *****
02 0 (  6-11 ) 
03 0 ( 11-16 ) 
04 1 ( 16-21 ) *

ランダム:
% perl -le 'print int(rand(100)) for(1..10)' | ./hist.pl 
01 3 (  8-26 ) ***
02 1 ( 26-44 ) *
03 2 ( 44-62 ) **
04 2 ( 62-80 ) **
05 2 ( 80-98 ) **
% perl -le 'print int(rand(100)) for(1..100)' | ./hist.pl 
01  6 (   0-13  ) ******
02 20 (  13-26  ) ********************
03  7 (  26-39  ) *******
04 13 (  39-52  ) *************
05 20 (  52-65  ) ********************
06  8 (  65-78  ) ********
07 13 (  78-91  ) *************
08 13 (  91-104 ) *************
% perl -le 'print int(rand(100)) for(1..1000)' | ./hist.pl
01 100 (  0-9  ) ************************************************************
02  90 (  9-18 ) ******************************************************
03  98 ( 18-27 ) **********************************************************
04  91 ( 27-36 ) ******************************************************
05  84 ( 36-45 ) **************************************************
06  95 ( 45-54 ) *********************************************************
07  81 ( 54-63 ) ************************************************
08 100 ( 63-72 ) ************************************************************
09  87 ( 72-81 ) ****************************************************
10  78 ( 81-90 ) **********************************************
11  96 ( 90-99 ) *********************************************************

0以上1未満のとき(補正済み):
% perl -le 'print int(rand(1)*100)/100 for(1..50)' | ./hist.pl
01  8 (  0-14 ) ********
02 10 ( 14-28 ) **********
03  5 ( 28-42 ) *****
04  8 ( 42-56 ) ********
05  6 ( 56-70 ) ******
06  8 ( 70-84 ) ********
07  5 ( 84-98 ) *****

とりあえず深く考えないで数字をつっこめばグラフを作ってくれるので結構便利かと思われ。

参考文献

- Wikipedia:ヒストグラム
- 東京大学教養学部統計学教室 (編集), "統計学入門", 東京大学出版会, 1991. (p.22)