「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。

MENU

統計的機械学習の数理100問(14/20)

 いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として

を用いることにする。

5. 正則化

  • 標本数nが変数の数[etx:p]よりも小さい場合、最小二乗法など通常の手段は用いることができない。このときに\mathrm{Lasso}または\mathrm{Ridge}を用いる。
  • 線形回帰において最小化する目的関数に係数の値が大きくならなくするための正則化項を加える。この正則化項を\mathcal{L}^1-ノルムの定数\lambda倍にしたものを\mathrm{Lasso}\mathcal{L}^2-ノルムの定数\lambda倍にしたものを\mathrm{Ridge}という。
  • \mathrm{Lasso}では\lambdaが大きくなるにつれて係数の推定値が0になる場合が出てきて、\lambda\rightarrow\inftyではすべての係数の推定値が0になる。このためモデル選択も担う。

 以下、線形回帰を考えるが、一般性を失うことなく切片項\beta_0=0とし、標本数をn、変数の数をpとして\boldsymbol{X}\in\mathbb{R}^{n\times p}とする。

5.1 Ridge

 行列{}^{t}\boldsymbol{X}\boldsymbol{X}が特異ではないとも、その行列式が小さい場合、信頼区間が大きくなるなど推定値の信頼性が下がる。そこで定数\lambda\geq0を用いて


\begin{aligned}
L=\displaystyle{\frac{1}{n}\|\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}\|^2}+\lambda\|\boldsymbol{\beta}\|_2^2
\end{aligned}

を最小化することで母数(ベクトル)\boldsymbol{\beta}を推定する*1。これを\mathrm{Ridge}という。
 目的関数L\boldsymbol{\beta}微分すると


\begin{aligned}
\displaystyle{\frac{\partial L}{\partial\boldsymbol{\beta}}}=-\displaystyle{\frac{2}{n}{}^{t}\boldsymbol{X}\left(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}\right)+2\lambda\boldsymbol{\beta}}
\end{aligned}

を得るから、{}^{t}\boldsymbol{X}\boldsymbol{X}+n\lambda Iが正則ならば、


\begin{aligned}
&\left(\displaystyle{\frac{2}{n}}{}^{t}\boldsymbol{X}\boldsymbol{X}+2\lambda I\right)\boldsymbol{\beta}=\displaystyle{\frac{2}{n}}{}^{t}\boldsymbol{X}\boldsymbol{Y},\\
\Leftrightarrow&\boldsymbol{\beta}=\left({}^{t}\boldsymbol{X}\boldsymbol{X}+n\lambda I\right)^{-1}\boldsymbol{X}\boldsymbol{Y}
\end{aligned}

が得られる。
 もし\lambda\gt0ならば{}^{t}\boldsymbol{X}\boldsymbol{X}が非負定値であることからその固有値p_1,\cdots,p_pはすべて非負である。したがってt\in\mathbb{R}に対して


\begin{aligned}
&\mathrm{det}\left({}^{t}\boldsymbol{X}\boldsymbol{X}+n\lambda I-tI\right)=0\\
\Leftrightarrow&t=\mu_1-n\lambda,\cdots,\mu_p+n\lambda\gt0
\end{aligned}

と、すべての固有値は正である。このため{}^{t}\boldsymbol{X}\boldsymbol{X}+n\lambda Iは正則であることが保証される。これはn,pの大小に依存しない。もしn\lt pならば{}^{t}\boldsymbol{X}\boldsymbol{X}は階数がn未満で特異である。したがって
 


\begin{aligned}
\lambda\gt\Leftrightarrow {}^{t}\boldsymbol{X}\boldsymbol{X}+n\lambda Iは正則
\end{aligned}

である。

5.2 劣微分

 f(x)=x^2+x+2|x|のように絶対値が入ることで微分可能でない点をもつ関数にも微分をできるように概念を拡張できれば、最適化をより柔軟に適用し得る。そこで微分の概念を拡張する。
 前提としてfが凸関数であることを仮定する。


微分 凸関数f:\mathbb{R}\rightarrow\mathbb{R}およびx_0\in\mathbb{R}{}^{\forall}x\in\mathbb{R}について


\begin{aligned}
f(x)\geq f(x_0)+z(x-x_0)
\end{aligned}

を満たすようなz\in\mathbb{R}の集合をfx=x_0における劣微分という。

 f微分可能な点においては対応する劣微分の値は微分の値と一致する。

5.3 Lasso

 定数\lambda\geq0を用いて


\begin{aligned}
L=\displaystyle{\frac{1}{2n}\|\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}\|^2}+\lambda\|\boldsymbol{\beta}\|^2
\end{aligned}

を最小化することで母数(ベクトル)\boldsymbol{\beta}を推定する方法を\mathrm{Lasso}という。
 \mathrm{Lasso}の解析には劣微分を用いる。簡単のため、


\begin{aligned}
\displaystyle{\frac{1}{n}\sum_{i=1}^{n}x_{i,j}x_{i,k}}=\begin{cases}
1,&j=k,\\
0,&j\neq k
\end{cases}
\end{aligned}

を仮定し、


\begin{aligned}
s_j:=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}x_{i,j}y_i}
\end{aligned}

とおく。
 このときL\beta_jに関する劣微分*2を求めれば、|\beta_j|\beta_j=0での劣微分[-1,1]であるから、


\begin{aligned}
\displaystyle{\frac{\partial L}{\partial\beta_j }}&=\begin{cases}
\displaystyle{-\frac{1}{n}\sum_{i=1}^{n}x_{i,j}\left(y_i-\sum_{k=1}^{p}x_{i,k}\beta_k\right)}+\lambda,&\beta_j\gt0,\\
\displaystyle{-\frac{1}{n}\sum_{i=1}^{n}x_{i,j}\left(y_i-\sum_{k=1}^{p}x_{i,k}\beta_k\right)}-\lambda,&\beta_j\lt0,\\
\displaystyle{-\frac{1}{n}\sum_{i=1}^{n}x_{i,j}\left(y_i-\sum_{k=1}^{p}x_{i,k}\beta_k\right)}+\lambda a,&\beta_j=0,{}^{\forall}a\in[-1,1]
\end{cases}\\
&=\begin{cases}
\ -s_j+\beta_j+\lambda,&\beta_j\gt0,\\
\ -s_j+\beta_j-\lambda,&\beta_j\lt0,\\
\ -s_j+\beta_j+\lambda a,&\beta_j=0,{}^{\forall}a\in[-1,1]
\end{cases}
\end{aligned}

である。したがって


\begin{aligned}
\beta_j=\begin{cases}
s_j-\lambda,&s_j\gt\lambda,\\
s_j+\lambda,&s_j\lt-\lambda,\\
0,&-\lambda\leq s_j\leq\lambda
\end{cases}
\end{aligned}

が得られる。
 もし規格化の仮定が満たされない場合は、


\begin{aligned}
\displaystyle{\frac{\partial L}{\partial\beta_j }}&=\begin{cases}
\displaystyle{-\frac{1}{n}\sum_{i=1}^{n}x_{i,j}\left(y_i-\sum_{k=1}^{p}x_{i,k}\beta_k\right)}+\lambda,&\beta_j\gt0,\\
\displaystyle{-\frac{1}{n}\sum_{i=1}^{n}x_{i,j}\left(y_i-\sum_{k=1}^{p}x_{i,k}\beta_k\right)}-\lambda,&\beta_j\lt0,\\
\displaystyle{-\frac{1}{n}\sum_{i=1}^{n}x_{i,j}\left(y_i-\sum_{k=1}^{p}x_{i,k}\beta_k\right)}+\lambda a,&\beta_j=0,{}^{\forall}a\in[-1,1]
\end{cases}
\end{aligned}

において、\beta_k,k\neq jを固定して\beta_jを更新する。これをj=1,\cdots,pまで何度も繰り返して収束させる。そして各段階で\beta_jを更新する際にy_iをその残差r_{i,j}=y_i-\displaystyle{\sum_{k\neq j}x_{i,k}\beta_k}で置き換える。また


\begin{aligned}
\begin{cases}
\displaystyle{-\frac{1}{n}\sum_{i=1}^{n}x_{i,j}(r_{i,j}-x_{i,j}\beta_j)}+\lambda,&\beta_j\gt0,\\
\displaystyle{-\frac{1}{n}\sum_{i=1}^{n}x_{i,j}(r_{i,j}-x_{i,j}\beta_j)}-\lambda,&\beta_j\lt0,\\
\displaystyle{-\frac{1}{n}\sum_{i=1}^{n}x_{i,j}(r_{i,j}-x_{i,j}\beta_j)}+\lambda a,&\beta_j=0,{}^{\forall}a\in[-1,1]
\end{cases}
\end{aligned}

においてs_j=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}r_{i,j}x_{i,j}}として、\beta_j=S_{\lambda}(s_j)を更新する。

5.4 RidgeとLassoの比較

 \mathrm{Ridge}および\mathrm{Lasso}はいずれも\lambdaが大きくなるにつれて係数の絶対値は減少し0に近づいていく。\mathrm{Lasso}\lambdaが一定値以上になると係数値が0となり、そのタイミングが変数ごとに相違するために変数選択に用いることができる。このように\mathrm{Lasso}は変数選択に用いることができる。
 簡単のためp=2とすれば、この最適化問題は、中心が(\hat{\beta}_1,\hat{\beta}_2)であるような


\begin{aligned}
\displaystyle{\frac{1}{n}\sum_{i=1}^{n}x_{i,j}x_{i,k}}=\begin{cases}
1,&j=k,\\
0,&j\neq k
\end{cases}
\end{aligned}

の値が等しい点の集合(楕円)と正則化項(\mathcal{L}^1正則化|\beta_1|+|\beta_2|\leq C,\mathcal{L}^2正則化\beta_1^2+\beta_2^2\leq C)との交点を求める問題に等しい。


LassoおよびRidgeの幾何学的イメージ(左:Lasso/右:Ridge)

     

5.5 λの設定方法

 \lambdaはクロスバリデーションを適用して決定するのが多い。Rの場合、\mathrm{glmnet}パッケージを用いると便利である。

###############
### λの設定 ###
###############

library("glmnet")

df <- read.table(file = "C:/Users/Julie/Downloads/crime.txt",header = F,sep = ",")

X <- as.matrix(df[,3:ncol(df)])
y <- df[,1]

cv_fit <- cv.glmnet(X,y) # クロスバリデーションによるλ設定

plot(cv_fit)


lambda_min <- cv_fit$lambda.min

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

5.6 Rでのシミュレーション

5.6.1 Ridge
#################
### Ridge回帰 ###
#################

library("ggplot2")

# Ridge
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(i in 1:p){
    X_bar[i] <- mean(X[,i])
    X[,i] <- X[,i] - X_bar[i]
  }
  
  y_bar <- mean(y)
  y <- y - y_bar
  
  beta <- drop(solve(t(X) %*% X + n* lambda * diag(p)) %*% t(X) %*% y)
  beta_0 <- y_bar - sum(X_bar * beta)
  
  return(list(beta = beta, beta_0 = beta_0))
}


# シミュレーション
vc_names <- c("警察への年間資金","25歳以上で高卒人数の割合","16-19歳で高校に通っていない人の割合",
              "18-24歳で大学生の割合","25歳以上で4年制大卒の割合")

df <- read.table(file = "(ファイルパス)",header = F,sep = ",")

X <- df[,3:ncol(df)]
y <- df[,1]

p <- ncol(X)

lambda_seq <- seq(0, 100, 0.1)
coef_seq <- lambda_seq

df_betas <- data.frame(matrix(NA,nrow = p, ncol = length(lambda_seq)))
colnames(df_betas) <- 1:ncol(df_betas)

# Ridge推定
for(lambda in lambda_seq){
  df_betas[,which(lambda_seq==lambda)] <- ridge(X,y,lambda)$beta
}

df_betas_graph <- NULL

for(i in 1:nrow(df_betas)){
  df_betas_graph <- rbind(df_betas_graph,
                          data.frame(lambda = lambda_seq,
                                     var = vc_names[i],
                                     estimate = unlist(df_betas[i,]))
  )
}

# 図示
g <- ggplot(df_betas_graph,aes(x = lambda,y = estimate, group = var, color = var)) + geom_line()
g <- g + xlab("λ") + ylab("推定値")+ labs(title = "Ridge回帰による各母数の推定値",color = "")
g <- g + theme_classic() + theme(plot.title = element_text(hjust = 0.5),
                                 legend.position = c(1, 1),
                                 legend.justification = c(1, 1),
                                 legend.title=element_text(size = 12),
                                 legend.text=element_text(size = 12))
g <- g + geom_hline(yintercept = 0, linetype = 1)
plot(g)
   
5.6.2 Lasso
#################
### Lasso回帰 ###
#################
soft_th <- function(lambda,x){
  return(sign(x) * pmax(abs(x) - lambda,0))
}

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.0001){
    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))
}

###
vc_names <- c("警察への年間資金","25歳以上で高卒人数の割合","16-19歳で高校に通っていない人の割合",
              "18-24歳で大学生の割合","25歳以上で4年制大卒の割合")

df <- read.table(file = "(ファイルパス)",header = F,sep = ",")

X <- df[,3:ncol(df)]
y <- df[,1]

p <- ncol(X)

lambda_seq <- seq(0, 100, 0.1)
coef_seq <- lambda_seq

df_betas <- data.frame(matrix(NA,nrow = p, ncol = length(lambda_seq)))
colnames(df_betas) <- 1:ncol(df_betas)

# Lasso推定
for(lambda in lambda_seq){
  df_betas[,which(lambda_seq==lambda)] <- lasso(X,y,lambda)$beta
}

df_betas_graph <- NULL

for(i in 1:nrow(df_betas)){
  df_betas_graph <- rbind(df_betas_graph,
                          data.frame(lambda = lambda_seq,
                                     var = vc_names[i],
                                     estimate = unlist(df_betas[i,]))
  )
}

# 図示
g <- ggplot(df_betas_graph,aes(x = lambda,y = estimate, group = var, color = var)) + geom_line()
g <- g + xlab("λ") + ylab("推定値")+ labs(title = "Lasso回帰による各母数の推定値",color = "")
g <- g + theme_classic() + theme(plot.title = element_text(hjust = 0.5),
                                 legend.position = c(1, 1),
                                 legend.justification = c(1, 1),
                                 legend.title=element_text(size = 12),
                                 legend.text=element_text(size = 12))
g <- g + geom_hline(yintercept = 0, linetype = 1)
plot(g)
   

5.7 Pythonによるシミュレーション

5.7.1 Ridge
import numpy as np
import copy
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib

def ridge(x,y, lam = 0):
    X = copy.copy(x)
    n, p = X.shape
    X_bar = np.zeros(p)
    s = np.zeros(p)
    for j in range(p):
        X_bar[j] = np.mean(X[:,j])
    for j in range(p):
        s[j] = np.std(X[:,j])
        X[:,j] = (X[:, j] - X_bar[j])/s[j]
    y_bar = np.mean(y)
    y = y - y_bar
    
    beta = np.linalg.inv(X.T@X + n * lam * np.eye(p))@X.T@y
    
    for j in range(p):
        beta[j] = beta[j] / s[j]
    beta_0 = y_bar - X_bar.T@beta
    return {'beta': beta, 'beta_0': beta_0}

# テストデータ
df = np.loadtxt("C:/Users/Julie/Downloads/crime.txt", delimiter = ",")
X = df[:,[i for i in range(2,7)]]
p = X.shape[1]
y = df[:,0]
lambda_seq = np.arange(0,50,0.5)
plt.xlim(0,50)
plt.ylim(-7.5,15)
plt.xlabel("lambda")
plt.ylabel("beta")
labels = ["警察への年間資金","25歳以上で高卒人数の割合","16-19歳で高校に通っていない人の割合",
          "18-24歳で大学生の割合","25歳以上で4年制大卒の割合"]

for j in range(p):
    coef_seq = []
    for l in lambda_seq:
        coef_seq.append(ridge(X,y,l)['beta'][j])
    plt.plot(lambda_seq, coef_seq, label = "{}".format(labels[j]))
plt.legend(loc = "upper right")
             
5.7.2 Lasso
#############
### Lasso ###
#############
import numpy as np
import copy

def soft_th(lam, x):
    return np.sign(x) * np.maximum(np.abs(x) - lam, 0)

def lasso(x,y,lam = 0):
    X = copy.copy(x)
    n, p = X.shape
    X_bar = np.zeros(p)
    s = np.zeros(p)
    
    for j in range(p):
        X_bar[j] = np.mean(X[:,j])
    for j in range(p):
        s[j] = np.std(X[:,j])
        X[:,j] = (X[:,j]-X_bar[j])/s[j]
    y_bar = np.mean(y)
    y = y - y_bar
    eps = 1
    
    beta = np.zeros(p)
    beta_old = np.zeros(p)
    
    while eps > 0.0001:
        for j in range(p):
            idx = list(set(range(p)) - {j})
            r = y - X[:,idx]@beta[idx]
            beta[j] = soft_th(lam,r.T@X[:,j]/n)
        eps = np.max(np.abs(beta - beta_old))
        beta_old = beta
    for j in range(p):
        beta[j] = beta[j]/s[j]
    beta_0 = y_bar - X_bar.T@beta
    return {'beta':beta, 'beta_0':beta_0}

# テストデータ
df = np.loadtxt("C:/Users/Julie/Downloads/crime.txt", delimiter = ",")
X = df[:,[i for i in range(2,7,1)]]
p = X.shape[1]
y = df[:,0]
lasso(X,y,20)

lambda_seq = np.arange(0,200,0.5)
plt.xlim(0,200)
plt.ylim(-10,20)
plt.xlabel("lambda")
plt.ylabel("beta")
labels = ["警察への年間資金","25歳以上で高卒人数の割合","16-19歳で高校に通っていない人の割合",
          "18-24歳で大学生の割合","25歳以上で4年制大卒の割合"]

for j in range(p):
    coef_seq = []
    for l in lambda_seq:
        coef_seq.append(lasso(X,y,l)['beta'][j])
    plt.plot(lambda_seq, coef_seq, label = "{}".format(labels[j]))
plt.legend(loc = "upper right")
plt.title("各λに対する各係数の推定値")

#####
from sklearn.linear_model import Lasso
from sklearn.linear_model import LassoCV

Las = Lasso(alpha = 20)
Las.fit(X,y)
Las.coef_

# alphasに指定した値をグリッドサーチする
Lcv = LassoCV(alphas = np.arange(0.1, 30, 0.1), cv = 10)
Lcv.fit(X,y)
Lcv.alpha_
Lcv.coef_
             

*1:\boldsymbol{x}=(x_1,\cdots,x_n)に対して、\mathcal{L}^1-ノルム\|\boldsymbol{x}\|および\mathcal{L}^2-ノルム\|\boldsymbol{x}\|_2\|\boldsymbol{x}\|=\sqrt{\displaystyle{\sum_{i=1}^{n}|x_i|}},\|\boldsymbol{x}\|_2=\sqrt{\displaystyle{\sum_{i=1}^{n}x_i^2}}で定義される。

*2:記号としては普通の微分と同じ書き方をする。

プライバシーポリシー お問い合わせ