第1章　線形回帰

1

inner.prod <- function(x, y) {
  return(sum(x * y))
}

linear <- function(X, y) {
  n <- nrow(X)
  p <- ncol(X)
  X <- as.matrix(X)
  x.bar <- array(dim = p)
  for (j in 1:p)
    x.bar[j] <- mean(X[, j])  
  for (j in 1:p)
    X[, j] <- X[, j] - x.bar[j]  # Xの中心化
  y <- as.vector(y)
  y.bar <- mean(y)
  y <- y-y.bar                   # yの中心化
  beta <- as.vector(## 空欄(1) ##)
  beta.0 <- ## 空欄(2) ##
  return(list(beta = beta, beta.0 = beta.0))
}

7(b)

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

8

linear.lasso <- function(X, y, lambda = 0) {
  X <- as.matrix(X)
  n <- nrow(X)
  p <- ncol(X)
  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
  scale <- array(dim = p)
  for (j in 1:p) {
    scale[j] <- sqrt(sum(X[, j] ^ 2) / n)
    X[, j] <- X[, j] / scale[j]
  }
  eps <- 1
  beta <- rep(0, p)
  beta.old <- rep(0, p)
  while (eps > 0.001) {
    for (j in 1:p) {
      r <- ## 空欄(1) ##
      beta[j] <- soft.th(lambda, sum(r * X[, j]) / n) / (sum(X[, j] * X[, j]) / n) 
    }
    eps <- max(abs(beta - beta.old))
    beta.old <- beta
  }
  for (j in 1:p) 
    beta[j] <- beta[j] / scale[j]
  beta.0 <- ## 空欄(2) ##
  return(list(beta = beta, beta.0 = beta.0))
}

crime <- read.table("crime.txt")
X <- crime[, 3:7]
y <- crime[, 1] 
linear.lasso(X, y, 10)
linear.lasso(X, y, 50)
linear.lasso(X, y, 100)

9

library(glmnet)

df <- read.table("crime.txt")
x <- as.matrix(df[, 3:7])
y <- df[, 1]
fit <- glmnet(x, y)
plot(fit, xvar = "lambda")

10

warm.start <- function(X, y, lambda.max = 100) {
  dec <- round(lambda.max / 50)
  lambda.seq <- seq(lambda.max, 1, -dec)
  r <- length(lambda.seq)
  p <- ncol(X)
  coef.seq <- matrix(nrow = r, ncol = p)
  coef.seq[1, ] <- linear.lasso(X, y, lambda.seq[1])$beta 
  for (k in 2:r) {
    coef.seq[k, ] <- ## 空欄 ##
  }
  return(coef.seq)
}

crime <- read.table("crime.txt")
X <- crime[, 3:7]
y <- crime[, 1]
coef.seq <- warm.start(X, y, 200)
p <- ncol(X)
lambda.max <- 200
dec <- round(lambda.max / 50)
lambda.seq <- seq(lambda.max, 1, -dec)
plot(log(lambda.seq), coef.seq[, 1],
     xlab = "log(lambda)", ylab = "係数",
     ylim = c(min(coef.seq), max(coef.seq)), type="n")
for (j in 1:p) 
  lines(log(lambda.seq), coef.seq[, j], col = j)
X <- as.matrix(X)
y <- as.vector(y)
cv <- cv.glmnet(X, y)
plot(cv)

13

crime <- read.table("crime.txt")
X <- crime[, 3:7]
y <- crime[, 1]
linear(X, y)
ridge(X, y)
ridge(X, y, 200)

14

df <- read.table("crime.txt")
x <- df[, 3:7]
y <- df[, 1]
p <- ncol(x)
lambda.max <- 3000
lambda.seq <- seq(1, lambda.max)
plot(lambda.seq,
     xlim = c(0, lambda.max), ylim = c(-12, 12),
     xlab = "lambda", ylab = "係数",        ## この部分
     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)      ## この部分
}
legend("topright",
       legend=c("警察への年間資金",
                "25 歳以上で高校を卒業した人の割合",
                "16-19 歳で高校に通っていない人の割合",
                "18-24 歳で大学生の割合",
                "25 歳以上で 4 年制大学を卒業した人の割合"),
       col = 1:p, lwd = 2, cex = .8)

17

n <- 500
x <- array(dim = c(n, 6))
z <- array(dim = c(n, 2))
for (i in 1:2)
  z[, i] <- rnorm(n)
y <- ## 空欄(1) ##
for (j in 1:3)
  x[, j] <- z[, 1] + rnorm(n) / 5
for (j in 4:6)
  x[, j] <- z[, 2] + rnorm(n) / 5
glm.fit <- glmnet(## 空欄(2) ##)
plot(glm.fit)

20

alpha <- seq(0.01, 0.99, 0.01)
m <- length(alpha)
mse <- array(dim = m)
for (i in 1:m){
  cvg <- cv.glmnet(x, y, alpha = alpha[i])
  mse[i] <- min(cvg$cvm)
}
best.alpha <- alpha[which.min(mse)]
best.alpha
cva <- cv.glmnet(x, y, alpha = best.alpha)
best.lambda <- cva$lambda.min
best.lambda