統計的機械学習の数理100問（08/20） - 「大人の教養・知識・気付き」を伸ばすブログ

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

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

を用いることにする。

前回
2.　分類
- 2.1.　ロジスティック回帰
- 2.2　Newton-Raphson法の適用
  - 2.2.1　2変数関数での適用例(R)
  - 2.2.2　2変数関数での適用例(Python)
次回

前回

power-of-awareness.com

2.　分類

　有限個の数値を目的変数とするような分類の問題を取り扱う。

2.1.　ロジスティック回帰

　目的変数が2個のいずれかの値を取り、 $p\geq1$ 個の説明変数からそのいずれであるのかを対応付ける問題を考える。すなわちデータ

$\begin{aligned} (\boldsymbol{X}_1,Y_1),\cdots,(\boldsymbol{X}_n,Y_n)\in\mathbb{R}^{p}\times\{-1,1\} \end{aligned}$

から決定則

$\begin{aligned} f:&\mathbb{R}^{p}\rightarrow\{-1,1\}\\ f:&\boldsymbol{X}\mapsto Y \end{aligned}$

を導く。
　ここでは

$\begin{aligned} P(Y=1)&=\displaystyle{\frac{{}^t\boldsymbol{X}\boldsymbol{\beta}}{1+e^{\beta_0+{}^t\boldsymbol{X}\boldsymbol{\beta}_1}}},\\ P(Y=-1)&=\displaystyle{\frac{1}{1+e^{\beta_0+{}^t\boldsymbol{X}\boldsymbol{\beta}_1}}} \end{aligned}$

となるような $\boldsymbol{\beta}=(\beta_0,\boldsymbol{\beta}_1)\in\mathbb{R}^{p+1}$ が存在すると仮定する。このような形態の回帰モデルをロジスティック回帰モデルという。

2.1.1　ロジスティック曲線のシミュレーション(R)

# ロジスティック関数
func_logis <- function(x)exp(num_beta_0 + num_beta * x)/(1 + exp(num_beta_0 + num_beta * x))

num_beta_0 <- 0
vc_beta <- c(0, 0.2, 0.5, 0.5, 1, 2, 10)

num_m <- length(vc_beta)
num_beta <- vc_beta[1]

plot(func_logis, xlim = c(-10,10), ylim = c(0,1), xlab = "x", ylab = "P(Y=1|x)", col = 1,
     main = "ロジスティック曲線")

for(i in 2:num_m){
  num_beta <- vc_beta[i]
  par(new = T)
  plot(func_logis, xlim = c(-10,10), ylim = c(0,1), xlab = "", ylab = "", axes = F, col = i)
}

legend("topleft", legend = vc_beta, col = 1:length(vc_beta), lwd = 2, cex = 0.8)

2.1.2　ロジスティック曲線のシミュレーション(Python)

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib

def logis(x):
    return np.exp(beta_0 + beta * x) / (1+np.exp(beta_0 + beta * x))

beta_0 = 0

beta_seq = np.array([0, 0.2, 0.5, 1, 2, 10])
x_seq = np.arange(-10,10,0.1)

# 図示
plt.xlabel("x")
plt.ylabel("P(Y=1|x)")
plt.title("ロジスティック曲線")

for i in range(beta_seq.shape[0]):
    beta = beta_seq[i]
    p = logis(x_seq)
    plt.plot(x_seq, p, label = '{}'.format(beta))
plt.legend(loc = 'upper left')

2.1.3　パラメータの推定

　パラメータの推定は最尤法で行うこととする。すなわち各 $(\boldsymbol{X}_i,Y_i)$ が独立かつ同一分布に従うとし、尤度 $L(\boldsymbol{\beta})$

$\begin{aligned} L(\boldsymbol{\beta})=\displaystyle{\prod_{i=1}^{n}\frac{1}{1+\exp\left\{-Y_i(\beta_0+{}^{t}\boldsymbol{X}_i\boldsymbol{\beta}_1)\right\}}}, \end{aligned}$

を最大にするような、さらにはそれと同値だが対数関数が単調増加であることを活用し対数尤度 $l(\boldsymbol{\beta})$

$\begin{aligned} l(\boldsymbol{\beta})=\log L(\boldsymbol{\beta})=\displaystyle{\sum_{i=1}^{n}\log\left\{1+e^{-Y_i(\beta_0+{}^{t}\boldsymbol{X}_i\boldsymbol{\beta}_1)}\right\}} \end{aligned}$

を最大にするような $\boldsymbol{\beta}\in\mathbb{R}^{p+1}$ を求める。ただしロジスティック回帰の最尤法は代数的に解けるわけではないため、数値的に解くことを考える。

2.2　Newton-Raphson法の適用

　 $C^{l},l\in\mathbb{N}$ 級の関数 $f(x)$ に対して $f(x)=0$ を解くことを考える。適当な初期値 $x=x_0$ が与えられているものとして、点 $(x_0,f(x_0))$ から $f(x)$ の接線を引いて $x$ 軸との交点を $x=x_1$ とし、次に点 $(x_1,f(x_1))$ において同様のことを行うという手順を繰り返す。こうして得られた数列 $\{x_i\}$ は $f(x)=0$ の解に近づいていく。
　点 $(x_k,f(x_k))$ における $f(x)$ の接線の方程式は一般に

$\begin{aligned} y-f(x_k)=f^{\prime}(x_k)(x-x_k) \end{aligned}$

で与えられる。したがって $y=0$ との交点 $x_{k+1}$ は

$\begin{aligned} x_{k+1}=x_k-\displaystyle{\frac{f(x_k)}{f^{\prime}(x_k)}} \end{aligned}$

で得られる。試行回数をあらかじめ決めておく、もしくは $|x_{k+1}-x_k|$ の上限値を決めてそれよりも絶対誤差が小さくなったときに終了させるなど、予め終了条件を決めておくことで数値解を得ることができる。
　一般の次元 $p$ に関して $\nabla l(\boldsymbol{\beta})=\boldsymbol{0}$ を満たすような $\boldsymbol{\beta}={}^{t}(\beta_0,\boldsymbol{\beta}_1)\in\mathbb{R}^{p+1}$ を求める問題に $\mathrm{Newton}$ - $\mathrm{Raphson}$ 法を適用する。すなわち

$\begin{aligned} (\beta_0,\boldsymbol{\beta}_1)\leftarrow(\beta_0,\boldsymbol{\beta}_1) -\{\nabla^{2}l(\beta_0,\boldsymbol{\beta}_1)\}^{-1}\nabla l(\beta_0,\boldsymbol{\beta}_1) \end{aligned}$

により逐次的に $(\beta_0,\boldsymbol{\beta}_1)$ を更新していく。ただし

$\begin{aligned} \nabla f(\boldsymbol{v})&=\begin{bmatrix}\displaystyle{\frac{\partial f}{\partial v_1}}\\ \vdots\\ \displaystyle{\frac{\partial f}{\partial v_{p+1}}}\end{bmatrix}\in\mathbb{R}^{p+1},\\ \nabla^2 f(\boldsymbol{v})&=\begin{bmatrix} \displaystyle{\frac{\partial^2 f}{\partial {v_1}^2}}&\displaystyle{\frac{\partial^2 f}{\partial v_1\partial v_2}}&\cdots&\displaystyle{\frac{\partial^2 f}{\partial v_1\partial v_{p+1}}}\\ \displaystyle{\frac{\partial^2 f}{\partial v_2\partial v_1}}&\displaystyle{\frac{\partial^2 f}{\partial {v_2}^2}}&\cdots&\displaystyle{\frac{\partial^2 f}{\partial v_2\partial v_{p+1}}}\\ \vdots&&\ddots&\vdots\\ \displaystyle{\frac{\partial^2 f}{\partial v_{p+1}\partial v_1}}&\displaystyle{\frac{\partial^2 f}{\partial v_{p+1}\partial v_2}}&\cdots&\displaystyle{\frac{\partial^2 f}{\partial v_{p+1}\partial v_{p+1}}} \end{bmatrix}\in\mathbb{R}^{(p+1)\times(p+1)},\\ \end{aligned}$

とする。
　ロジスティック回帰にこれを適用する。負の対数尤度 $-l(\beta_0,\boldsymbol{\beta}_1)$ を微分すると $v_i=e^{Y_i(\beta_0+\boldsymbol{X}_i\boldsymbol{\beta}_1)},i=1,2,\cdots,n$ として

$\begin{aligned} \nabla L&=\begin{bmatrix} \displaystyle{\frac{\partial L}{\partial\beta_1}}\\ \vdots\\ \displaystyle{\frac{\partial L}{\partial\beta_n}} \end{bmatrix} =-{}^{t}\boldsymbol{X}\boldsymbol{u}\in\mathbb{R}^{p+1},\\ \boldsymbol{u}&=\begin{bmatrix} \displaystyle{\frac{Y_1 v_1}{1+v_1}}\\ \vdots\\ \displaystyle{\frac{Y_n v_n}{1+v_n}} \end{bmatrix} \end{aligned}$

と書ける。さらに $y_i^2=1$ であることに注意すれば、

$\begin{aligned} \nabla^2 L&=\begin{bmatrix} \displaystyle{\frac{\partial^2 L}{\partial\beta_1^2}}&\displaystyle{\frac{\partial^2 L}{\partial\beta_1\partial\beta_2}}&\cdots&\displaystyle{\frac{\partial^2 L}{\partial\beta_1\partial\beta_{p+1}}}\\ \displaystyle{\frac{\partial^2 L}{\partial\beta_2\partial\beta_1}}&\displaystyle{\frac{\partial^2 L}{\partial\beta_2^2}}&\cdots&\displaystyle{\frac{\partial^2 L}{\partial\beta_2\partial\beta_{p+1}}}\\ \vdots&\ddots&&\vdots\\ \displaystyle{\frac{\partial^2 L}{\partial\beta_{p+1}\partial\beta_1}}&\displaystyle{\frac{\partial^2 L}{\partial\beta_{p+1}\partial\beta_2}}&\cdots&\displaystyle{\frac{\partial^2 L}{\partial\beta_{p+1}\partial\beta_{p+1}}} \end{bmatrix}={}^{t}\boldsymbol{X}W\boldsymbol{X}\in\mathbb{R}^{(p+1)\times(p+1)},\\ W&=\begin{bmatrix} \displaystyle{\frac{v_1}{(1+v_1)^2}}&\cdots&0\\ \vdots&\ddots&\vdots\\ 0&\cdots&\displaystyle{\frac{v_n}{(1+v_n)^2}} \end{bmatrix} \end{aligned}$

である。そのような $W,\boldsymbol{u}$ を用いると、

$\begin{aligned} \boldsymbol{\beta}\leftarrow\boldsymbol{\beta}+\left({}^{t}\boldsymbol{X}W\boldsymbol{X}\right)^{-1}{}^{t}\boldsymbol{X}\boldsymbol{u} \end{aligned}$

と書ける。また $z=\boldsymbol{X}\boldsymbol{\beta}+W^{-1}\boldsymbol{u}\in\mathbb{R}^{n}$ とおけば

$\begin{aligned} \boldsymbol{\beta}\leftarrow({}^{t}\boldsymbol{X}W\boldsymbol{X}{}^{t}\boldsymbol{X}W\boldsymbol{z}) \end{aligned}$

とも書き表すことができる。

2.2.1　2変数関数での適用例(R)

　関数

$\begin{aligned} \begin{eqnarray}\left\{ \begin{array}{l} f(x,y)=x^2+y^2-1=0,\\ g(x,y)=x+y=0 \end{array} \right. \end{eqnarray} \end{aligned}$

を満たす解 $(x,y)$ を初期値として $(x_0,y_0)=(3,4)$ を与えて $\mathrm{Newton}$ - $\mathrm{Raphson}$ 法で解く。

##############################
### Newton-Raphson法の検討 ###
##############################

# f(x,y) = x^2 + y^2 -1, g(x,y) = x+yについて初期値(x,y) = (3,4)から
# f(x,y)=0かつg(x,y)=0を満たす(x,y)を数値的に解く

# 解く対象となる関数およびその1次偏微分をそれぞれ定義する
func_f <- function(vc) vc[1]^2 + vc[2]^2 - 1
func_f_pd_x <- function(vc) 2 * vc[1]
func_f_pd_y <- function(vc) 2 * vc[2]

func_g <- function(vc) vc[1] + vc[2]
func_g_pd_x <- function(vc) 1
func_g_pd_y <- function(vc) 1

# 初期値の設定
vc <- c(3,4)

eps <- 10e-14 # 絶対誤差がeps以下になったら終了

err <- Inf # 絶対誤差
num_loop <- 0 # ループ回数

while(err>eps){
  vc_bf <- vc
  vc <- vc_bf - solve(matrix(c(func_f_pd_x(vc_bf),func_f_pd_y(vc_bf),
                               func_g_pd_x(vc_bf),func_g_pd_y(vc_bf)),
                          ncol = 2, byrow = T)) %*% c(func_f(vc_bf),func_g(vc_bf))
  num_loop <- num_loop + 1
  err <- sqrt(t(vc_bf - vc) %*% (vc_bf - vc))
}

# 出力
vc
num_loop

2.2.2　2変数関数での適用例(Python)

import numpy as np

# 関数f(x,y),g(x,y)およびその1次偏微分を定義
def f(z):
    return z[0]**2 + z[1]**2-1
def f_pd_x(z):
    return 2*z[0]
def f_pd_y(z):
    return 2*z[1]
def g(z):
    return z[0] + z[1]
def g_pd_x(z):
    return 1
def g_pd_y(z):
    return 1

# Newton-Raphson法を適用する
z = np.array([3,4])

eps = 10**(-14)
err = float('inf')
num_loop = 0

while(err>eps):
    z_bf = z
    J = np.array([[f_pd_x(z),f_pd_y(z)],[g_pd_x(z),g_pd_y(z)]]) # Jacobian
    z = z_bf - np.linalg.inv(J)@np.array([f(z),g(z)]) 
    err = np.sqrt((z_bf - z).T@(z_bf - z))
    num_loop += 1

# 出力
print(z)
print(num_loop)

次回

power-of-awareness.com

前回

2. 分類

2.1. ロジスティック回帰

2.1.1 ロジスティック曲線のシミュレーション(R)

2.1.2 ロジスティック曲線のシミュレーション(Python)

2.1.3 パラメータの推定

2.2 Newton-Raphson法の適用

2.2.1 2変数関数での適用例(R)

2.2.2 2変数関数での適用例(Python)

次回