応用統計学・統計解析法レポート問７

main.R

#!/usr/bin/env Rscript

# 統計データがなかったので、お正月番組に触発されてテーブルクロス引きでコップが
# ずれた距離を測りました。コップの中の水の容量は
#   (x1, x2,..., x5) = (0, 50,.., 200)
# です。それぞれの容量で５回ずつ測定。単位はcmとccです。

x1 <- c(1.5, 2.2, 2.4, 2.2, 2.7)
x2 <- c(1.8, 1.6, 1.9, 1.6, 1.6)
x3 <- c(1.7, 2.1, 1.6, 1.5, 2.0)
x4 <- c(1.7, 2.3, 1.6, 1.9, 1.8)
x5 <- c(1.6, 2.0, 2.1, 1.8, 1.7)

mu_hat <- mean(c(mean(x1), mean(x2), mean(x3), mean(x4), mean(x5)))

st <- sum(c(sum((x1 - mu_hat) ^ 2), sum((x2 - mu_hat) ^ 2), sum((x3 - mu_hat) ^ 2), sum((x4 - mu_hat) ^ 2), sum((x5 - mu_hat) ^ 2)))
se <- sum(c(sum((x1 - mean(x1)) ^ 2), sum((x2 - mean(x2)) ^ 2), sum((x3 - mean(x3)) ^ 2), sum((x4 - mean(x4)) ^ 2), sum((x5 - mean(x5)) ^ 2)))

# モデル１
# 水の容量は無関係ですべて同じ
# Xij = mu + Eij, Eij ~ N(mu, sigma)
aic1 <- 5 * 5 * (log(st / (5 * 5)) + log(pi * exp(1))) + 2 * 2
aic1

# モデル２
# 水の容量に依る
# Xij = mu + Aij + Eij, Eij ~ N(mu, sigma)
aic2 <- 5 * 5 * (log(se / (5 * 5)) + log(pi * exp(1))) + 2 * (5 + 1)
aic2

# aic1 > aic2なのでaic2の方が優れたモデル。つまり、テーブルクロス引きでコップの
# 移動量は水の容量に依存する。

# 水の容量が0ccの時だけ値が外れている気がしたのでx1を外して同様に解析する。

mu_hat <- mean(c(mean(x2), mean(x3), mean(x4), mean(x5)))

st <- sum(c(sum((x2 - mu_hat) ^ 2), sum((x3 - mu_hat) ^ 2), sum((x4 - mu_hat) ^ 2), sum((x5 - mu_hat) ^ 2)))
se <- sum(c(sum((x2 - mean(x2)) ^ 2), sum((x3 - mean(x3)) ^ 2), sum((x4 - mean(x4)) ^ 2), sum((x5 - mean(x5)) ^ 2)))

aic1 <- 5 * 4 * (log(st / (5 * 4)) + log(pi * exp(1))) + 2 * 2
aic1

aic2 <- 5 * 4 * (log(se / (5 * 4)) + log(pi * exp(1))) + 2 * (4 + 1)
aic2

# モデル１の方がAICが小さくなる。よって、水の容量が一定を超えた後は水の容量に
# 依らなくなる。