度数分布の問題

度数分布（frequency distribution）、度数分布表（frequency table）、度数分布図（histogram）について勉強します。“数学チュートリアル やさしく語る 確率統計 西岡康夫（著）”の1.4節 度数分布の練習問題をRで解きます。

練習問題
ある職場における女性社員の身長（cm: 小数点以下は四捨五入）は以下のようであった。
165 157 149 158 159 162 161 172 163 155
161 159 153 159 163 160 168 157 164 158
このデータについてヒストグラムを作れ。

階級分け

度数分布はデータを階級（class）に分けて、それぞれの階級に属するデータの個数（度数 frequency）をカウントしたリストです。階級幅（class interval）はデータの個数\(n\)の平方根\(\sqrt{n}\)を目安にします。

変数dataにデータを代入します。

data <- c(165, 157, 149, 158, 159, 162, 161, 172, 163, 155, 161, 159, 153, 159, 163, 160, 168, 157, 164, 158)

階級幅の目安を計算します。

sqrt(length(data))
## [1] 4.472136

4.47ですので、切りが良い5 (cm)を階級幅とします。

そして、データの最小値と最大値を考えて、すべてのデータを含む6個の階級に分けることにします。

min(data)
## [1] 149
max(data)
## [1] 172
# 階級
1 145cm以上150cm未満
2 150cm以上155cm未満
3 155cm以上160cm未満
4 160cm以上165cm未満
5 165cm以上170cm未満
6 170cm以上175cm未満

階級に属するデータをそれぞれ抽出します。

sd <- sort(data)
class1 <- subset(sd, 145 <= sd & sd < 150)
class2 <- subset(sd, 150 <= sd & sd < 155)
class3 <- subset(sd, 155 <= sd & sd < 160)
class4 <- subset(sd, 160 <= sd & sd < 165)
class5 <- subset(sd, 165 <= sd & sd < 170)
class6 <- subset(sd, 170 <= sd & sd < 175)
class1
## [1] 149
class2
## [1] 153
class3
## [1] 155 157 157 158 158 159 159 159
class4
## [1] 160 161 161 162 163 163 164
class5
## [1] 165 168
class6
## [1] 172

階級別に昇順にソートしたデータをまとめます。

# 階級 データ
1 145cm以上150cm未満 149
2 150cm以上155cm未満 153
3 155cm以上160cm未満 155, 157, 157, 158, 158, 159, 159, 159
4 160cm以上165cm未満 160, 161, 161, 162, 163, 163, 164
5 165cm以上170cm未満 165, 168
6 170cm以上175cm未満 172

度数分布表の作成

階級にまとめたデータを度数分布表にします。度数分布表はデータを階級別に分類して、階級値と度数を表で表したものです。階級値（class value）は、各階級の最大値と最小値の相加平均（階級幅の中心値）です。

# 1 2 3 4 5 6
階級 [145, 150) [150, 155) [155, 160) [160, 165) [165, 170) [170, 175)
階級値 147.5 152.5 157.5 162.5 167.5 172.5
度数 1 1 8 7 2 1

度数は、length(class#)で求まります。

ヒストグラムのプロット

度数分布表をグラフにしたものが度数分布図（ヒストグラム）です。

グラフをプロットするためにパッケージggplot2をインストールします1 2。

if (!require("ggplot2")) {
  install.packages(c("ggplot2"), repos="http://cran.rstudio.com/")
  library(ggplot2)
}
## Loading required package: ggplot2

上記で求めた度数分布表をヒストグラムにします3。

# Make data frame for ggplot
cv <- c(147.5, 152.5, 157.5, 162.5, 167.5, 172.5)
freq <- c(length(class1), length(class2), length(class3), length(class4), length(class5), length(class6))
df <- data.frame(cv, freq)
# Add plot layers
g <- ggplot(df, aes(x = cv, y = freq))
g <- g + geom_bar(width = 5, stat = "identity") # set bar chart
g <- g +  scale_x_continuous(breaks = seq(147.5, 172.5, by=5), limits=c(145, 175)) # ticks
g <- g + xlab("Class value (cm)") # x label
g <- g + ylab("Frequency")        # y label 
g <- g + ggtitle("Histogram by manual")     # title
# Plot
plot(g)

度数のカウントを自動で行うこともできます。

df2 <- data.frame(data)
g <- ggplot(df2, aes(x = df2$data))
g <- g + geom_histogram(binwidth = 5) # set histogram chart
g <- g + xlab("Class value (cm)") # x label
g <- g + ylab("Frequency")        # y label 
g <- g + ggtitle("Histogram by ggplot") # title
plot(g)

このように階級の取り方で度数分布は変わります。

hist関数だけでもヒストグラムをプロットすることができます。

hist(data)

  1. 公式ドキュメント: Index. ggplot2 2.1.0↩

  2. r - Elegant way to check for missing packages and install them? - Stack Overflowにあるif (!require("package")) {}を真似しました。このif文がないとR Markdownをhtmlに変換する度にパッケージがダウンロードされます。↩

  3. ggplot2で棒グラフや積み上げ棒グラフを描く方法 - biostatisticsを参考にしました。↩