www.crosshyou.info

政府統計の総合窓口のデータや、OECDやUCIやのデータを使って、Rの練習をしています。ときどき、読書記録も載せています。

HistDataパッケージのJevons

Photo by Setu Chhaya on Unsplash 

HistDataパッケージのJevonsというデータは、1871年のNatureという雑誌に掲載されたW. Stanley Jevonsの実験のデータです。

黒いビーズを複数個、パッと見せて、何個だったか答えさせるという実験です。人間が一度に認識できる数がどのくらいか、というのを調べる実験です。一般に7個前後、と言われていますね。

早速データを見てみます。

actual が実際のビーズの数、estimatedが答えた数、frequencyは度数、errorは誤差ですね。

一番始めの行は、actual:3, estimated:3, frquency:23, error:0 なので、実際のビーズが3のとき、3と答えたのが23回あって、誤差は0だった、という意味です。

ヘルプのコードを実行していきましょう。

はじめは、xtabs()関数で表を作っています。

ビーズが、3つ、4つのときは全部正解していますが、5つのときは102回正解で、5回間違っていますね。15個のときは正解は2回で9回間違っています。

次のヘルプのコードは、サンフラワープロット、sunflowerplot()関数です。

lm()関数で回帰直線モデルを作って、abline()関数でその直線を描いています。

次は、gplotパッケージのbaloonplot()関数を使っています。

そして、最後は、平均の誤差をグラフにしています。

reshapeパッケージのuntable()関数でfrequencyの数だけ実際の観測を作って、テーブルのデータを生データにしています。例えば、actual:3, estimated:3はfrequencyが23でしたので、actual:3, frequency:3の行を23行作っています。

そして、plyrパッケージのddply()関数でactualごとのerrorの平均を計算して、plot()関数でグラフにしたようです。

reshapeパッケージやplyrパッケージを使わなくてもできそうな気がします。

やってみましょう。

rep()関数とデータフレームの行を指定する方法を応用してuntable()と同じことができますし、tapply()関数がddply()関数と同じことができました。

今回はこれでおしまいです。

以下が今回のコードです。

# 2022-05-01
# HIstData_Jevons
#
library(HistData)
data("Jevons")
str(Jevons)
#
# show as tables
xtabs(frequency ~ estimated + actual, data = Jevons, subset = NULL)
xtabs(frequency ~ error + actual, data = Jevons, subset = NULL)
#
# show an sunflowerplot with regression line
with(Jevons,
     sunflowerplot(actual, estimated, frequency,
                   main = "Jevons data on numerical estimation"))
Jmod <- lm(estimated ~ actual, data = Jevons, weights = frequency)
abline(Jmod, col = "green", lwd = 2, lty = 2)
#
# show as balloonplots
library(gplots)
with(Jevons,
     balloonplot(actual, estimated, frequency, xlab = "actual",
                 ylab = "estimated",
                 main = "Jevons data on numerical estimation: Errros 
                 Bubble are area propotional to frequency",
                 text.size = 0.8))
with(Jevons,
     balloonplot(actual, error, frequency, xlab = "actual",
                 ylab = "error",
                 main = "Jevons data on numerical estimation:
                 Errors  Bubble are propotional to frequency",
                 text.size = 0.8))
#
# plot average error
library(reshape)
unJevons <- untable(Jevons, Jevons$frequency)
str(unJevons)
str(Jevons)
library(plyr)
mean_error <- function(df) mean(df$error, na.rm = TRUE)
Jmean <- ddply(unJevons, .(actual), mean_error)
with(Jmean,
     plot(actual, V1, ylab = "Mean error", xlab = "Actual number",
          type = "b", main = "Jevons data"))
abline(h = 0, col = "red", lty = 2, lwd = 2)
#
# reshapeパッケージ、plyrパッケージを使わずに再現する
n <- nrow(Jevons) # Jevonsの行数
index <- rep(c(1:n), Jevons$frequency) # Jevonsのx番目の行を何回繰り返すか
Jevons2 <- Jevons[index, ] # untable()と同じことをしている
Jmean2 <- with(Jevons2, tapply(error, actual, mean, na.rm = TRUE)) # actual
# ごとのerrorの平均値を計算
plot(Jmean2, type = "b", pch = 16, col = "blue", cex = 2,
     xlab = "Actual number", ylab = "Mean error",
     main = "Jevons data") # グラフ作成
abline(h = 0, col = "red", lty = 2, lwd = 2) # 0の水平線を追加