情報量統計学（共立出版） - R入門

はじめに

「情報量統計学」の計算を再現します。

目次

第8章 回帰モデル

pp.133-135

データは、p.128の表8.1に示されている。以下のとおり。これをカレントディレクトリにtable8_1.csvというファイルで保存しておく。

i,  xi,     yi
 1, 0.0,  0.012
 2, 0.1,  0.121
 3, 0.2, -0.097
 4, 0.3, -0.061
 5, 0.4, -0.080
 6, 0.5,  0.037
 7, 0.6,  0.196
 8, 0.7,  0.077
 9, 0.8,  0.343
10, 0.9,  0.448
11, 1.0,  0.434

ファイルからデータの読み込み、AICを計算する関数を作成。予備計算。

> dtf <- read.csv("table8_1.csv", header = TRUE)
> xi <- dtf$xi
> yi <- dtf$yi
> daic <- function(n, m, d) return(n * log(2 * pi) + n * log(d) + n + 2 * (m + 2))
> n <- nrow(dtf)
> ssx <- sum(xi)
> ssx2 <- sum(xi ^ 2)
> ssy <- sum(yi)
> ssy2 <- sum(yi ^ 2)
> ssxy <- sum(xi * yi)

平均0の正規分布を計算。書籍では計算結果を途中で丸めているため（例えば以下の例では書籍ではd(-1)=0.0533だが、実際には0.05333981818…）、本計算ではそのようにしていないため、計算結果はぴったり一致しない。以下、同じ。

> m <- -1
> d <- ssy2 / n
> cat(sprintf("m = %d, d = %f, AIC = %f\n", m, d, daic(n, m, d)))
m = -1, d = 0.053340, AIC = 0.974854

平均a0の正規分布を計算。

> m <- 0
> a0 <- ssy / n
> d <- (ssy2 - a0 * ssy) / n
> cat(sprintf("m = %d, d = %f, AIC = %f\n", m, d, daic(n, m, d)))
m = 0, d = 0.036440, AIC = -1.216377

直線回帰モデルを計算。書籍（初版第10刷、p.134）ではAICの計算結果は10.31となっているが、誤植と思われる（p.135の表8.2では-10.31となっている）。

> m <- 1
> a0 <- (ssx2 * ssy - ssx * ssxy) / (n * ssx2 - ssx ^ 2)
> a1 <- (n * ssxy - ssx * ssy) / (n * ssx2 - ssx ^ 2)
> d <- (ssy2 - a0 * ssy - a1 * ssxy) / n
> cat(sprintf("m = %d, d = %f, AIC = %f\n", m, d, daic(n, m, d)))
m = 1, d = 0.013312, AIC = -10.292953

2次の多項式回帰モデルを計算。

> m <- 2
> c1 <- c(n, sum(xi), sum(xi ^ 2))
> c2 <- c(sum(xi), sum(xi ^ 2), sum(xi ^ 3))
> c3 <- c(sum(xi ^ 2), sum(xi ^ 3), sum(xi ^ 4))
> mxxx <- matrix(c(c1, c2, c3), ncol = 3, byrow = TRUE)
> mxy <- matrix(c(ssy, ssxy, sum(xi ^ 2 * yi)), ncol = 1)
> mxa <- solve(t(mxxx) %*% mxxx) %*% mxxx %*% mxy
> d <- (ssy2 - sum(t(mxa) %*% mxy)) / n
> cat(sprintf("m = %d, d = %f, AIC = %f\n", m, d, daic(n, m, d)))
m = 2, d = 0.005929, AIC = -17.191032

最後に、表8.2（p.135）を作成する。自由パラメーターが1の場合は計算が特殊のため、最初に別途計算。あとはfor文を使用して繰り返し計算。

> m = -1
> d <- ssy2 / n
> cat(sprintf("m = %d, d = %.5f AIC = %.2f\n", m, d, daic(n, m, d)))
m = -1, d = 0.05334 AIC = 0.97
> mxy <- matrix(yi, ncol = 1)
> for (i in 0:5) {
+ mxxx <- matrix(0., nrow = n, ncol = (i + 1))
+ for (j in 1:n) {
+ mxxx[j, ] <- xi[j] ^ (0:i)
+ }
+ mxa <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
+ d <- sum((mxy - mxxx %*% mxa) ^ 2) / n
+ cat(sprintf("m = %d, d = %.5f AIC = %.2f\n", i, d, daic(n, i, d)))
+ }
m = 0, d = 0.03644 AIC = -1.22
m = 1, d = 0.01331 AIC = -10.29
m = 2, d = 0.00593 AIC = -17.19
m = 3, d = 0.00516 AIC = -16.72
m = 4, d = 0.00440 AIC = -16.46
m = 5, d = 0.00425 AIC = -14.86

---