第4章　情報量基準（問題40～48）

46

RSS.min <- function(X, y, T) {
  m <- ncol(T)
  S.min <- Inf
  for (j in 1:m) {
    q <- T[, j]
    S <- sum((lm(y ~ X[, q])$fitted.values - y) ^ 2) / n
    if (S < S.min) {
      S.min <- S
      set.q <- q
    }
  }
  return(list(value = S.min, set = set.q))
}

library(MASS)

df <- Boston
X <- as.matrix(df[, c(1, 3, 5, 6, 7, 8, 10, 11, 12, 13)])
y <- df[[14]]
p <- ncol(X)
n <- length(y)
AIC.min <- Inf
for (k in 1:p) {
  T <- combn(1:p, k)
  res <- RSS.min(X, y, T)
  AIC <- ## 空欄(1) ##
  if (AIC < AIC.min) {
    AIC.min <- ## 空欄(2) ##
    set.min <- ## 空欄(3) ##
  }
}
AIC.min
set.min

48

library(MASS)

df <- Boston
X <- as.matrix(df[, c(1, 3, 5, 6, 7, 8, 10, 11, 12, 13)])
y <- df[[14]]
n <- nrow(X)
p <- ncol(X)

IC <- function(k) {
  T <- combn(1:p, k)
  res <- RSS.min(X, y, T)
  AIC <- ## 空欄(1) ##
  BIC <- ## 空欄(2) ##
  return(list(AIC = AIC, BIC = BIC))
}

AIC.seq <- NULL
BIC.seq <- NULL
for (k in 1:p) {
  AIC.seq <- c(AIC.seq, ## 空欄(3) ##)
  BIC.seq <- c(BIC.seq, ## 空欄(4) ##)
}
plot(1:p, ylim = c(min(AIC.seq), max(BIC.seq)), type = "n",
     xlab = "変数の個数", ylab = "AIC/BICの値")
lines(AIC.seq, col = "red")
lines(BIC.seq, col = "blue")
legend("topleft", legend = c("AIC", "BIC"), col = c("red", "blue"),
       lwd = 1, cex = .8)