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

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

MENU

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

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

を用いることにする。

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)
プライバシーポリシー お問い合わせ