R Lasso regression (ラッソ回帰)

線形回帰(重回帰含む):http://upo-net.ouj.ac.jp/tokei/contents/sub_contents/c01_05_00.xml がわかりやすい。以下、そちらより抜粋。

library(lars)

library(glmnet)

 

data(diabetes)

attributes(diabetes$x)

attributes(diabetes$x2)

attributes(diabetes$y)

 

Linear <- lm(diabetes$y ~ diabetes$x)

summary(Linear)

 

LASSO.net  <- glmnet(diabetes$x, diabetes$y, alpha = 1)

plot(LASSO.net)

 

#乱数を発生させ、y=2x+εを満たす(x,y)を1000組構成

set.seed(1)

x <- rnorm(1000)

set.seed(2)

 

#x1はyとは無相関のダミー変数

x1 <- rnorm(1000)

y <- 2*x + rnorm(1000)

 

plot(x,y)

 

#線形回帰の最小2乗法を実施

LinearRegression <- function(beta){

    (y - beta[1] - beta[2]*x) %*% (y - beta[1] - beta[2]*x)

}

optim(c(0,0), LinearRegression)$par

lm(y ~ x)

 

#LASSO回帰の関数を作成

LassoRegression <- function(beta){

 

    #Pα(β)(ペナルティ項)の第2項(α=1なので、第1項は0となる)

    LassoRes <- sum(abs(beta))

 

    #線形回帰の最小2乗法部分

    RSS <- (y - beta[1] - beta[2]*x - beta[3]*x1) %*% (y - beta[1] - beta[2]*x - beta[3]*x1)/(2*length(x))

 

    #両者の合計で目的関数作成

    RSS + lambda * LassoRes

}

 

#LASSOなので、Pα(β)の第2項の係数(λ)は1

lambda <- 3

 

#LASSO回帰を実施。初期値をβ=(0, 0, 0)として最適化

optim(c(0, 0, 0), LassoRegression)$par

 

#グラフを描く際の変数を-5から10まで0.1刻みで生成

BetaX <- seq(-5, 10, by = 0.1)

 

#線形回帰式のy切片(β0)をBetaXの各点について計算

BetaInt <- mean(y) - BetaX*mean(x)

 

lambda <- 3

LassoRes <- lambda * (abs(BetaInt) + abs(BetaX))

 

#BetaXの点の数だけ値を格納するRSSを準備

RSS <- rep(NA, length(BetaX))

 

#誤差の2乗をBetaXの点の数だけ計算し、RSSに格納

for(i in 1:length(BetaX)) {

    RSS[i] <- (y - BetaInt[i] - BetaX[i]*x) %*% (y - BetaInt[i] - BetaX[i]*x)/(2*length(x))

}

 

#誤差の2乗項とペナルティ項の合計を最小化するようなLassoBetaを求値

LassoBeta <- BetaX[(RSS + LassoRes) == min(RSS + LassoRes)]

 

#BetaXの各点に対して、RSS項の値をプロット

plot(BetaX, RSS, type = "n", ylab = "objective function", main = paste("Lambda = ", lambda, sep = ""))

lines(BetaX, RSS)

 

#BetaXの各点に対して、ペナルティ項をプロット

lines(BetaX, LassoRes, col = "blue")

 

#BetaXの各点に対して、RSS項とペナルティ項の合計をプロット

lines(BetaX, RSS + LassoRes, col = "red")

 

abline(v = LassoBeta, lty = 2)

text(LassoBeta + 0.1, 30, paste("Beta = ", LassoBeta, sep = ""), adj = 0)

legend("bottomright", c("RSS", "LassoRes", "ObjectiveFunction"), col = c("black", "blue", "red"), lty = 1)