第5章　正則化（問題49～56）

51

curve(soft.th(5, x), -10, 10)

53

lasso <- function(X, y, lambda = 0) {
  X <- as.matrix(X)
  X <- scale(X)
  p <- ncol(X)
  n <- length(y)
  X.bar <- array(dim = p)
  for (j in 1:p) {
    X.bar[j] <- mean(X[, j])
    X[, j] <- X[, j] - X.bar[j]
  }
  y.bar <- mean(y)
  y <- y - y.bar
  eps <- 1
  beta <- array(0, dim = p)
  beta.old <- array(0, dim = p)
  while (eps > 0.001) {
    for (j in 1:p) {
      r <- ## 空欄(1) ##
      beta[j] <- ## 空欄(2) ##
    }
    eps <- max(abs(beta - beta.old))
    beta.old <- beta
  }
  beta.0 <- ## 空欄(3) ##
  return(list(beta = beta, beta.0 = beta.0))
}

df <- read.table("crime.txt")
x <- df[, 3:7]
y <- df[, 1]
p <- ncol(x)
lambda.seq <- seq(0, 100, 0.1)
coef.seq <- lambda.seq
plot(lambda.seq, coef.seq, xlim = c(0, 100), ylim = c(-40, 40),
     xlab = "lambda", ylab = "beta", main = "各lambdaについての各係数の値", 
     type = "n", col = "red")
for (j in 1:p) {
  coef.seq <- NULL
  for (lambda in lambda.seq)
    coef.seq <- c(coef.seq, ## 空欄(4) ##)
  par(new = TRUE)
  lines(lambda.seq, coef.seq, col = j)
}

54

beta <- drop(solve(t(X) %*% X + n * lambda * diag(p)) %*% t(X) %*% y)

55

library(glmnet)

df <- read.table("crime.txt")
X <- as.matrix(df[, 3:7])
y <- df[,1]
cv.fit <- cv.glmnet(X, y)
lambda.min <- cv.fit$lambda.min
fit <- glmnet(X, y, lambda = lambda.min)