略解 (Rプログラム)

51

soft.th = function(lambda,x){
  sign(x) * pmax(abs(x) - lambda, 0 )
}

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 <- y-X[,-j]%*%beta[-j]
      beta[j] <- soft.th(lambda,sum(r*X[,j])/n)
    }
    eps <- max(abs(beta - beta.old))
    beta.old <- beta
  }
  beta.0 <- y.bar - sum(X.bar*beta)
  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, lasso(x,y,lambda)$beta[j])
  par(new = TRUE)
  lines(lambda.seq, coef.seq, col = j)
}

54

ridge <- 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) {
    beta <- drop(solve(t(X) %*% X + n * lambda * diag(p)) %*% t(X) %*% y)
    eps <- max(abs(beta - beta.old))
    beta.old <- beta
  }
  beta.0 <- y.bar - sum(X.bar*beta)
  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, ridge(x,y,lambda)$beta[j])
  par(new = TRUE)
  lines(lambda.seq, coef.seq, col = j)
}

55

library(glmnet)

## Loading required package: Matrix

## Loaded glmnet 4.1-4

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
lambda.min

## [1] 16.63647

fit <- glmnet(X, y, lambda = lambda.min)
fit$beta

## 5 x 1 sparse Matrix of class "dgCMatrix"
##           s0
## V3  9.933671
## V4 -2.906792
## V5  3.264911
## V6  .       
## V7  .