基本的なグラフとして理解しておくべきは、
棒グラフ、円グラフ、折れ線グラフの三種類であるが、
グラフには様々なものがあり、状況に合わせて適切に使うことが重要。
さらに、本ブログにおける一連の投稿記事では、上記の三種類に加えて、
ヒストグラム、箱ひげ図、散布図についても解説する。
今回は、ヒストグラムについて整理する。
ヒストグラムは、どうも、誤解の多いグラフである。
見た目が棒グラフに似ているので、棒グラフと混同している人も多い。
しかしながら、ヒストグラムと棒グラフは似て非なる物である。
少なくとも、Excel 棒グラフでヒストグラムの代用はできない。
出来る場合もあるが、色々と機能が足りないのである。
とりあえず、いつものように...
# Windows と Linux の人は次のコマンド
X <- read.table("clipboard", sep=",")
# Mac の人は次のコマンド
X <- read.table(pipe("pbpaste") , sep=",")
データを読み込む。
Aji,Saba,Buri,Hirame,Karei,Kasago,Mebaru,Kawahagi,Tai,Isaki,Ika,Tako
Hokkaido,0,2,18,12,1210,7,452,0,0,0,660,3
Aomori,0,0,52,14,111,2,66,0,77,0,25,7
Iwate,0,3,2,0,14,0,230,0,0,0,12,2
Miyagi,1,75,25,33,2140,0,116,0,0,0,1,1
Akita,21,1,33,4,0,1,76,0,55,0,4,5
Yamagata,110,0,19,6,0,0,14,0,24,0,7,0
Fukushima,0,0,0,0,72,0,60,0,0,0,0,0
Iabaraki,14,57,39,45,36,20,112,1,79,0,109,46
Chiba,381,96,92,100,1,29,65,34,154,135,328,2
Tokyo,38,6,0,0,3,1,0,2,1,0,0,2
Kanagawa,683,185,823,1,1,32,56,29,168,21,256,11
Nigata,194,25,56,5,2,1,96,1,161,0,483,0
Toyama,1,0,21,0,0,0,0,0,3,0,9,0
Ishikawa,88,17,25,6,0,5,36,5,32,0,0,0
Fukui,160,5,307,4,0,6,3,23,170,1,157,0
Shizuoka,56,109,58,11,0,65,28,5,28,205,8,0
Aichi,20,67,98,1,1,20,8,3,150,22,5,1
Mie,19,15,106,0,0,18,13,1,210,77,61,0
Kyoto,32,5,81,1,5,15,6,3,91,4,9,0
Osaka,32,0,2,0,0,5,3,0,14,0,4,0
Hyogo,93,23,81,9,0,125,68,3,20,0,15,49
Wakayama,22,80,99,6,0,6,15,4,116,121,59,0
Tottori,16,0,3,0,0,0,2,0,29,14,41,0
Shimane,24,3,15,7,0,11,6,17,25,16,155,0
Okayama,1,0,0,0,0,34,2,0,40,0,1,11
Hiroshima,64,7,1,0,1,45,48,1,17,0,10,24
Yamaguchi,274,27,366,49,3,52,94,66,308,16,1,0
Tokushima,16,10,44,5,0,21,14,1,29,6,96,0
Kagawa,0,0,0,0,0,0,0,0,2,0,0,0
Ehime,58,0,9,0,0,25,9,5,55,2,0,2
Kochi,0,2,1,0,0,0,0,0,2,7,0,0
Fukuoka,0,0,0,0,0,0,0,0,0,48,82,0
Nagasaki,0,0,384,22,0,131,0,0,306,155,0,0
Kumamoto,55,0,3,0,0,53,0,0,123,2,0,0
Oita,544,97,4,17,0,62,2,217,75,12,32,0
Miyazaki,49,27,3,0,0,0,0,1,35,63,0,0
Kagoshima,109,52,130,1,0,56,0,1,140,23,7,25
Okinawa,20,0,21,0,0,0,0,4,7,3,0,0
今回のデータは、農林水産省の「都道府県別魚種別遊漁採捕量」のデータ。
私の大好きな魚介のデータ満載！今度、じっくりとデータを見よう。
今回のサンプル、元のデータの属性（変数）の数が多すぎたので適当に編集してある。
編集の仕方に問題があるので、間違っても、このデータを何かに使わないように！
さて、とりあえず、ヒストグラムを書いてみましょう。説明は後で。
Rでは、hist()関数を使ってヒストグラムを描く。描き方は以下の通り。
# 鯵のヒストグラム
hist(X$Aji)
この図、棒グラフに良く似ているように思える。
一体、何が異なっているのか？自分で考えてみるのも重要。
棒グラフの話では、どのように説明であったか、今一度思い出してみる。
ある変数について、ある対象の数量を可視化するのが棒グラフである、
と説明した。そのように説明したハズである。
では、この図では何が表現されているか、よく観察してみる。
まず、y軸ラベルには、「Frequency」と書かれている。
日本語では、「頻度」と訳される。
次に、x軸の目盛りを見てみる。何やら数字が書いてある。
これは、ある変数における「値の幅」を表してる。
棒グラフでは、ある変数のある対象の数量の大きさを棒の長さで表現したが、
ヒストグラムでは、ある「値の幅」の中に含まれる「対象の個数」を、
「帯の面積」で表しているのである。
多くの人は、この説明で理解できるハズなのだが。
一応、念の為に棒グラフとヒストグラムを比較してみる。
# グラフを並べて描くためのオマジナイ
par(mfrow=c(1,2))
# まずは、棒グラフ
barplot(X$Aji,main="Barplot for Aji", names=rownames(X), las=2, ylab="Number of fish catches (tones)")
このように、データの全体を把握するために用いるのがヒストグラムである。
あまり、有用性を理解できない人もいるかもしれないが、
実は、非常に良く利用する。棒グラフよりも。
例えば、大学の成績を付ける際には、ヒストグラムを最初に描いて、
成績がどのように分布しているのかを観察し、合格ラインを考えたりする。
また、ある集落における立地環境を調べるために、
集落周辺の標高や土地利用を集計し、集落の立地特性を調べるときにも利用する。
ちょっと、特殊な使い方ではあるが。
とにかく、利用価値の高いのがヒストグラムである。
ヒストグラムについての基本的なここまで。
ここからは、ヒストグラムの描き方について話。
ヒストグラムを描く際に最初に気になるのは、
一体、何個の区切りを設けるべきかという問題である。
この区切りのことを、「階級（break）」と呼び、
その階級の値を「階級値（break point）」と呼ぶ。
何も考えずに、コンピュータの既定値を使っている人も居るが、それはNG。
ちゃんと、分け方というのも決まっている。
とは言え、階級数を決めるための理想的な方法は存在しない。
Rでは、階級数を決定するための方法がいくつか存在し、
これらを使用する以外にも、自分で階級を設定することもできる。
状況に合わせて、説明しやすいものを選ぶ。この辺は適当に。
何も説明できないのはダメであるが、説明できるのであればどの方法でも良い。
ただし、データの性質と分類方法の相性は存在する。
さて、階級数を決定する方法の中で、最もシンプルな方法は「平方根選択」。
以下の方法で選択する。説明は不要であろう。
Rでの実行方法は以下の通り。
(k <- round(sqrt(nrow(X))))
今回のデータの場合には、「6」が適切な結果となるようだ。
> (k <- round(sqrt(nrow(X))))
[1] 6
階級数は、整数で求める必要があるので、round()関数でシンプルに四捨五入。
丸め込みには他の方法も。それは、次の説明で。
ところで、この方法、数が増えるほど、階級数が大きくなってしまう。
例えば、10000の対象があると、「階級数は100」にもなってしまう。
データの全体を俯瞰するには、少々、大き過ぎはしなか？
平方根選択は、非常にシンプルな方法であるが、
最も一般的なものに、「スタージェスの公式」と呼ばれるものがある。
hist()関数の既定値。何も指定しない場合は、この公式が適用される。
なお、一般的ということで、この方法が正しいという意味ではない！
とりあえず、スタージェスの公式。
括弧の形に注目！大カッコではない。何だこれは？
大カッコに見えるが、括弧の上だけ内側に向いている。これ、「天井関数」の記号。
小難しい定義では、ある実数以上の最小の整数と定義される。
要するに、四捨五入は、0.5以上なら繰り上がり、0.4以下なら繰り下がるのだが、
天井関数では、小数点以下が0以上なら、とにかく繰り上がる。
例えば、「2.1」という数字があれば「3」となる。全部、繰り上がる。
ちなみに、「床関数」というのもあって、これは天井関数の逆。全て繰り下がる。
なお、床関数の記号は、括弧の下だけ内側に向く。
括弧の中は、二進対数となっている。Nは対象の数（総数）。
対数の話をし始めると終わらないのでいずれ。
とにかく、Rでは以下のように計算する。
# ちゃんと公式通りに計算する場合
(k <- ceiling(log2(nrow(X))+1))
ここで、ceiling()関数は、天井関数のこと。訳そのまま。
log2()関数は、二進対数を計算するための関数。
なお、Rには、nclass.Sturges()関数というものがあり、通常はそれを用いる。
# 楽をしてRのStruges()関数を使う場合
(k <- nclass.Sturges(X$Aji))
もちろん、結果は同じで「7」となる。以下は実行結果。
> (k <- ceiling(log2(nrow(X))+1))
[1] 7
> (k <- nclass.Sturges(X$Aji))
[1] 7
スタージェスの公式は、非常に良く利用されるのであるが、
データ数が200を超えるような大きなデータには不向きであり、
また、データが釣り鐘状に分布（正規分布）していない場合にも向かない。
次に、「スコットの公式」。
以下の公式で計算される。これも、それほど難しくない。
右辺の方は、大丈夫であろうか？
分子は、ある変数の最大値から最小値を引いたもの。
分母の「σ（シグマ）」は標準偏差で、Nは対象の数（総数）。
これが解っていれば、この計算は簡単。Rでは以下のように計算する。
# ちゃんと公式通りに計算する場合。
(k <- ceiling((max(X$Aji)-min(X$Aji))/(3.5*sd(X$Aji)*(nrow(X)^(-1/3)))))
nclass.scott() 関数というのもある。通常はこれを使う。
# 楽をしてR のscott()関数を使う場合
(k <- nclass.scott(X$Aji))
以下が実行結果。当然、結果は同じ。「5」階級を推奨する。
> # ちゃんと計算する場合。
> (k <- ceiling((max(X$Aji)-min(X$Aji))/(3.5*sd(X$Aji)*(nrow(X)^(-1/3)))))
[1] 5
> # R のscott()関数を使う場合
> (k <- nclass.scott(X$Aji))
[1] 5
スコットの公式は、スタージェスの公式と同様に、
理論的には釣り鐘状の分布（正規分布）を前提としているが、
実用的には、スタージェスの公式と比較してその制約の度合いが小さい。
階級数を決定する方法は、他にもある。
今度は、「フリードマン＝ダイアコニスの公式（FD法）」。
ふむふむ。
これは、何のことであろうか？
実は、中央値の話を見れば理解できるハズ。そうそう、分位数の話。
中央値の場合は、全データの半分番目であり、二分位に相当するが、
ここでは、Qの添字の分母が「4」となっているので、
全データを四等分した値、四分位。
つまり、3/4番目の値から、1/4番目の値を引いている。
ちにみに、この範囲のことをIQR（Interquantile Range）と呼ぶ。
日本語では「四分位範囲」と呼ぶ。頻繁に出てくる単語。
この範囲には、中央値を中心にデータの50%が収まる。
四分位は、データを観察する上で、最も基本的な指標であり、頻繁に登場する。
この話は、次の回で詳しく述べる。
さて、これも、Rで実際に計算してみる。
# ちゃんと計算する場合。
(k=ceiling((max(X$Aji)-min(X$Aji))/(2*(quantile(X$Aji)[[4]]-quantile(X$Aji)[[2]])*(nrow(X)^(-1/3)))))
四分位範囲は、quantile()関数を使って求めている。
この関数は少々特殊な形式で値を扱っているので大カッコ二つで呼び出している。
それ以外は、特に説明することもあるまい。この方法に関してもRの関数がある。
> # R のFD()関数を使う場合
(nclass.FD(X$Aji))
もちろん結果は同じ。最初の3つの方法比べて数が大きい。「15」を推奨。
> (k=ceiling((max(X$Aji)-min(X$Aji))/(2*(quantile(X$Aji)[[4]]-quantile(X$Aji)[[2]])*(nrow(X)^(-1/3)))))
[1] 15
> (nclass.FD(X$Aji))
[1] 15
エラーが発生するが、以下のコマンドを実行する。
エラーの原因は後で解説する。
つまりは、どういうことであるか？
何やら、先程と様子が異なる。ここで、y軸ラベルとその目盛りに注目。
軸ラベルが「density（密度）」となっていて、値が非常に小さくなっている。
棒グラフというのは、棒の長さで値の大きさを表していたが、
ヒストグラムでは、棒の長さよりも、面積によって値の大きさを表現する。
したがたって、階級幅を個別に変えて、頻度を選択すると、
ヒストグラムの帯の面積比と対応する値の間に矛盾が生じてしまう。
これがエラーの原因であった。
なるほど。確かに納得できる。頻度で描いたヒストグラムでは、
直感的に、100t以上の帯の高さがどうも不自然であった。
このように、ヒストグラムが面積によって表されていることを意識するのは重要。
では、この値は一体どのような値となっているのか？
何らかの方法で基準化されているのは確かのようである。
とりあえず、次の式を見てみることにする。
この式において
は、階級値であり、今回のデータでは、
0, 20, 40, 60, 80, 100, 200, 400, 600, 700 の9つの値から成るベクトルが該当する。
添字の意味は...大丈夫であろう。「i+1番目」と「i番目」のこと。
は、その階級に該当する密度である。この式を見て、意味が解ったらセンスが良い。
実は、この式は、面積の和を表していて、その和が「1」となるようになっている。
うん？疑う人も居るかもしれない。
数式に慣れていない人には辛いかもしれない。
とりあえず、四角形の面積の計算はどようなものであったか？
これは大丈夫であろう。
面積＝縦×横
ふむふむ。これは納得できる。では、このヒストグラムにおける帯の面積はどうなるか？
とりあえず、一つ目の帯だけで考えることにする。
横・・・x軸なので、階級値の二番目の「20」から一番目の「0」を引いて「20」
縦・・・y軸なので、おおよそ「0.22」くらい。今は、グラフを見て目測。
したがって、一つ目の帯の面積は、おおよそ「20×0.22=4.4」くらいだろうか？
さて、これは少し大雑把だったかもしれない。ちゃんとした計算をしてみる。
次のコマンドでは、グラフを描画しない。慌ててはいけない。
# hist()関数の結果を、変数hに格納する。図は出さない。
h <- hist(X$Aji, main=title, xlab="number of fish catches", breaks=breaks, plot=FALSE)
# 変数hに格納されている中身を確認する。
str(h)
今回は、ヒストグラムの結果を変数「h」に格納にしている。
plot = FALSE は、グラフを描画しないためのパラメータ。
二つ目のstr()関数は、変数の中身を確認するための関数。
実行結果は次のようになる。
> # 変数hに格納されている中身を確認する。
> str(h)
List of 7
$ breaks : num [1:10] 0 20 40 60 80 100 200 400 600 700
$ counts : int [1:9] 17 6 4 1 2 4 2 1 1
$ intensities: num [1:9] 0.02237 0.00789 0.00526 0.00132 0.00263 ...
$ density : num [1:9] 0.02237 0.00789 0.00526 0.00132 0.00263 ...
$ mids : num [1:9] 10 30 50 70 90 150 300 500 650
$ xname : chr "X$Aji"
$ equidist : logi FALSE
- attr(*, "class")= chr "histogram"
とりあえず、ここで必要となるのは、breaks と density の値。他は、今回は無視。
変数「h」の中身を呼び出す方法は次の通り。
# 階級の値をベクトルとして取り出す
h$breaks
# 密度の値をベクトルで取り出す
h$density
実行結果は次の通り。
> h$breaks
[1] 0 20 40 60 80 100 200 400 600 700
> h$density
[1] 0.0223684211 0.0078947368 0.0052631579 0.0013157895 0.0026315789
[6] 0.0010526316 0.0002631579 0.0001315789 0.0002631579
ふむふむ。先程は、目測であったが、
最初の帯の密度の値は「0.0223684211」である。
なるほど、我ながら、良い目測であった。
では、面積の合計がちゃんと「1」になるかを確認する。
# まずは、結果を初期化する。
s2 = 0
# For文を回して、各帯の面積を足していく。
for(i in seq(1:(length(h$breaks)-1))){
# i番目までの総面積＝i-1番目までの総面積+（i番目の帯の高さ×i番目の帯の横幅）
s2 = s2 + h$density[i] * (h$breaks[i+1]-h$breaks[i])
# i番目までの総面積
print(s2)
}
実行結果は以下の通り。面積の計算をFor文で繰り返してて、結果を足しているだけ。
print()関数は、毎回の計算結果を逐次表示させるためのもの。無くても良い。
結果は、次のようになる。
> # まずは、結果を初期化する。
> s2 = 0
>
> # For文を回して、各帯の面積を足していく。
> for(i in seq(1:(length(h$breaks)-1))){
+ # i番目までの総面積＝i-1番目までの総面積+（i番目の帯の高さ×i番目の帯の横幅）
+ s2 = s2 + h$density[i] * (h$breaks[i+1]-h$breaks[i])
+
+ # i番目までの総面積
+ print(s2)
+ }
[1] 0.4473684
[1] 0.6052632
[1] 0.7105263
[1] 0.7368421
[1] 0.7894737
[1] 0.8947368
[1] 0.9473684
[1] 0.9736842
[1] 1
少々、回りくどい処理だったかもしれないが、
最後番目（正確には9番目）の総面積の結果が、ちゃんと「1」になっている。
なるほど、確かに面積となっている。
さて、今回も、すっかり長い話になってしまったが、
たかが、ヒストグラムに詳しい説明をしてくれる教科書は意外に少ない。
特に、階級数の決定方法は、スタージェスの公式で済ませるものも多い。
また、ヒストグラムは、多くの場合は、頻度で描いているが、密度で描く場合もある。
面積という観点から考えると、実は、密度で描いた方が良いのである。
基準化されていて、ヒスグラムとして直感的に理解しやすい。
最後に、スコットの公式を当てはめ、密度で描いた図を示す。
# 最終的に描かれたヒストグラム
# 並べて描くためのオマジナイ。全部で12個
par(mfrow=c(4,3))
for(i in seq(1:ncol(X))){
# タイトルラベルの調整
title <- paste("Histogram of fish catches of ", colnames(X)[i], sep="")
title <- paste(title, "Scott's normal reference rule", sep="\n")
# ヒストグラムを描く。色は、薄い青色で
hist(X[,i], main=title, xlab="number of fish catches",
breaks=nclass.scott(X[,i]), freq=FALSE, col="lightblue")
}
棒グラフ、円グラフ、折れ線グラフの三種類であるが、
グラフには様々なものがあり、状況に合わせて適切に使うことが重要。
さらに、本ブログにおける一連の投稿記事では、上記の三種類に加えて、
ヒストグラム、箱ひげ図、散布図についても解説する。
今回は、ヒストグラムについて整理する。
ヒストグラムは、どうも、誤解の多いグラフである。
見た目が棒グラフに似ているので、棒グラフと混同している人も多い。
しかしながら、ヒストグラムと棒グラフは似て非なる物である。
少なくとも、Excel 棒グラフでヒストグラムの代用はできない。
出来る場合もあるが、色々と機能が足りないのである。
とりあえず、いつものように...
# Windows と Linux の人は次のコマンド
X <- read.table("clipboard", sep=",")
# Mac の人は次のコマンド
X <- read.table(pipe("pbpaste") , sep=",")
データを読み込む。
Aji,Saba,Buri,Hirame,Karei,Kasago,Mebaru,Kawahagi,Tai,Isaki,Ika,Tako
Hokkaido,0,2,18,12,1210,7,452,0,0,0,660,3
Aomori,0,0,52,14,111,2,66,0,77,0,25,7
Iwate,0,3,2,0,14,0,230,0,0,0,12,2
Miyagi,1,75,25,33,2140,0,116,0,0,0,1,1
Akita,21,1,33,4,0,1,76,0,55,0,4,5
Yamagata,110,0,19,6,0,0,14,0,24,0,7,0
Fukushima,0,0,0,0,72,0,60,0,0,0,0,0
Iabaraki,14,57,39,45,36,20,112,1,79,0,109,46
Chiba,381,96,92,100,1,29,65,34,154,135,328,2
Tokyo,38,6,0,0,3,1,0,2,1,0,0,2
Kanagawa,683,185,823,1,1,32,56,29,168,21,256,11
Nigata,194,25,56,5,2,1,96,1,161,0,483,0
Toyama,1,0,21,0,0,0,0,0,3,0,9,0
Ishikawa,88,17,25,6,0,5,36,5,32,0,0,0
Fukui,160,5,307,4,0,6,3,23,170,1,157,0
Shizuoka,56,109,58,11,0,65,28,5,28,205,8,0
Aichi,20,67,98,1,1,20,8,3,150,22,5,1
Mie,19,15,106,0,0,18,13,1,210,77,61,0
Kyoto,32,5,81,1,5,15,6,3,91,4,9,0
Osaka,32,0,2,0,0,5,3,0,14,0,4,0
Hyogo,93,23,81,9,0,125,68,3,20,0,15,49
Wakayama,22,80,99,6,0,6,15,4,116,121,59,0
Tottori,16,0,3,0,0,0,2,0,29,14,41,0
Shimane,24,3,15,7,0,11,6,17,25,16,155,0
Okayama,1,0,0,0,0,34,2,0,40,0,1,11
Hiroshima,64,7,1,0,1,45,48,1,17,0,10,24
Yamaguchi,274,27,366,49,3,52,94,66,308,16,1,0
Tokushima,16,10,44,5,0,21,14,1,29,6,96,0
Kagawa,0,0,0,0,0,0,0,0,2,0,0,0
Ehime,58,0,9,0,0,25,9,5,55,2,0,2
Kochi,0,2,1,0,0,0,0,0,2,7,0,0
Fukuoka,0,0,0,0,0,0,0,0,0,48,82,0
Nagasaki,0,0,384,22,0,131,0,0,306,155,0,0
Kumamoto,55,0,3,0,0,53,0,0,123,2,0,0
Oita,544,97,4,17,0,62,2,217,75,12,32,0
Miyazaki,49,27,3,0,0,0,0,1,35,63,0,0
Kagoshima,109,52,130,1,0,56,0,1,140,23,7,25
Okinawa,20,0,21,0,0,0,0,4,7,3,0,0
今回のデータは、農林水産省の「都道府県別魚種別遊漁採捕量」のデータ。
私の大好きな魚介のデータ満載！今度、じっくりとデータを見よう。
今回のサンプル、元のデータの属性（変数）の数が多すぎたので適当に編集してある。
編集の仕方に問題があるので、間違っても、このデータを何かに使わないように！
さて、とりあえず、ヒストグラムを書いてみましょう。説明は後で。
Rでは、hist()関数を使ってヒストグラムを描く。描き方は以下の通り。
hist(X$Aji)
この図、棒グラフに良く似ているように思える。
一体、何が異なっているのか？自分で考えてみるのも重要。
棒グラフの話では、どのように説明であったか、今一度思い出してみる。
ある変数について、ある対象の数量を可視化するのが棒グラフである、
と説明した。そのように説明したハズである。
では、この図では何が表現されているか、よく観察してみる。
まず、y軸ラベルには、「Frequency」と書かれている。
日本語では、「頻度」と訳される。
次に、x軸の目盛りを見てみる。何やら数字が書いてある。
これは、ある変数における「値の幅」を表してる。
棒グラフでは、ある変数のある対象の数量の大きさを棒の長さで表現したが、
ヒストグラムでは、ある「値の幅」の中に含まれる「対象の個数」を、
「帯の面積」で表しているのである。
多くの人は、この説明で理解できるハズなのだが。
一応、念の為に棒グラフとヒストグラムを比較してみる。
par(mfrow=c(1,2))
# まずは、棒グラフ
barplot(X$Aji,main="Barplot for Aji", names=rownames(X), las=2, ylab="Number of fish catches (tones)")
# 次に、ヒストグラム
hist(X$Aji, main="Histogram of fish catches of Aji", xlab="number of fish catches")
今回は、ヒストグラムのパラメータで、mainとxlabが指定されているが、
棒グラフと折れ線グラフにも登場したパラメータなので改めて説明はしない。
さて、この二つのグラフを横に並べて表示させると違いが良く解る。
棒グラフでは、各都道府県という対象に対し、鯵の漁獲量が棒の長さで表現されている。
一方、ヒストグラムでは、全体としての漁獲量が帯の面積で表現されている。
二つのグラフの「見た目」についても観察してみる。
少し、見難いが、棒グラフの場合には、各棒の間に「隙間」が存在するが、
ヒストグラムでは、各棒の間に「隙間」が存在しない。
「それで？そんなことが、何で重要やねん！ふざけんな〜！」
と言わないように。実は、重要なこと。
棒グラフは、各棒は、各対象の状況を示しているので、
隣の棒とは、独立した関係を持っている。したがって、隙間は必要。
一方、ヒストグラムは、各帯は、連続しているため、離してはならない。
あくまで、任意の区間で区切っているだけであって、実際のデータは連続している。
それゆえに、各帯の間には隙間を作らない。
それゆえに、各帯の間には隙間を作らない。
このように、一見すると、似ているように見えても、
何を表現したいのか、どのような背景を持っているかを理解すると、
両者が全く異なっていることが解る。
前置きが長くなったが、今度は、ヒストグラムの見方を考えることにする。
要するに、このグラフによって、どような情報を読み取れるか、という問題である。
ヒストグラムから読み取れる情報は、全体的なデータの偏りである。
つまり、データがどのように、分布しているかを読み取ることが目的である。
別の考え方をすると、分散を可視化する方法の一つとも考えることができる。
別の考え方をすると、分散を可視化する方法の一つとも考えることができる。
今回の場合は、鯵の都道府県別の漁獲量を示している。
このヒストグラムを見てみると、0〜100t の幅の帯が最も広く、
その頻度を見てみると、30 となっている。
したがって、全国の都道府県の内、30の都道府県では、
鯵の漁獲量が100t以内であることが解る。
一方、600〜700tもの漁獲量がある都道府県も存在していることが解る。
これらは、データ全体としては外れ値のようにも思える。
その頻度を見てみると、30 となっている。
したがって、全国の都道府県の内、30の都道府県では、
鯵の漁獲量が100t以内であることが解る。
一方、600〜700tもの漁獲量がある都道府県も存在していることが解る。
これらは、データ全体としては外れ値のようにも思える。
あまり、有用性を理解できない人もいるかもしれないが、
実は、非常に良く利用する。棒グラフよりも。
例えば、大学の成績を付ける際には、ヒストグラムを最初に描いて、
成績がどのように分布しているのかを観察し、合格ラインを考えたりする。
また、ある集落における立地環境を調べるために、
集落周辺の標高や土地利用を集計し、集落の立地特性を調べるときにも利用する。
ちょっと、特殊な使い方ではあるが。
とにかく、利用価値の高いのがヒストグラムである。
ヒストグラムについての基本的なここまで。
ここからは、ヒストグラムの描き方について話。
ヒストグラムを描く際に最初に気になるのは、
一体、何個の区切りを設けるべきかという問題である。
この区切りのことを、「階級（break）」と呼び、
その階級の値を「階級値（break point）」と呼ぶ。
何も考えずに、コンピュータの既定値を使っている人も居るが、それはNG。
ちゃんと、分け方というのも決まっている。
とは言え、階級数を決めるための理想的な方法は存在しない。
Rでは、階級数を決定するための方法がいくつか存在し、
これらを使用する以外にも、自分で階級を設定することもできる。
状況に合わせて、説明しやすいものを選ぶ。この辺は適当に。
何も説明できないのはダメであるが、説明できるのであればどの方法でも良い。
ただし、データの性質と分類方法の相性は存在する。
さて、階級数を決定する方法の中で、最もシンプルな方法は「平方根選択」。
以下の方法で選択する。説明は不要であろう。
Rでの実行方法は以下の通り。
(k <- round(sqrt(nrow(X))))
今回のデータの場合には、「6」が適切な結果となるようだ。
> (k <- round(sqrt(nrow(X))))
[1] 6
階級数は、整数で求める必要があるので、round()関数でシンプルに四捨五入。
丸め込みには他の方法も。それは、次の説明で。
ところで、この方法、数が増えるほど、階級数が大きくなってしまう。
例えば、10000の対象があると、「階級数は100」にもなってしまう。
データの全体を俯瞰するには、少々、大き過ぎはしなか？
平方根選択は、非常にシンプルな方法であるが、
最も一般的なものに、「スタージェスの公式」と呼ばれるものがある。
hist()関数の既定値。何も指定しない場合は、この公式が適用される。
なお、一般的ということで、この方法が正しいという意味ではない！
とりあえず、スタージェスの公式。
括弧の形に注目！大カッコではない。何だこれは？
大カッコに見えるが、括弧の上だけ内側に向いている。これ、「天井関数」の記号。
小難しい定義では、ある実数以上の最小の整数と定義される。
要するに、四捨五入は、0.5以上なら繰り上がり、0.4以下なら繰り下がるのだが、
天井関数では、小数点以下が0以上なら、とにかく繰り上がる。
例えば、「2.1」という数字があれば「3」となる。全部、繰り上がる。
ちなみに、「床関数」というのもあって、これは天井関数の逆。全て繰り下がる。
なお、床関数の記号は、括弧の下だけ内側に向く。
括弧の中は、二進対数となっている。Nは対象の数（総数）。
対数の話をし始めると終わらないのでいずれ。
とにかく、Rでは以下のように計算する。
# ちゃんと公式通りに計算する場合
(k <- ceiling(log2(nrow(X))+1))
ここで、ceiling()関数は、天井関数のこと。訳そのまま。
log2()関数は、二進対数を計算するための関数。
なお、Rには、nclass.Sturges()関数というものがあり、通常はそれを用いる。
# 楽をしてRのStruges()関数を使う場合
(k <- nclass.Sturges(X$Aji))
もちろん、結果は同じで「7」となる。以下は実行結果。
> (k <- ceiling(log2(nrow(X))+1))
[1] 7
> (k <- nclass.Sturges(X$Aji))
[1] 7
スタージェスの公式は、非常に良く利用されるのであるが、
データ数が200を超えるような大きなデータには不向きであり、
また、データが釣り鐘状に分布（正規分布）していない場合にも向かない。
以下の公式で計算される。これも、それほど難しくない。
右辺の方は、大丈夫であろうか？
分子は、ある変数の最大値から最小値を引いたもの。
分母の「σ（シグマ）」は標準偏差で、Nは対象の数（総数）。
これが解っていれば、この計算は簡単。Rでは以下のように計算する。
# ちゃんと公式通りに計算する場合。
(k <- ceiling((max(X$Aji)-min(X$Aji))/(3.5*sd(X$Aji)*(nrow(X)^(-1/3)))))
nclass.scott() 関数というのもある。通常はこれを使う。
# 楽をしてR のscott()関数を使う場合
(k <- nclass.scott(X$Aji))
以下が実行結果。当然、結果は同じ。「5」階級を推奨する。
> # ちゃんと計算する場合。
> (k <- ceiling((max(X$Aji)-min(X$Aji))/(3.5*sd(X$Aji)*(nrow(X)^(-1/3)))))
[1] 5
> # R のscott()関数を使う場合
> (k <- nclass.scott(X$Aji))
[1] 5
スコットの公式は、スタージェスの公式と同様に、
理論的には釣り鐘状の分布（正規分布）を前提としているが、
実用的には、スタージェスの公式と比較してその制約の度合いが小さい。
階級数を決定する方法は、他にもある。
今度は、「フリードマン＝ダイアコニスの公式（FD法）」。
ふむふむ。
これは、何のことであろうか？実は、中央値の話を見れば理解できるハズ。そうそう、分位数の話。
中央値の場合は、全データの半分番目であり、二分位に相当するが、
ここでは、Qの添字の分母が「4」となっているので、
全データを四等分した値、四分位。
つまり、3/4番目の値から、1/4番目の値を引いている。
ちにみに、この範囲のことをIQR（Interquantile Range）と呼ぶ。
日本語では「四分位範囲」と呼ぶ。頻繁に出てくる単語。
この範囲には、中央値を中心にデータの50%が収まる。
四分位は、データを観察する上で、最も基本的な指標であり、頻繁に登場する。
この話は、次の回で詳しく述べる。
さて、これも、Rで実際に計算してみる。
# ちゃんと計算する場合。
(k=ceiling((max(X$Aji)-min(X$Aji))/(2*(quantile(X$Aji)[[4]]-quantile(X$Aji)[[2]])*(nrow(X)^(-1/3)))))
四分位範囲は、quantile()関数を使って求めている。
この関数は少々特殊な形式で値を扱っているので大カッコ二つで呼び出している。
それ以外は、特に説明することもあるまい。この方法に関してもRの関数がある。
> # R のFD()関数を使う場合
(nclass.FD(X$Aji))
もちろん結果は同じ。最初の3つの方法比べて数が大きい。「15」を推奨。
> (k=ceiling((max(X$Aji)-min(X$Aji))/(2*(quantile(X$Aji)[[4]]-quantile(X$Aji)[[2]])*(nrow(X)^(-1/3)))))
[1] 15
> (nclass.FD(X$Aji))
[1] 15
フリードマン＝ダイアコニスの公式は、確率として求められる分布に対して、
誤差が最小となるように考えられた方法であり、前提を必要としない。
他の階級決定法と比較して階級数が多く設定される。
なるほど、色々とあるようだ。いやいや、実はもっとあるのだが、
ここでは、以上の方法のみを挙げておく。
さて、では、これらの結果を見比べてみよう。
# 並べて描くためのオマジナイ
par(mfrow=c(2,2))
# 平方根選択
title <- paste("Histogram of fish catches of Aji", "Square-root choice", sep="\n")
hist(X$Aji, main=title, xlab="number of fish catches", breaks=k <- round(sqrt(nrow(X))))
# スタージェスの公式
title <- paste("Histogram of fish catches of Aji", "Sturges's formula", sep="\n")
hist(X$Aji, main=title, xlab="number of fish catches", breaks=nclass.Sturges(X$Aji))
# スコットの公式
title <- paste("Histogram of fish catches of Aji", "Scott's normal reference rule", sep="\n")
hist(X$Aji, main=title, xlab="number of fish catches", breaks=nclass.scott(X$Aji))
# フリードマン＝ダイアコニスの公式
title <- paste("Histogram of fish catches of Aji", "Freedman-Diaconis' choice", sep="\n")
hist(X$Aji, main=title, xlab="number of fish catches", breaks=nclass.FD(X$Aji))
結局のところ、どれを選択するかは難しい問題であるが、
全体的なバランスで考えると、スコットの公式が無難のようである。
FD法に関しては、帯が多くなりすぎるという点に難がある。
次に、手動で設定する方法も考え得る。
それが理にかなっているのであれば構わない。
それが理にかなっているのであれば構わない。
今回の場合、実を言うと、対象数が38であり、分析するにはデータサイズが小さい。
また、データの偏りも大きいが、多くは100t以内に収まっている。
100t以上は、外れ値のようにも考えられる。
そこで、100tまでを20刻みで分け、それ以上を200刻みで分ける。
Rのhist()関数では、そのような分け方も可能である。
エラーが発生するが、以下のコマンドを実行する。
エラーの原因は後で解説する。
# 階級を設定する
b <- c(0, 20, 40, 60, 80, 100, 200, 400, 600, 700)
# タイトルをつけて、ヒストグラムを描く
title <- paste("Histogram of fish catches of Aji", "Manual breaks", sep="\n")
hist(X$Aji, main=title, xlab="number of fish catches", breaks=b, freq=TRUE)
棒グラフの場合、棒の横幅は一定であったが、
ヒストグラムでは、このように階級幅によって帯の横幅が変わる。
さて、では、エラーの原因について。
最初にも述べたが、ヒストグラムは面積の大きさが重要であった。
このヒストグラムは、面積比に問題があって、エラーの原因となっている。
最初にも述べたが、ヒストグラムは面積の大きさが重要であった。
このヒストグラムは、面積比に問題があって、エラーの原因となっている。
つまりは、どういうことであるか？
今回のhist()関数では、パラメータに freq = TRUE を指定している。
これは、ヒストグラムを頻度で描くという意味である。
要するに、ある階級幅における値の発生回数。
今度は、パラメータを明示しない場合。
どのような結果となるかを実験してみる。
今度は、エラーは発生しない。
どのような結果となるかを実験してみる。
今度は、エラーは発生しない。
# 階級を設定する
b <- c(0, 20, 40, 60, 80, 100, 200, 400, 600, 700)
# タイトルをつけて、ヒストグラムを描く
title <- paste("Histogram of fish catches of Aji", "Manual breaks", sep="\n")
hist(X$Aji, main=title, xlab="number of fish catches", breaks=b)
何やら、先程と様子が異なる。ここで、y軸ラベルとその目盛りに注目。
軸ラベルが「density（密度）」となっていて、値が非常に小さくなっている。
棒グラフというのは、棒の長さで値の大きさを表していたが、
ヒストグラムでは、棒の長さよりも、面積によって値の大きさを表現する。
したがたって、階級幅を個別に変えて、頻度を選択すると、
ヒストグラムの帯の面積比と対応する値の間に矛盾が生じてしまう。
これがエラーの原因であった。
なるほど。確かに納得できる。頻度で描いたヒストグラムでは、
直感的に、100t以上の帯の高さがどうも不自然であった。
このように、ヒストグラムが面積によって表されていることを意識するのは重要。
では、この値は一体どのような値となっているのか？
何らかの方法で基準化されているのは確かのようである。
とりあえず、次の式を見てみることにする。
この式において
は、階級値であり、今回のデータでは、0, 20, 40, 60, 80, 100, 200, 400, 600, 700 の9つの値から成るベクトルが該当する。
添字の意味は...大丈夫であろう。「i+1番目」と「i番目」のこと。
は、その階級に該当する密度である。この式を見て、意味が解ったらセンスが良い。実は、この式は、面積の和を表していて、その和が「1」となるようになっている。
うん？疑う人も居るかもしれない。
数式に慣れていない人には辛いかもしれない。
とりあえず、四角形の面積の計算はどようなものであったか？
これは大丈夫であろう。
面積＝縦×横
ふむふむ。これは納得できる。では、このヒストグラムにおける帯の面積はどうなるか？
とりあえず、一つ目の帯だけで考えることにする。
横・・・x軸なので、階級値の二番目の「20」から一番目の「0」を引いて「20」
縦・・・y軸なので、おおよそ「0.22」くらい。今は、グラフを見て目測。
したがって、一つ目の帯の面積は、おおよそ「20×0.22=4.4」くらいだろうか？
さて、これは少し大雑把だったかもしれない。ちゃんとした計算をしてみる。
次のコマンドでは、グラフを描画しない。慌ててはいけない。
# hist()関数の結果を、変数hに格納する。図は出さない。
h <- hist(X$Aji, main=title, xlab="number of fish catches", breaks=breaks, plot=FALSE)
# 変数hに格納されている中身を確認する。
str(h)
今回は、ヒストグラムの結果を変数「h」に格納にしている。
plot = FALSE は、グラフを描画しないためのパラメータ。
二つ目のstr()関数は、変数の中身を確認するための関数。
実行結果は次のようになる。
> # 変数hに格納されている中身を確認する。
> str(h)
List of 7
$ breaks : num [1:10] 0 20 40 60 80 100 200 400 600 700
$ counts : int [1:9] 17 6 4 1 2 4 2 1 1
$ intensities: num [1:9] 0.02237 0.00789 0.00526 0.00132 0.00263 ...
$ density : num [1:9] 0.02237 0.00789 0.00526 0.00132 0.00263 ...
$ mids : num [1:9] 10 30 50 70 90 150 300 500 650
$ xname : chr "X$Aji"
$ equidist : logi FALSE
- attr(*, "class")= chr "histogram"
変数「h」の中身を呼び出す方法は次の通り。
# 階級の値をベクトルとして取り出す
h$breaks
# 密度の値をベクトルで取り出す
h$density
実行結果は次の通り。
> h$breaks
[1] 0 20 40 60 80 100 200 400 600 700
> h$density
[1] 0.0223684211 0.0078947368 0.0052631579 0.0013157895 0.0026315789
[6] 0.0010526316 0.0002631579 0.0001315789 0.0002631579
ふむふむ。先程は、目測であったが、
最初の帯の密度の値は「0.0223684211」である。
なるほど、我ながら、良い目測であった。
では、面積の合計がちゃんと「1」になるかを確認する。
# まずは、結果を初期化する。
s2 = 0
# For文を回して、各帯の面積を足していく。
for(i in seq(1:(length(h$breaks)-1))){
# i番目までの総面積＝i-1番目までの総面積+（i番目の帯の高さ×i番目の帯の横幅）
s2 = s2 + h$density[i] * (h$breaks[i+1]-h$breaks[i])
# i番目までの総面積
print(s2)
}
実行結果は以下の通り。面積の計算をFor文で繰り返してて、結果を足しているだけ。
print()関数は、毎回の計算結果を逐次表示させるためのもの。無くても良い。
結果は、次のようになる。
> # まずは、結果を初期化する。
> s2 = 0
>
> # For文を回して、各帯の面積を足していく。
> for(i in seq(1:(length(h$breaks)-1))){
+ # i番目までの総面積＝i-1番目までの総面積+（i番目の帯の高さ×i番目の帯の横幅）
+ s2 = s2 + h$density[i] * (h$breaks[i+1]-h$breaks[i])
+
+ # i番目までの総面積
+ print(s2)
+ }
[1] 0.4473684
[1] 0.6052632
[1] 0.7105263
[1] 0.7368421
[1] 0.7894737
[1] 0.8947368
[1] 0.9473684
[1] 0.9736842
[1] 1
少々、回りくどい処理だったかもしれないが、
最後番目（正確には9番目）の総面積の結果が、ちゃんと「1」になっている。
なるほど、確かに面積となっている。
さて、今回も、すっかり長い話になってしまったが、
たかが、ヒストグラムに詳しい説明をしてくれる教科書は意外に少ない。
特に、階級数の決定方法は、スタージェスの公式で済ませるものも多い。
また、ヒストグラムは、多くの場合は、頻度で描いているが、密度で描く場合もある。
面積という観点から考えると、実は、密度で描いた方が良いのである。
基準化されていて、ヒスグラムとして直感的に理解しやすい。
最後に、スコットの公式を当てはめ、密度で描いた図を示す。
# 最終的に描かれたヒストグラム
# 並べて描くためのオマジナイ。全部で12個
par(mfrow=c(4,3))
for(i in seq(1:ncol(X))){
# タイトルラベルの調整
title <- paste("Histogram of fish catches of ", colnames(X)[i], sep="")
title <- paste(title, "Scott's normal reference rule", sep="\n")
# ヒストグラムを描く。色は、薄い青色で
hist(X[,i], main=title, xlab="number of fish catches",
breaks=nclass.scott(X[,i]), freq=FALSE, col="lightblue")
}
ここまでヒストグラムを詳しく説明していただきありがとうございました。大変勉強になりました。
返信削除