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)