文系のための「数の可視化」（７）

ヒストグラムと箱ヒゲ図は、両方とも分散あるいは標準偏差の状況を
視覚的に解りやすく表現したものである。

ヒストグラムは、全体的な分布の形を視覚化し、
データの偏り方を視覚化する。

一方、箱ヒゲ図は、全体的な分布の範囲を四分位で表し、
データの偏り方を視覚化する。

分散や標準偏差というのは、全体のバラツキの程度を指標化しているだけ。
それに対して、箱ヒゲ図やヒストグラムは、全体のデータの偏り方を観察できる。

また、標準偏差からの乖離の程度を見たり、
あるいは、対応する分布を考えるためにも用いる。

これらは、統計においては、非常に重要なことである。

さて、今回は、この辺りの話には深入りせず、
もう少し、データのバラツキの可視化方法について考えてみる。

「箱ヒゲ図」と「ヒストグラム」の両方の特性を併せ持った可視化方法がある。

ここで紹介する方法は、かなり魅力的。利用価値も高い。
それにも関わらず、あまり、出回っていないから不思議である。
未だに、箱ヒゲ図を使った図しか見ない。

悪しきソフトウェア依存の習慣が原因だろうが。
まぁ、そういったことはどうでも良いことである。

さて、今回は以下の３つの方法を紹介する。

• ヴァイオリン・プロット（violin plot）
• ビーンズ・プロット（bean plot）
• ビースウォーム（beeswarm）

基本的には、箱ヒゲ図であるが、ヒストグラムの特徴も持っている。

慣れると、箱ヒゲ図よりも得られる情報量が多く、可読性も高い。

これらは、Rの初期状態では使えない。パッケージのインストールから始める。
Rコンソール上で、次のようなコマンドを入力する。

ダウンロード先の指定の画面が出てきたら、最も近い場所を選ぶ。
別に、どこを選んでも良いのだが、遠いとダウンロードに時間がかかる。

そうそう。ダウンロードだから、インターネットの接続環境が必要。

install.packages("vioplot", dep=TRUE)        # ヴァイオリンプロット
install.packages("beanplot", dep=TRUE)     # ビーンプロット
install.packages("beeswarm", dep=TRUE)   # ビースウォーム

インストール作業は、最初の一回のみ。install.packages()関数を用いる。
引数にdep=TRUEとあるが、これは依存関係の確認の有無。

Rでは、様々なパッケージを組み合わせるので、
目的のパッケージだけでは動かない場合がある。

TRUEにしておくと、足りないパッケージは自動的にインストール。

さて、インストールができたら、今度は、ライブラリを読み込む。
ライブラリの読み込みは、毎回必要。
読み込まないと使えない。ヘルプも見れない。

library(vioplot)
library(beanplot)
library(beeswarm)

では、さっそく。と言いたいところであるが...。

新しく、パッケージをインストールしたら、最初にヘルプを確認すること。

R使いは、ヘルプを見るのが当たり前。いや、見ないと何もできない。
教科書やホームページに頼ってばかりではいけない。

最初は、パッケージついての情報を確認する。

help(package="vioplot")

ここで重要なのは、「Index:」 と書かれている部分から下の部分。
左端に「vioplot」とあり、空白の右側に「violin plot」と書いてある。

実は、左端が関数名で、右側はその説明。
ふむ。vioplotの説明は、非常に不親切なようである。説明が少ない。

vioplotパッケージには、一つの関数しか含まれていないようであるが、
パッケージによっては、数十もの関数を含む場合もある。

これで、関数名が解ったので、今度はvioplot()関数のヘルプを確認する。

?vioplot        # あるいは、help(vioplot)

最も重要な部分は、「Usage:」の部分。パラメータの確認はここで。
その次に、「Arguments:」の部分で、各パラメータの説明を確認する。
最後に、「Example:」を確認して、使い方を確認する。

まぁ、理解しにくいものや、無理に理解する必要の無いものもあるのだけれど、
その辺りの判断には、ある程度の慣れが必要。

同様に、他のパッケージについても確認。これは各自で。
ひと通り、確認したところで、いよいよ、実際に使ってみる。

今回のデータは、何するか？迷うが、発電所のデータにしてみるか。

このデータは、平成22年度の「地域別発電電力量」のデータ。
総務省統計局からダウンロードしたものを編集したもの。詳細は、棒グラフの話。

これをいつものように、

# Windows と Linux の人は次のコマンド
X <- read.table("clipboard", sep=",") 

# Mac の人は次のコマンド
X <- read.table(pipe("pbpaste") , sep=",")

入力して、その後に以下をコピーして、コマンドを実行する。

pnum,max,hydro_pnum,hydro_max,thermal_pnum,thermal_max,atomic_pnum,atomic_max
Hokkaidou,66,7419,53,1234,11,4065,1,2070
Touhoku,228,17206,209,2423,13,11286,2,3274
Tokyo,192,64988,162,8981,25,38696,3,17308
Chubu,197,32828,183,5219,11,23969,1,3617
Hokuriku,138,8057,127,1904,6,4400,1,1746
Kansai,165,34877,149,8196,12,16907,3,9768
Chugoku,110,11986,97,2906,12,7801,1,1280
Shikoku,65,6963,58,1141,4,3797,1,2022
Kyushu,194,20330,139,3279,45,11577,2,5258
Okinawa,22,1919,0,0,21,1919,0,0

前処理として、色を付けるための色の定義を。
これは、棒グラフの話で出てきた色の定義方法。

col.h <- rgb(50, 100, 255, maxColorValue = 255)  # 水力発電所のイメージ色
col.t <- rgb(255, 100, 50, maxColorValue = 255)  # 火力発電所のイメージ色
col.a <- rgb(200, 255, 100, maxColorValue = 255)  # 原子力発電所のイメージ色??

まずは、通常の箱ヒゲ図から。今回は発電所種別ごとの総出力で見てみる。

boxplot(X[,c(4,6,8)], names=c("Hydro","Thermal","Atomic"), col=c(col.h, col.t, col.a))

title(main="Maximum output of electronic plants in Japan", ylab="Maximum output (Kw)")

これは大丈夫のはず。これで、何の変哲もない、普通の箱ヒゲ図が描けた。

次に、ヴァイオリン・プロット。描き方はかなり異なる。

vec.h <- as.vector(X[,4])    # 元のデータから水力発電の総出力をベクトルで抽出
vec.t <- as.vector(X[,6])    # 元のデータから火力発電の総出力をベクトルで抽出
vec.a <- as.vector(X[,8])    # 元のデータから原子力発電の総出力をベクトルで抽出

vioplot(vec.h, vec.t, vec.a, names=c("Hydro","Thermal","Atomic"), col="lightpink")
title(main="Maximum output of electronic plants in Japan", ylab="Maximum output (Kw)")

violot()関数は、ベクトルデータを並べて指定する。
そのため、vec.h、vec.t、vec.a の３つのベクトルを作成し、
それぞれ、水力発電、火力発電、原子力発電の総出力を代入。

実は、色を個別に付けるのは無理なので、ここでは、薄いピンク色。
ソースコードが公開されているので、修正することは難しくないが、
後々になって、説明することが面倒なので、そのままで。

中心の白い「○」が中央値。
真ん中の黒い「帯」がIQR（四分位範囲）。
真ん中を突っ切る「線」がIQRの±1.5倍のヒゲ。

左右に広がる薄いピンクの「山」は、カーネルと呼ばれる密度推定の山。
これが、ヒストグラムに相当する。カーネルの話はいずれ。

vioplot()関数は、ヒストグラムの特性と箱ヒゲ図の特性を合わせた特性を持ち、
重要なこととして、箱ヒゲ図で表現されている全て情報も表現している。
問題は、データの指定方法で、この関数だけデータの指定方法が特殊である。

今度は、ビーンプロットを試してみる。

beanplot(X[,c(4,6,8)], names=c("Hydro","Thermal","Atomic"), col="lightblue", 

    main="Maximum output of electronic plants in Japan", ylab="Maximum output (Kw)", cut=0)

beanplot()関数でのデータの選択は、箱ヒゲ図と同じ。

このグラフにおいて、真ん中の太い線は、中央値を表し、

複数の細い「線」は、データが存在している場所を示している。

横に広がる「山」は、vioplot()関数と同様にカーネルの山。

vioplot()関数と比較すると beanplot()関数の山は細かく波を打っている。

これは、カーネルのパラメータの違い。カーネルは、様々な描き方ができる。

beaswarm()関数は、他の二つの方法と比較して、

分布に重点を置いているため、より詳細な密度推定が可能となる。

その一方で、箱ヒゲ図と比較して、外れ値やIQRの表現力は低い。

注意するべきは、cut パラメータ。この値は「0」にする。

既定のままだと、データが存在しない所まで山裾が広がってしまう。

色に関しては、複数のbean（マメ）の色を変えることはできず、

title()関数を使って、後から、タイトルを入れることもできない。

最後はビースウォーム。

beeswarm(X[,c(4,6,8)], labels=c("Hydro","Thermal","Atomic"), col=c(col.h, col.t, col.a), pch=16)

title(main="Maximum output of electronic plants in Japan", ylab="Maximum output (Kw)")

beeswarm()関数は、分布をカーネルの山で表現せず、実際の点を布置する。

今回の場合は、対象数が少ないので、見難いが、

対象数がある程度あれば、この方法が最も誤解無く情報を伝えることができる。

boxplot()関数と異なるのは、ラベルの指定パラメータ。

boxplot()関数がnames であるのに対し、beeswarm()関数では labels となっている。

beeswarm()関数は、他の二つの方法と比較して、推定を用いていないため、

より実態に即した情報を伝えることができるが、外れ値やIQRの表現はできない。

ただし、色の指定は可能なので、色分けをしたい場合に使えるというメリットがある。

さてさて。では、今回紹介した方法のうちで、最も適切なものはどれか？

結論から言うと、ケース・バイ・ケースであるが、無難な方法はビースウォーム。

うん？色々と欠点があるのに？

実は、ビースウォーム・プロットの欠点は、

箱ヒゲ図と重ね合わせることで克服できる。

boxplot(X[,c(4,6,8)], names=c("Hydro","Thermal","Atomic"))

beeswarm(X[,c(4,6,8)], labels=c("Hydro","Thermal","Atomic"), col="red", pch=16, add=TRUE)

title(main="Maximum output of electronic plants in Japan", ylab="Maximum output (Kw)")

beeswarm()関数のパラメータの最後に add=TRUE としてある。

このようにすることで、前に描画したプロットに、新しいプロットを追加できる。

なお、この方法は、他の方法にも適用できるが、

ヴァイオリン・プロットやビーンプロットを箱ヒゲ図に重ねるメリットは無い。

文系のための「数の可視化」（６）

さて、今回は「箱ひげ図」について話をする。
他のグラフと比較すると馴染みが無いかもしれない。

箱ヒゲ図は、四分位あるいはヒンジに対応する可視化の方法であり、
ヒストグラムと同様に、バラツキに対応する可視化の方法である。

分散や標準偏差との対応という点では、
ヒストグラムよりも、より、バラツキの状況を良く表している。

では、早速、データの準備を準備する。
今回も、ヒストグラムの話で用いたデータと同じデータを用いる。
このデータは、農林水産省の「都道府県別魚種別遊漁採捕量」のデータを、

本ブログでの説明用に加工したものであり、他の用途には適していない。

いつものように、

# 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

では、箱ひげ図を書いてみよう。まずは、何の設定もしない。

boxplot(X)

とにかく、何かが表示されたが、色々と欠けている。
何が足りないか、あるいは、不適切であるか、解っているハズ。

boxplot(X, main="Distribution of fish cathes", ylab="Number of fish cathces (tones)", las=2)

ふむ。最初のグラフには、タイトル、y軸ラベルとy軸の単位が記載されていない。
そのようなグラフを描いてはいけない。

欠けているわけではないが、項目名が全て表示されていないのも問題である。
そこで、las=2 というパラメータを設定し、
項目ラベルの方向をx軸に対して垂直になるように再配置している。

ところで、この状態でも悪くは無いのだが、項目ラベルが横向きであると、
少々、見難い。首を90度に曲げてみることになる。実は、横方向に描くことができる。
そのようにするには、horizontal パラメータを TRUE に設定すれば良い。

# 項目ラベルが画面からはみ出るので余白を調整
par(mar=c(5,5,3,1))

# 箱ひげ図を描画する。
boxplot(X, main="Distribution of fish cathes", xlab="Number of fish cathces (tones)", las=1, horizontal=TRUE)

今度は、xとyの軸が入れ替わるので、ylabとlasのパラメータを調整している。
先程の箱ひげ図と比較すると、かなり見易い。

実を言うと、描画に関してはここまで。後は、色を付けるくらいであろうか？

さて、冒頭にも述べたように、箱ひげ図は、四分位あるいはヒンジに対応する。
どちらを採用するかは、ソフトウェアに依存するが、Rではヒンジを使っている。

今度は、説明のために、鯛の箱ひげ図だけを描画してみる。

boxplot(X$Tai, main="Distribution of fish cathes (The case of Aji)", ylab="Number of fish cathces(tones)")

今回は、変数が一つしか存在しないので、対象ラベルを表示できない。

仕方が無いので、グラフタイトルに鯛の箱ひげ図であることを示す。

では、この図を見ながら、箱ひげ図を見ていくことにする。

中央に箱があって、その真ん中に太い横線が入っている。
箱からは、上下に点線（ヒゲ）が伸びていて、点線の端には横棒が描いている。
あと、上の方には丸い点がある。

では、真ん中の箱から。これは四分位範囲（IQR）に対応し、
箱の下の辺は第1四分位数を表し、箱の上の辺は第3四分位数を表している。
なお、Rの場合には、下側ヒンジと上側ヒンジの範囲である。

真ん中の太い横棒は、中央値を表している。

下向きのヒゲの端にある横線は、外れ値を除いた最小値、
上向きのヒゲの端にある横線は、外れ値を覗いた最大値を示している。

そして、上の方についてある丸が、上方向の外れ値となる。
このデータには、存在しないが、下方向の外れ値も取り得る。

うん？ここで、疑問に感じた人はセンスが良い。
外れ値は、どのようにして決められているのか？以外に難しい問題である。
実際に、色々な方式があって、ソフトウェアにも依存する。

ここでは、一般的によく用いられている説明を紹介する。

その方法とは、「IQR」を基準にする方法である。
この範囲を基準に、外れ値の範囲を決める。

まずは、軽度の外れ値。通常は「○」記号で表す。

• 軽度の負の外れ値：
• 軽度の正の外れ値：

次は、軽度の外れ値。通常は「＊」記号で表す。

• 重度の負の外れ値：
• 重度の正の外れ値：

以上の外れ値の定義から、ヒゲの長さは、最大でもIQRの1.5倍よりも短い。
Rの場合、軽度と重度の外れ値を分けず、IQRからの離れ具合は、
rangeパラメータで決める。このパラメータは既定では「1.5」である。

以上が、箱ヒゲ図を読み解く上で、必要最低限の情報である。

さて、この鯛の箱ヒゲ図から、どのような情報が読み取れるかを検討してみる。
図だけだと、色々と具体的な数値が見えないので、
一度、箱ヒゲ図の結果を別の変数に格納し、結果を表示する。

# まずは、変数に結果を格納する。プロットは不要。
b.Tai <- boxplot(X$Tai, plot=FALSE)

# 格納されているデータを確認する。
str(b.Tai)

一連の処理結果は以下の通り。

> # まずは、変数に結果を格納する。プロットは不要。
> b.Tai <- boxplot(X$Tai, plot=FALSE)
>
> # 格納されているデータを確認する。
> str(b.Tai)
List of 6
 $ stats: integer [1:5, 1] 0 7 33.5 123 210
 $ n    : num 38
 $ conf : num [1:2, 1] 3.77 63.23
 $ out  : num [1:2] 308 306
 $ group: num [1:2] 1 1
 $ names: chr "1"

必要なデータは、$stat という所に入っているようなので、

b.Tai$stat

中身は以下の通り。

> b.Tai$stat
      [,1]
[1,]   0.0
[2,]   7.0
[3,]  33.5
[4,] 123.0
[5,] 210.0
attr(,"class")
        1
"integer"

ふむふむ。なるほど。最後の結果は、外れ値を除いた状態になっている。
では、箱ヒゲ図に描かれている情報を整理する。

まず、この箱ヒゲ図では、IQRが全体的に下方向に偏っていて、
また、中央値もIQRの下側に寄っていることから、

全体的に低い値を取る対象が多く、大きい値を取るものは少ないことが解る。

外れ値を除くと、データは0〜200t の間に収まり、
全データの約50%は、7〜123t の幅に収まっている。

この状況をヒストグラムと比較してみる。
両者を並べて表示させると、状況がより明確に見えてくる。

# ヒストグラムと箱ヒゲ図のレイアウト枠を設定（5:1）
lay <- matrix(c(1,1,1,1,1,2), ncol=1)
layout(lay)

# ヒストグラムをグラフ領域の上に描く。
par(mar=c(4,5,3,1), cex=1.1)
hist(X$Tai, breaks=nclass.FD, freq=FALSE, col="pink", main="", xlab="")
title(main="Density of fish catches for Tai")

# 箱ヒゲ図をヒストグラムの下に描く
par(mar=c(0,5,0,1), bty="n")
boxplot(X$Tai, main="", xlab="", ylim=c(0,350), horizontal=TRUE, xaxt="n", col="pink")
mtext(side=3, "Number of fish catches")


layout()関数は、描画領域の番号を成分とした行列を指定することで、
グラフ領域に複数のグラフを表示することができる関数である。

今回の場合は、最初に「6行✕1列」の行列を作成し、
その成分には、上から5行までは1を、最終行には2を入れてある。

したがって、グラフ領域は、縦方向に6等分され、

上の5つ分のスペースには1つ目のプロット（ヒストグラム）が、
下の1つ分のスペースには2つ目のプロッタ（箱ヒゲ図）が、それぞれ描かれる。

こうして、上の図のように、二つの図を並べて書くことができる。
以前に、par()関数のmfrowパラメータを用いたが、
layout()関数のように、描画領域を比率を変えて設定することはできない。

さて、このように二つのグラフを並べて表示することで、
データのバラツキの状況が良く解る。

実は、箱ヒゲ図とヒストグラムの両方の特性を持ったグラフも存在するが、
その話は、次の話に置いておくことにする。

文系のための「データのバラツキ」（２）

ようやく、データ分析の入り口に辿り着いた。
ここまでの話は、基礎的な話が多かった。

今日の話は、データの「要約」について。

元のデータを分析する前に、データの全容を把握することは重要である。
本来、データの可視化は、全体的な状況を把握してから行うべきであって、
可視化してから考えるのではない。

まず、基本的な情報として考えられるものは何だろうか？
これまでに、やってきた情報を整理すると、

• 行数＝対象の数
• 列数＝属性（変数）の総数
• 各変数の平均と中央値
• 各変数の分散と標準偏差
• 各変数間の共分散と相関係数

と言ったところか。

これだけの情報が得られるだけで、様々なことが解るのだが、
今回は、これらの指標に加えて、「分位数（quantile）」について考えてみる。

分位数は、データのバラツキを、分散とは違った方法で観察する方法とも言える。

とにかく、説明するよりも、手を動かした方が理解しやすい。
とりあえず、いつものように...

# 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

このデータは、ヒストグラムの話で用いたデータと同じ。
農林水産省の「都道府県別魚種別遊漁採捕量」のデータを編集して作ったデータ。
架空とまでは言わないが、現実の分析には適さない。

では、早速。分位数を出してみる。まずは、鯵の分位数を計算。

quantile(X$Aji, probs=seq(0,1,by=0.25))

ここでは、quantile()関数を用いて分位数を計算している。
quantile は、日本語では「分位数」を意味する。

この何分割するかを決めているのが、probs パラメータである。
四分位の場合には、probs=seq(0,1,by=0.25) となっているということは、
0〜1までを、0.25刻みで「四等分」している。

このように、4つに分ける分位数のことを「四分位（quartiles）」と呼ぶ。

つまり、0.00, 0.25, 0.50, 0.75, 1.00 というベクトルによって指定しているわけである。
ここで、注意するべきは、0番目の値が存在すること。したがって、
四分位の場合は、5つの値が返ってくる。

中央値の話でも述べたが、様々な分位数が存在し、それぞれに名前が付いている。

• 二分位（median）＝中央値
• 三分位（tertiles） 
• 四分位（quartiles）
• 五分位（quintiles）
• 六分位（sextiles）
• 十分位（deciles） 
• 十二分位（duo-deciles）
• 二十分位（vigintiles）
• 百分位（percentiles）
• 千分位（permiles）

Rの quantile()関数で、これらを計算することは難しくない。

quantile(X$Aji, probs=seq(0,1,by=0.5))    # 中央値
quantile(X$Aji, probs=seq(0,1,by=0.2))    # 五分位
quantile(X$Aji, probs=seq(0,1,by=0.1))    # 十分位

実行結果は以下の通り。

> quantile(X$Aji, probs=seq(0,1,by=0.5))    # 中央値
  0%  50% 100%
   0   23  683
> quantile(X$Aji, probs=seq(0,1,by=0.2))    # 五分位
   0%   20%   40%   60%   80%  100%
  0.0   0.4  19.8  40.2 102.6 683.0
> quantile(X$Aji, probs=seq(0,1,by=0.1))    # 十分位
   0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100%
  0.0   0.0   0.4  14.2  19.8  23.0  40.2  57.8 102.6 218.0 683.0

色々とあるが、四分位が最も一般的であるので、
quantile()関数の既定は四分位。したがって、パラメータ指定も不要

quantile(X$Aji)

パラメータ指定無しで、このようにすると、以下のように四分位が計算される。

> quantile(X$Aji)
  0%  25%  50%  75% 100%
   0    1   23   82  683

ここで、四分位の考え方を整理してみる。

• 全体の下から0/4（0%）番目の数：第0四分位数（最小値）
• 全体の下から1/4（25%）番目の数：第1四分位数
• 全体の下から2/4（50%）番目の数：第2四分位数（中央値）
• 全体の下から3/4（75%）番目の数：第3四分位数
• 全体の下から4/4（100%）番目の数：第4四分位数（最大値）

となる。ここで、第1四分位数〜第3四分位数の範囲のことを四分位範囲と呼ぶ。

これは、データの下から25%番目目と75%番目の間の範囲である。
別の言い方をするならば、中央値を中心に全データの50%が収まる範囲である。

英語では、Inter-Quantile Range と呼び、IQRと表記する。
Rには、IQR()関数が存在し、四分位範囲も容易に計算できる。

IQR(X$Aji)

実行結果は、以下の通り。

> IQR(X$Aji)
[1] 81

今回のデータの場合は、1〜82までの範囲なので、
四分位範囲は、「82-1=81」となり、正常に計算されている。

やはり、数理的な話も必要であろう。ここからは、数学的な話。

実は、色々と前提となる知識を必要とするので、
色々と正しくない表記も含まれる。表記方法は変則的。厳密な説明ではない。
Wikipedia の解説では、分位数を関数として定義し、実数空間への投影で説明。

意味の解る人は、Wikipedia へ。
意味の解らない人は、以下の説明へ。

まず、以下のようなベクトルがあったとして、



これを昇順に並べ変えたものを次のように定義する。



総数を「n」とし、求めたい分位数を「q」とすると、
次のように表すことができる。



うん？添字の部分が...何番目？落ち着いて考えてみる。

通常は、1から数えるが、コンピュータでは0番目を基準とする場合が多い。
そこで、「(n-1)q」のところで、一番目を「0」からの順番に書き直している。

少し、混乱するかも。要するに、元のデータが以下のようになっていて、



これを、「0」からの順番で書きなおし、



「1」ずらした分を後で、正しい位置に戻すために「+1」としている。
少々、回りくどい表現であるが仕方あるまい。

さて、これで計算できるかと言うと、実は以下のような条件が付随する。
添字部分が、「(n-1)q+1」だと解りにくいので、とりあえず「t」と置く。



見慣れない記号が一杯出てくる。とりあえず、順番に。

まず、「」の記号は、ある集合に含まれていることを表す記号で、
逆が、「」 の記号。つまり、ある集合に含まれていない状況を表す。

では、一体何の集合なのかというと、 という記号が教えてくれている。
これは、Natural の「N」だから、自然数を表している。自然数は1以上の整数。

要するに、上の条件というのは、「整数」だったらそのままで良いけれど、
「小数点」の値が出てきた場合には、再配分し直しましょう、と言っている。

では、早速、この式に当てはめてRで確認する。
仮に、第1四分位を求めたいのであれば、q = 1/4 = 0.25 となる。

# まずは、データをソートする。
x <- sort(X$Aji)

# 第1四分位数を算出する
(t <- (length(x)-1)*0.25+1)

実行結果は以下の通り。

> (t <- (length(x)-1)*0.25+1)
[1] 10.25

なるほど。見事に、自然数では無い。
ということで、厄介な自然数で無い場合の処理。

((ceiling(t)-t)*x[floor(t)]+(t-floor(t))*x[ceiling(t)])

そして、その結果は、以下の通り。

> ((ceiling(t)-t)*x[floor(t)]+(t-floor(t))*x[ceiling(t)])
[1] 1

確かに、合っているようだが、「1」というのは特殊すぎて自身が無い。
そこで、第1四分位〜第4四分位を連続して計算してみる。

x <- sort(X$Aji)    # まずは、データをソートし、ベクトルxに格納
n <- length(x)        # ベクトルxの数を数える

# For文による繰り返し処理
for(q in seq(0, 1, by=0.25)){
    # 四分位数を算出する
    t <- (n-1)*q+1

    # 自然数かをチェックする。
    if (t > 0 && identical(round(t),t)==TRUE) {
        # 自然数の場合はそのまま
        print(x[t])
    } else {
        # 自然数で無い場合は調整する
        print((ceiling(t)-t)*x[floor(t)]+(t-floor(t))*x[ceiling(t)])
    }
}

さて、ここでは、length()関数とif文、identical()関数が初登場である。

まず、length()関数は、ベクトルに含まれるスカラーの個数を返す関数である。
行列のnrow()関数やnrow()関数に似ているが、行列に対してlength()関数を使うと、
その行列に含まれる全成分の総数を返す。

 identical()関数は、二つの値が同じあるかを確認するための関数であり、
同じ場合には、「TRUE」を返し、異なる場合には「FALSE」を返す。
ここでは、四捨五入した値と元の数字を比較するために用いている。

if文は、for文と同様に、プログラミングの基本となる構文であり、
ある処理の条件分岐を定義する。

ここでは、「t」が自然数であるか、否か、の分岐を示している。
括弧内の処理が、TRUEの場合は、そのままの値を出力し、
それ以外の場合（else）の時には、値を調整して出力している。

この分岐処理の内部では、二つの条件が含まれている。
一つ目は、t > 0 という条件（1以上の正の整数であるという条件）
二つ目は、t が小数点を含んでいるか否か？という条件である。

この二つの条件は、「&&」で結ばれている。これをAND演算と呼ぶ。
すなわち、両方の条件が揃った場合を「真」とする条件が与えられている。

さて、この処理の実行結果は以下の通り。

[1] 0
[1] 1
[1] 23
[1] 82
[1] 683

なるほど、Rの quantile()関数 の結果と見事に同じである。
以上が、分位数の算出方法である。

ところで、今回は、分位数の算出方法を、



で説明した。Excel や R の既定（type 7）はこの方法に従っている。

しかしながら、分位数の分割方法は9通りも存在し、
この辺りは実に厄介な話である。実を言うと、良く理解できない。

この四分位に対して、ヒンジという考え方がある。
これは、四分位と異なり、考え方は至ってシンプル。

中央値より上側の中央値を上側ヒンジと呼び、
中央値より下側の中央値を下側ヒンジと呼ぶ。

R では、fivenum()関数によって計算することができる。

fivenum(X$Aji)

実行結果は以下の通り。

> fivenum(X$Aji)
[1]   0   1  23  88 683

この方法は、五数要約と呼ばれ、fivenum()関数は、
最小値、下側ヒンジ、中央値、上側ヒンジ、最大値、の5つの値を返す。

実用上は、四分位と五数要約を意識することは無いが、
箱ひげ図の場合には、五数要約を用いる。

最後に、今回のデータの四分位をまとめて出力してみる。
これは、summary()関数を用いることで、簡単に計算できる。

summary(X)

行列データから、直接的計算してくれる。非常に楽。
よく見ると、四分位だけではなく、平均値も計算されている。

> summary(X)
      Aji              Saba             Buri           Hirame       
 Min.   :  0.00   Min.   :  0.00   Min.   :  0.0   Min.   :  0.000  
 1st Qu.:  1.00   1st Qu.:  0.00   1st Qu.:  3.0   1st Qu.:  0.000  
 Median : 23.00   Median :  5.00   Median : 23.0   Median :  1.000  
 Mean   : 84.08   Mean   : 26.21   Mean   : 79.5   Mean   :  9.447  
 3rd Qu.: 82.00   3rd Qu.: 27.00   3rd Qu.: 81.0   3rd Qu.:  8.500  
 Max.   :683.00   Max.   :185.00   Max.   :823.0   Max.   :100.000  

     Karei             Kasago           Mebaru          Kawahagi     
 Min.   :   0.00   Min.   :  0.00   Min.   :  0.00   Min.   :  0.00  
 1st Qu.:   0.00   1st Qu.:  0.00   1st Qu.:  0.50   1st Qu.:  0.00  
 Median :   0.00   Median :  6.50   Median : 11.00   Median :  1.00  
 Mean   :  94.74   Mean   : 22.32   Mean   : 44.74   Mean   : 11.24  
 3rd Qu.:   1.75   3rd Qu.: 31.25   3rd Qu.: 63.75   3rd Qu.:  4.00  
 Max.   :2140.00   Max.   :131.00   Max.   :452.00   Max.   :217.00  

      Tai             Isaki             Ika              Tako       
 Min.   :  0.00   Min.   :  0.00   Min.   :  0.00   Min.   : 0.000  
 1st Qu.:  8.75   1st Qu.:  0.00   1st Qu.:  0.25   1st Qu.: 0.000  
 Median : 33.50   Median :  2.00   Median :  8.50   Median : 0.000  
 Mean   : 72.26   Mean   : 25.08   Mean   : 69.39   Mean   : 5.026  
 3rd Qu.:121.25   3rd Qu.: 19.75   3rd Qu.: 60.50   3rd Qu.: 2.000  
 Max.   :308.00   Max.   :205.00   Max.   :660.00   Max.   :49.000  

四分位は、分散と同様にデータ全体のバラツキを観察する上で重要であるが、
この状況を、眺めても、いまいち解りにくい。やはり、標準偏差の方が解かりやすいか？

しかしながら、この状況を可視化することで、非常に重要な参考資料となり得る。
その可視化方法こそが、「箱ひげ図」なのである。この話については、次の話で。

文系のための「数の可視化」（５）

基本的なグラフとして理解しておくべきは、
棒グラフ、円グラフ、折れ線グラフの三種類であるが、
グラフには様々なものがあり、状況に合わせて適切に使うことが重要。

さらに、本ブログにおける一連の投稿記事では、上記の三種類に加えて、
ヒストグラム、箱ひげ図、散布図についても解説する。

今回は、ヒストグラムについて整理する。

ヒストグラムは、どうも、誤解の多いグラフである。
見た目が棒グラフに似ているので、棒グラフと混同している人も多い。
しかしながら、ヒストグラムと棒グラフは似て非なる物である。

少なくとも、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もの漁獲量がある都道府県も存在していることが解る。
これらは、データ全体としては外れ値のようにも思える。

このように、データの全体を把握するために用いるのがヒストグラムである。
あまり、有用性を理解できない人もいるかもしれないが、
実は、非常に良く利用する。棒グラフよりも。

例えば、大学の成績を付ける際には、ヒストグラムを最初に描いて、
成績がどのように分布しているのかを観察し、合格ラインを考えたりする。

また、ある集落における立地環境を調べるために、
集落周辺の標高や土地利用を集計し、集落の立地特性を調べるときにも利用する。
ちょっと、特殊な使い方ではあるが。

とにかく、利用価値の高いのがヒストグラムである。

ヒストグラムについての基本的なここまで。
ここからは、ヒストグラムの描き方について話。

ヒストグラムを描く際に最初に気になるのは、
一体、何個の区切りを設けるべきかという問題である。
この区切りのことを、「階級（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"

とりあえず、ここで必要となるのは、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")
}