Photo by Robert Thiemann on Unsplash
RのHistDataパッケージの中のChestSizesです。
Queteletという人が調べたスコットランドの軍人のチェストサイズ、胸囲ですかね？のデータセットです。
このデータをつかって自然のデータが正規分布になっていることをデモンストレーションしたそうです。
データを読み込みます。
str()関数で確認すると、16行2列のデータフレームでしたので全部表示してみました。
chestがサイズでその人数がcountですね。
では、ヘルプのコードを実行してみます。
plot()関数をつかって横軸がchest, 縦軸がcountのグラフを描きました。type = "b" としているので、線と丸印の両方が描かれています。2列のデータフレームだと1列目が横軸で2列目が縦軸として描かれるんですね。
試しに適当なデータでやってみます。
予想どおり2列のデータフレームだと、1列目が横軸、2列目が縦軸になりました。
では、3列のデータフレームではどうでしょうか？
3列だと散布図マトリックスになるんですね。
脱線してしまいました。ヘルプのコードに戻ります。
barplot()関数で棒グラフを描いています。正規分布っぽいかたちですよね。
次のヘルプのコードは、
です。cestが正規分布するとしたらどういう分布になるか、期待値を計算しています。
ここまででヘルプのコードは終わっています。
この期待値と実際の値をくらべましょう。
黒の実線が実際の値で赤の点線が期待値です。非常によく似た形状ですよね。
もう一つ、barplotに期待値を赤の点線で重ねてみましょう。
barplot()関数で生成された棒グラフは、横軸の座標がよくわからないので、いったんbplotという名前のオブジェクトとして保存しています。
このbplotというオブジェクトをstr()関数、class()関数で調べると、16行1列のマトリックスであることがわかりました。そしてこの値が横軸の座標でした。
そこで、lines()関数で横軸をbplot, 縦軸をexpectedにして棒グラフに赤い点線を載せました。
これでおしまいです。
以下に今回のコードを記載しておきます。
#
library(HistData)
data("ChestSizes")
#
str(ChestSizes)
ChestSizes
#
# frequency polygon
plot(ChestSizes, type = "b")
#
# 適当なデータで実験
plot(data.frame(number = 1:16, count = rnorm(16)), type = "b")
#
# 3列で実験
plot(data.frame(r1 = 1:16, r2 = rnorm(16), r3 = rep(1:4, 4)))
#
# barplot
barplot(ChestSizes[ , 2], ylab = "Frequency", xlab = "Cest size")
#
# calculate expected frequencies under normality,
# chest ~ N(xbar, std)
n_obs <- sum(ChestSizes$count)
n_obs
xbar <- with(ChestSizes, weighted.mean(chest, count))
xbar
std <- with(ChestSizes, sd(rep(chest, count)))
std
expected <- with(ChestSizes,
diff(pnorm(c(32, chest) + 0.5, xbar, std))*sum(count))
expected
#
# plot actual data and expected data
plot(ChestSizes, type = "b", lwd = 1)
lines(ChestSizes$chest, expected, lty = 2, col = 2)
#
# barplot objects
bplot <- barplot(ChestSizes[ , 2], xlab = "chest", ylab = "frequency")
str(bplot)
class(bplot)
bplot
lines(bplot, expected, col = "red", lwd = 3, lty = 2)
