第５章 正則化¶

51¶

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy import stats
from numpy.random import randn
import copy
%matplotlib inline
# import japanize_matplotlib　# Google Colabの場合
In [ ]:
# Anacondaの場合は下記
import matplotlib
from matplotlib import font_manager
matplotlib.rc("font", family="BIZ UDGothic")
In [ ]:
def soft_th(lam, x):
    return np.sign(x) * np.maximum(np.abs(x)-lam, 0)
In [ ]:
x_seq = np.arange(-10, 10, 0.1)
plt.plot(x_seq, soft_th(5, x_seq))
plt.plot([-5, -5], [4, -4], c="black", linestyle="dashed", linewidth=0.8)
plt.plot([5, 5], [4, -4], c="black", linestyle="dashed", linewidth=0.8)
plt.title("soft_th(lam,x)")
plt.text(-1.5, 1, "λ=5", fontsize=15)

53¶

In [ ]:
def lasso(x, y, lam=0):  # lamはlambdaの略
    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.001:
        for j in range(p):
            index = list(set(range(p))-{j})
            r = y - X[:, index]@beta[index]
            beta[j] = soft_th(lam, r.T@X[:, j] / n)
        eps = np.max(np.abs(beta - beta_old))
        beta_old = copy.copy(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}
In [ ]:
df = np.loadtxt("crime.txt", delimiter="\t")
X = df[:, [i for i in range(2, 7)]]
p = X.shape[1]
y = df[:, 0]
In [ ]:
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歳以上で４年制大学を卒業した人の割合"]
for j in range(p):
    coef_seq = []
    for lam in lambda_seq:
        coef_seq.append(lasso(X, y, lam)["beta"][j])
    plt.plot(lambda_seq, coef_seq, label="{}".format(labels[j]))
plt.legend(loc="upper right")
plt.title("各λについての各係数の値")

54¶

In [ ]:
def ridge(x, y, lam=0):  # lamはlambdaの略
    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}
In [ ]:
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歳以上で４年制大学を卒業した人の割合"]
for j in range(p):
    coef_seq = []
    for lam in lambda_seq:
        coef_seq.append(ridge(X, y, lam)["beta"][j])
    plt.plot(lambda_seq, coef_seq, label="{}".format(labels[j]))
plt.legend(loc="upper right")
plt.title("各λについての各係数の値")

55¶

In [ ]:
from sklearn.linear_model import Lasso
from sklearn.linear_model import LassoCV
In [ ]:
Las = Lasso(alpha=20)
Las.fit(X, y)
Las.coef_
In [ ]:
# alphasに指定した値をグリッドサーチする, cvはKの数
Lcv = LassoCV(alphas=np.arange(0.1, 30, 0.1), cv=10)
Lcv.fit(X, y)
Lcv.alpha_
Lcv.coef_