先日、簡単にヒストグラムを作成できるサイトを作ったわけですが、内部は R の hist() 関数を参考にさせてもらいました。
ヒストグラムの作成なんて簡単そう、と思われる方もいらっしゃるかもしれませんが、結構複雑です。
今日は、この複雑な hist() 関数を、コアの部分だけ取り出して、簡単に説明しようと思います。
階級数を決める
まず、ヒストグラムを作成するときの問題点として、階級数をどうするかというのが問題になります。
つまり、データをいくつの領域に分割するか、ということです。
ヒストグラムの形状は、階級の分け方によって様々に変わります。
サンプルの数が少ないのに階級を多く分けてしまうと、縦がつぶれたヒストグラムができます。
逆に、階級数が少なすぎると横がつぶれてしまいます。
hist() 関数では、自動的に階級数を決める方法として、
- スタージェスの方法
- フリードマン-ディアコニスの方法
- スコットの方法
の3つが選択できます。
これらはそれぞれ、
nclass.Sturges() nclass.FD() nclass.scott()
という関数によって計算することができます。
例えば、スタージェスの方法を使うと、
x <- rnorm(100) nclass <- nclass.Sturges(x) nclass
[1] 8
となり、100 個の正規乱数に対する階級数(分割数)は 8 と決まるわけです。
階級幅を決める
次に、階級幅を決めます。
まず単純に、データの最小値から最大値までの等差数列を作ってみましょう。
breaks <- seq(min(x), max(x), length = nclass + 1) breaks
[1] -3.27601062 -2.61502035 -1.95403008 -1.29303981 -0.63204954 0.02894073 [7] 0.68993100 1.35092127 2.01191154
この方法では、階級の境目がきりのいい数値にはなりません。
そこで、pretty() という関数を使ってみることにします。
breaks <- pretty(x, n = nclass) breaks
[1] -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5
pretty() は、データをきりのいい階級幅で分けてくれる関数です。*1
pretty() 探訪
ちょっとここで寄り道して、関数 pretty() を探訪してみましょう。
といっても、関数 pretty() の内部は、C言語で書かれています。
そこで、pretty() を単純化した、simple.pretty() という関数を R で書いてみたので、それを見ていきましょう。
simple.pretty <- function(x, n) { lower <- min(x) upper <- max(x) cell <- (upper - lower) / n base <- 10^(floor(log10(cell))) if(cell <= 1.4*base) { unit <- base } else if(cell < 2.8*base) { unit <- 2 * base } else if(cell < 7*base) { unit <- 5 * base } else { unit <- 10 * base } lower.per.unit <- floor(lower / unit) upper.per.unit <- ceiling(upper / unit) lower <- lower.per.unit * unit upper <- upper.per.unit * unit ndiv <- upper.per.unit - lower.per.unit s <- seq.int(lower, upper, length.out = ndiv + 1) return(s) }
まず、lower, upper にデータの最小値、最大値をそれぞれ入れます。
それらの差を階級数で割り、階級幅の基本となる値 cell を出します。
求めたい階級幅は、この cell に近い値で、きりのいい数値にすればいいということです。
そのために、この cell の基本単位 base を求めます。
base <- 10^(floor(log10(cell)))
下記表のように base は10進数表現での基本単位になります。
| cell | base |
|---|---|
| 0.06 | 0.01 |
| 0.15 | 0.1 |
| 3.14 | 1 |
| 1024 | 1000 |
次は、cell が基本単位の何倍ぐらいの範囲にあるかで unit を求めます。
| cell が baseの | unit |
|---|---|
| 1倍 〜 1.4倍 | base |
| 1.4倍 〜 2.8倍 | 2 * base |
| 2.8倍 〜 7倍 | 5 * base |
| 7倍 〜 10倍 | 10 * base |
この unit が求めたかった「きりのいい階級幅」になります。
あとは、lower, upper を、この unit で表現できる値に直し、unit を幅とした等差数列を返せばいいだけです。
breaks <- simple.pretty(x, n = nclass) breaks
[1] -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5
きりのいい数値が返ってきました。
階級に入るデータ数のカウント
階級数、階級幅が決まったので、次は階級に入るデータ数をカウントします。
hist() 関数内では、この処理も C言語で書かれた関数を使っています。
というわけで、ここも同じく単純化した simple.bincount() という関数を R で書いてみましたので、それを見てみましょう。
simple.bincount <- function(x, breaks) { nx <- length(x) nbreaks <- length(breaks) counts <- integer(nbreaks - 1) for(i in 1:nx) { lo <- 1 hi <- nbreaks if(breaks[lo] <= x[i] && x[i] <= breaks[hi]) { while(hi - lo >= 2) { new <- (hi + lo) %/% 2 if(x[i] > breaks[new]) lo <- new else hi <- new } counts[lo] <- counts[lo] + 1 } } return(counts) }
データの入る階級を求めるために二分探索をしています。
二分探索については ここ(Wikipedia) をご参照ください。
データが入るかどうかを、単純に階級を一つ一つ調べていくよりも、二分探索した方が早く見つかります。
これを実行すると次のようになります。
counts <- simple.bincount(x, breaks) counts
[1] 1 2 1 1 8 16 24 20 14 9 3 1
それぞれの階級に、データがいくつ入るかが求まりました。
グラフィックス
あとは、上で求めたデータをもとにヒストグラムを描けばいいだけです。
r <- structure(list(breaks = breaks, counts = counts), class = "histogram") plot(r)
R のグラフィックスについてはあまり詳しくないので説明は勘弁して下さい。
ちなみに、簡単にヒストグラムを作成できるサイト では、ヒストグラム描画には Google Chart Tools を使っています。
まとめ
以上をまとめると、次のような関数が出来上がります。
my.hist <- function(x) { nclass <- nclass.Sturges(x) # スタージェス # nclass <- nclass.FD(x) # フリードマン-ディアコニス # nclass <- nclass.scott(x) # スコット breaks <- simple.pretty(x, n = nclass) counts <- simple.bincount(x, breaks) r <- structure(list(breaks = breaks, counts = counts), class = "histogram") plot(r, main = paste("Histogram of", substitute(x))) }
実行してみます。
x <- rnorm(100) my.hist(x)
ちゃんとヒストグラムが描かれましたね。
以上です。
