略解 (Rプログラム)

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 <- n * log(res$value / n) + 2 * k  
  if (AIC < AIC.min) {
    AIC.min <- AIC
    set.min <- res$set
  }
}
AIC.min

## [1] -1530.84

set.min

## [1]  1  3  4  6  8  9 10

47

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)
BIC.min <- Inf
for (k in 1:p) {
  T <- combn(1:p, k)
  res <- RSS.min(X, y, T)
  BIC <- n * log(res$value / n) + k * log(n)  
  if (BIC < BIC.min) {
    BIC.min <- AIC
    set.min <- res$set
  }
}
BIC.min

## [1] -1526.553

set.min

## [1] 10

AR2 <- function(X,y){
  p <- ncol(X)
  n <- length(y)
  y.bar <- mean(y)
  TSS <- sum((y - y.bar) ^ 2)
  D.max <- -Inf
  for (k in 1:p) {
    T <- combn(1:p, k)
    res <- RSS.min(X, y, T)
    D <- 1 - res$value / (n - k - 1) / TSS * (n - 1)
    if (D > D.max) {
      D.max <- D;
      set.max <- res$set
    }
  }
  print(D.max)
  set.max
}
AR2(X,y)

## [1] 0.9994328

## [1]  1  3  4  6  8  9 10

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 <- n * log(res$value / n) + 2 * k
  BIC <- n * log(res$value / n) + k * log(n)
  return(list(AIC = AIC, BIC = BIC))
}

AIC.seq <- NULL
BIC.seq <- NULL
for (k in 1:p) {
  AIC.seq <- c(AIC.seq, IC(k)$AIC)
  BIC.seq <- c(BIC.seq, IC(k)$BIC)
}
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)