『スパース推定100問 with R』 問題コード

第3章　グループLasso

37

gr <- function(X, y, lambda) {
  nu <- 1 / max(2 * eigen(t(X) %*% X)$values)
  p <- ncol(X)
  beta <- rnorm(p)
  beta.old <- rnorm(p)
  while (max(abs(beta - beta.old)) > 0.001) {
    beta.old <- beta
    gamma <- ## 空欄(1) ##
    beta <- ## 空欄(2) ##
  }
  return(beta)
}

## データの生成
n <- 100
p <- 3
X <- matrix(rnorm(n * p), ncol = p)
beta <- rnorm(p)
epsilon <- rnorm(n)
y <- 0.1 * X %*% beta + epsilon

## 係数の値の変化を表示
lambda <- seq(1, 50, 0.5)
m <- length(lambda)
beta <- matrix(nrow = m, ncol = p)
for (i in 1:m) {
  est <- gr(X, y, lambda[i]) 
  for (j in 1:p)
    beta[i, j] <- est[j]
}
y.max <- max(beta)
y.min <- min(beta)
plot(lambda[1]:lambda[m], ylim = c(y.min, y.max),
     xlab = "lambda", ylab = "係数の値", type = "n")
for (j in 1:p)
  lines(lambda, beta[, j], col = j + 1)
legend("topright", legend = paste("係数", 1:p), lwd = 2, col = 2:(p + 1))
segments(lambda[1], 0, lambda[m], 0)

38

n <- 100
p <- 1  # p <- 3 
X <- matrix(rnorm(n * p), ncol = p)
beta <- rnorm(p)
epsilon <- rnorm(n)
y <- 0.1 * X %*% beta + epsilon
lambda <- 0.01
nu <- 1 / max(2 * eigen(t(X) %*% X)$values)
p <- ncol(X)
m <- 10
beta <- rep(1, p)
beta.old <- rep(0, p)
t <- 0
val <- matrix(0, m, p)
while (t < m) {
  t <- t + 1
  val[t, ] <- beta
  beta.old <- beta
  gamma <- ## 空欄(1) ##
  beta <- ## 空欄(2) ##
}
eval <- array(dim = m) 
val.final <- val[m, ]
for (i in 1:m)
  eval[i] <- norm(val[i, ] - val.final, "2")
plot(1:m, ylim = c(0, eval[1]), type = "n",
     xlab = "回数", ylab = "L2 誤差", main = "FISTA と ISTA の比較")
lines(eval, col = "blue")
beta <- rep(1, p)
beta.old <- rep(0, p)
alpha <- 1
gamma <- beta
t <- 0
val <- matrix(0, m, p)
while (t < m) {
  t <- t + 1
  val[t, ] <- beta
  beta.old <- beta
  w <- ## 空欄(3) ##
  beta <- ## 空欄(4) ##
  alpha.old <- alpha
  alpha <- ## 空欄(5) ##
  gamma <- ## 空欄(6) ##
}
val.final <- val[m, ]
for (i in 1:m)
  eval[i] <- norm(val[i, ] - val.final, "2")
lines(eval, col = "red")
legend("topright", c("FISTA", "ISTA"), lwd = 1, col = c("red", "blue"))

41(b)

fista <- function(X, y, lambda) {
  nu <- 1 / max(2 * eigen(t(X) %*% X)$values)
  p <- ncol(X)
  alpha <- 1
  beta <- rnorm(p)
  beta.old <- rnorm(p)
  gamma <- beta
  while (max(abs(beta - beta.old)) > 0.001) {
    print(beta)
    beta.old <- beta
    w <- ## 空欄(1) ##
    beta <- ## 空欄(2) ##
    alpha.old <- alpha
    alpha <- (1 + sqrt(1 + 4 * alpha ^ 2)) / 2
    gamma <- beta + (alpha.old - 1) / alpha * (beta - beta.old)
  }
  return(beta)
}

42

group.lasso <- function(z, y, lambda = 0) {
  J <- length(z)
  theta <- list()
  for (j in 1:J)
    theta[[j]] <- rep(0, ncol(z[[j]]))
  for (m in 1:10) {
    for (j in 1:J) {
      r <- y
      for (k in 1:J) {
        if (k != j)
          r <- r - ## 空欄(1) ##
      }
      theta[[j]] <- ## 空欄(2) ##
    }
  }
  return(theta)
}
  
## データの生成
N <- 100
J <- 2
u <- rnorm(n)
v <- u + rnorm(n)
s <- 0.1 * rnorm(n)
t <- 0.1 * s + rnorm(n)
y <- u + v + s + t + rnorm(n)
z <- list()
z[[1]] <- cbind(u, v)
z[[2]] <- cbind(s, t)  

## 係数の値の変化を表示
lambda <- seq(1, 500, 10)
m <- length(lambda)
beta <- matrix(nrow = m, ncol = 4)
for (i in 1:m) {
  est <- group.lasso(z, y, lambda[i])
  beta[i, ] <- c(est[[1]][1], est[[1]][2], est[[2]][1], est[[2]][2])
}
y.max <- max(beta)
y.min <- min(beta)
plot(lambda[1]:lambda[m], ylim = c(y.min, y.max),
     xlab = "lambda", ylab = "係数の値", type = "n")
lines(lambda, beta[, 1], lty = 1, col = 2)
lines(lambda, beta[, 2], lty = 2, col = 2)
lines(lambda, beta[, 3], lty = 1, col = 4)
lines(lambda, beta[, 4], lty = 2, col = 4)
legend("topright",
       legend = c("グループ1", "グループ1", "グループ2", "グループ2"),
       lwd = 1, lty = c(1, 2), col = c(2, 2, 4, 4))
segments(lambda[1], 0, lambda[m], 0)

44

sparse.group.lasso <- function(z, y, lambda = 0, alpha = 0) {
  J <- length(z)
  theta <- list()
  for (j in 1:J)
    theta[[j]] <- rep(0, ncol(z[[j]]))
  for (m in 1:10) {
    for (j in 1:J) {
      r <- y
      for (k in 1:J) {
        if (k != j)
          r <- r - z[[k]] %*% theta[[k]]
      }
      theta[[j]] <- ## 空欄(1) ## 
    }
  }
  return(theta)
}

sparse.gr <- function(X, y, lambda, alpha = 0) {
  nu <- 1 / max(2 * eigen(t(X) %*% X)$values)
  p <- ncol(X)
  beta <- rnorm(p)
  beta.old <- rnorm(p)
  while (max(abs(beta - beta.old)) > 0.001) {
    beta.old <- beta
    gamma <- beta + nu * t(X) %*% (y - X %*% beta)
    delta <- ## 空欄(2) ##
    beta <- ## 空欄(3) ##
  }
  return(beta)
}

46

## データの生成
n <- 100
J <- 2
x <- rnorm(n)
y <- x + cos(x)
z[[1]] <- cbind(rep(1, n), x)
z[[2]] <- cbind(## 空欄(1) ##)

## 係数の値の変化を表示
lambda <- seq(1, 200, 5)
m <- length(lambda)
beta <- matrix(nrow = m, ncol = 5)
for (i in 1:m) {
  est <- ## 空欄(2) ##
  beta[i, ] <- c(est[[1]][1], est[[1]][2], est[[2]][1], est[[2]][2], est[[2]][3])
}
y.max <- max(beta)
y.min <- min(beta)
plot(lambda[1]:lambda[m], ylim = c(y.min, y.max),
     xlab = "lambda", ylab = "係数の値", type = "n")
lines(lambda, beta[, 1], lty = 1, col = 2)
lines(lambda, beta[, 2], lty = 2, col = 2)
lines(lambda, beta[, 3], lty = 1, col = 4)
lines(lambda, beta[, 4], lty = 2, col = 4)
lines(lambda, beta[, 5], lty = 3, col = 4)
legend("topright", legend = c("1", "x", "cos x", "cos 2x", "cos 3x"),
       lwd = 1, lty = c(1, 2, 1, 2, 3), col = c(2, 2, 4, 4, 4))
segments(lambda[1], 0, lambda[m], 0)

i <- 5  # lambda[5]の値を用いた

f.1 <- function(x) {
  return(beta[i, 1] + beta[i, 2] * x)
}

f.2 <- function(x) {
  return(beta[i, 3] * cos(x) + beta[i, 4] * cos(2 * x) + beta[i, 5] * cos(3 * x))
}

f <- function(x) {
  return(f.1(x) + f.2(x))
}

curve(f.1(x), -5, 5, col = "red", ylab = "関数の値")
curve(f.2(x), -5, 5, col = "blue", add = TRUE)
curve(f(x), -5, 5, add = TRUE)
legend("topleft", legend = c("f = f.1 + f.2", "f.1", "f.2"),
       col = c(1, "red", "blue"), lwd = 1)