はじめに
ロジスティック回帰は、入力変数$x$が実数全域の値を取り、出力変数$y$が$0<y<1$を満たすような場合において近似曲線を求めるために使われる。$0<y<1$という特徴から、0と1の判定つまり、『yes』『no』で表現される分類問題にも使用される。この例としては、電子メールのスパム判別が挙げられる。今回は、そのようなロジスティック回帰をPythonにおける最急降下法を用いて表現する。
ロジスティック曲線
$a,b$を実数とし、$x$を実数域を取る入力変数としたときの出力$y$が以下で与えられているような曲線をロジスティック曲線という。
y=\frac{1}{1+e^{ax+b}}
以下の回帰曲線から明らかなとおり、$0<y<1$の範囲内にどのような実数$x$でも必ず収まる。
損失関数
ロジスティック回帰曲線によって得られた入力変数$x$のときの近似出力$y_{opt}$と実際の出力変数である$y$との二乗偏差の和を損失関数と考え、その値が最小となるように最急降下法を用いて最適化する。
具体的には、サンプル点数が$n$であるとき、
Error=\frac{1}{n}\Sigma_{k=1}^{n}(y_k-y_{opt}(x_k))^2
ただし、
y_{opt}(x_k)=\frac{1}{1+e^{ax_k+b}}
とする。
このような二乗偏差の和である損失関数を最小化するときの$(a,b)$を最急降下法を用いて求める。
最急降下法
解空間の最小値を取る極値が一つ解空間の谷が一つの場合に有効な最小値の推定手法である。$(a,b)$に初期値を与え、以下のような式を用いて$(a,b)$を更新していく。つまり、$(a_k,b_k)$を算出する。
a_{k+1}=a_{k}-\alpha(\frac{\delta }{\delta a_{i}}Error(a_{k},b_{k}))
b_{k+1}=b_{k}-\beta(\frac{\delta }{\delta b_{k}}Error(a_k,b_k))
ただし$\alpha,\beta$は$(a,b)$の学習率である。
これは、以下のグラフのようにボールが坂に下り落ちる様に似ている。偏微分係数は各方向の坂の傾きと捉えることができる。
プログラム
上記のような考察を行うことで以下の様なプログラムを作成した。
近似曲線
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
num=10
x=np.linspace(-2,2,num)
y=[]
for i in range(num):
if x[i]<0:
y.append(0)
else:
y.append(1)
#plt.scatter(x,y)
h=1.0*10**-5
def error(a,b):
err=0
for i in range(num):
err=err+(1/num)*(y[i]-1/(1+np.exp((a*x[i]+b))))**2
return err
iterations=100
#初期値
a=1
b=1
alpha=2
beta=2
a_pre=[]
b_pre=[]
error_pre=[]
for i in range(iterations):
d_err_a=(error(a+h,b)-error(a,b))/h
d_err_b=(error(a,b+h)-error(a,b))/h
a=a-alpha*d_err_a
b=b-beta*d_err_b
a_pre.append(a)
b_pre.append(b)
error_pre.append(error(a,b))
plt.scatter(x,y,label="データ",color="blue")
plt.plot(x,1/(1+np.exp(a*x+b)),label="ロジスティック回帰予測",color="red")
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.savefig("ロジスティック回帰予測.png")
plt.show()
解の探索推移
一方で解の探索状態の推移は以下の様になプログラムを用いて調査した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
num=10
x=np.linspace(-2,2,num)
y=[]
for i in range(num):
if x[i]<0:
y.append(0)
else:
y.append(1)
#plt.scatter(x,y)
h=1.0*10**-5
def error(a,b):
err=0
y_opt=0
for i in range(num):
y_opt=1/(1+np.exp(a*x[i]+b))
err=err+(1/num)*(y[i]-y_opt)**2
return err
iterations=100
#初期値
a=1
b=1
alpha=2
beta=2
a_pre=[]
b_pre=[]
error_pre=[]
n=10
a_ary=np.linspace(-5,5,n)
b_ary=np.linspace(-5,5,n)
A,B=np.meshgrid(a_ary,b_ary)
Errors=np.zeros((n,n))
for i in range(n):
for k in range(n):
Errors[i][k]=error(A[i][k],B[i][k])
for i in range(iterations):
d_err_a=(error(a+h,b)-error(a,b))/h
d_err_b=(error(a,b+h)-error(a,b))/h
a=a-alpha*d_err_a
b=b-beta*d_err_b
a_pre.append(a)
b_pre.append(b)
error_pre.append(error(a,b))
plt.contourf(A,B,Errors,cmap="jet",alpha=0.5)
plt.scatter(a_pre,b_pre,c=error_pre,cmap="jet",label="最急降下法")
plt.legend()
plt.colorbar(label="誤差")
plt.xlabel("a")
plt.ylabel("b")
plt.savefig("ロジスティック回帰2乗誤差と最急降下法.png")
plt.show()
まとめ
今回は、線形単回帰に引き続き、非線形回帰の一種であるロジスティック回帰の実装を最急降下法を用いて行った。一般的には、対数関数で表される尤度関数を最大化させることでフィッテングを行うが、今回は線形単回帰と同様に、損失関数を最小化するという方針からフィッティングを行った。このような実装で気を付けるべき点は、線形単回帰とはことなりパラメータ$a,b$の重みが全く違うということである。したがって、そのようなパラメータの初期値に最適化した曲線の精度は大きく依存するため、適切な方法で初期パラメータを設定することが重要である。