ラプラス近似の説明をします。その簡単な応用として、ベイズ流のロジスティック回帰を定式化します。ついでにロジスティック回帰の説明もします。
最後に、pythonで実装してみます。PRMLの4章やとにかくすごいパンダ本の4章を参考にしてます。
記事で使っているソースコードはgithub に置いてあります。
https://github.com/msamunetogetoge
ラプラス近似
ラプラス近似は、結構荒っぽい近似ですが、計算が単純というメリットがあります。一方で、計算量がデータ数\(^2\)のオーダーで必要というデメリットがあります。
単純に言うと、任意の連続変数で、定義域が\((-\infty , \infty )\)である確率分布を正規分布で近似するのがラプラス近似です。
\(p(X) \)を確率分布とします。1対数尤度を最大化する\(X\)を\(X_{max} \)とします。つまり、\( \nabla \log p(X) |_{X=X_{max} }=0 \)を満たす点です。
\(X_{max} \)の周りで対数尤度をテーラー展開し、3次以降の項を無視した式を書き下します。
$$\begin{eqnarray}
\log p(X) &\sim &\log p(X_{max} ) – \frac{1}{2} (X-X_{max}) ^{T} \Lambda ^{-1} (X-X_{max}) \\
\Lambda &=& – \nabla \otimes \nabla \log p(X) |_{X=X_{max} }
\end{eqnarray}$$
右辺は、正規分布\(\mathcal{N}(X|X_{max} , \Lambda ^{-1} )\)の\( \log \)を取った式です。
つまり、対数尤度を最大にする点の周りでは、近似的に確率分布を正規分布だと思う事が出来ます。
\(p(X) \sim \mathcal{N}(X|X_{max} , \Lambda ^{-1} ) \)と近似する手法を、ラプラス近似(Laplace Approximation)と呼びます。
ロジスティック回帰
ロジスティック回帰は、分類問題に対する回帰分析みたいなものです。分類したい対象が2クラスの場合を考えます。\(t= (t_1 , \cdots , t_N ) . t_i \in \{ 0 ,1 \} \)がターゲットで、\(t=0 \)がクラス1, \(t=1 \)がクラス2に対応します。\(x \)をm次元の特徴量とすると、モデルは以下のように書けます。
$$\begin{eqnarray}
y &=& \sigma ( w^T x )\\
\sigma(a) &=& \frac{1}{1+e^{-a} }
\end{eqnarray}$$
\(y \)が0.5以上なら\(t=0 \), それ以外では\(t=1 \)というように予測値を決定します。適切な誤差関数を定めて、最適な\(w \)を決定するのがロジスティック回帰です。解説記事があるので良く分からない人は読んでみてください。
https://masamunetogetoge.com/classification-logisticregression
ベイズロジスティック回帰
ベイズ線形回帰のように、ロジスティック回帰もベイズ流で行ってみましょう。
ベイズ線形回帰については以下の記事を読んでみてください。
https://masamunetogetoge.com/byaes-linear-regression
ロジスティック回帰のように、パラメーター\(w\) と\(x\)を与えたときのラベルの確率分布は、シグモイド関数を使ったベルヌーイ分布です。
$$\begin{eqnarray}
p(t|w,x)= \prod_{n=1}^{N} \sigma (w^{T}x_n )^{t_n} (1-\sigma (w^{T}x_n ))^{1-t_n }
\end{eqnarray}$$
今は、2パラメーターを推定するのが目標なので、教師データのラベルと特徴量のセットを与えたときの事後分布\(p(w|t,x) \)を求める事が出来れば良いです。
ベイズの定理より、
$$\begin{eqnarray}
p(w|t,x ) \propto p(t|w,x) p(w|x)
\end{eqnarray}$$
が成り立ちます。事前分布\(p(w|x ) \)は適当に仮定してやる必要があります。ラプラス近似で事後分布を求めたいので、正規分布を仮定しておきます。
$$\begin{eqnarray}
p(w|x) = \mathcal{N}(w| w_0 , S_0 )
\end{eqnarray}$$
ラプラス近似で事後分布を求めるため、対数尤度を書き下しましょう。
$$\begin{eqnarray}
\log p(w|t) &=& -\frac{1}{2}(w-w_0)^{T} S_0 ^{-1} (w-w_0) \\
&+&\sum \left( t_n \log \sigma(w^T x_n ) +(1-t_n)(\sigma (w^T x_n ) )\right)
\end{eqnarray}$$
シグモイド関数の微分は、以下のようになります。
$$\begin{eqnarray}
\sigma^{,} (x)=\sigma (x) (1-\sigma (x) )
\end{eqnarray}$$
これを使って、
$$\begin{eqnarray}
\Lambda (w) &=&-\Delta _w \log p(w|t) \\
&=&-S_0 ^{-1} +\sum \sigma(w^T x_n ) (1- \sigma(w^T x_n ) ) x_n \cdot x_n ^T +{\rm const}
\end{eqnarray}$$
と計算出来ます。3 この結果と、対数尤度を最大にする \(w=w_{max} \)を使って、
$$\begin{eqnarray}
p(w|t) \sim \mathcal{N}(w|w_{max} , \Lambda(w_{max} ))
\end{eqnarray}$$
となります。この分布からパラメーター\(w\)をサンプリングして、
$$\begin{eqnarray}
p(t|w,x)= \prod_{n=1}^{N} \sigma (w^{T}x_n )^{t_n} (1-\sigma (w^{T}x_n ))^{1-t_n }
\end{eqnarray}$$
を計算することで2値分類を行います。
irisデータセットを使ってベイズロジスティック回帰を試してみます。
python による実装
初めに対数尤度の計算と、確率勾配法を実装しておきます。
import numpy as np
from scipy.special import expit as sigmoid
def likelifood(w,w_0,S,t,x):
w_1 = w-w_0
y=sigmoid(np.dot(w.T , x)).reshape(-1,)
lf = -1/2 * np.dot(np.dot(w_1.T , np.linalg.inv(S)) ,w_1) +np.sum(t*np.log(y) +(1-t)*np.log(1-y))
return lf
class SGD():
def __init__(self,w_0, S_0,x,t,rate,th):
self.r=rate
self.w_0 =w_0
self.S_0=S_0
self.x=x
self.t=t
self.th =th
def diff(self,w):
Nd = -np.dot((w-self.w_0).T,np.linalg.inv(self.S_0.T))
Br=np.dot(t.T-sigmoid(np.dot(w.T,self.x)),self.x.T)
d=Nd+Br
d=d.reshape(-1,1)
return d
def learn(self,N):
self.w=self.w_0
for i in range(N):
self.d=self.diff(self.w)
self.w_new =self.w-self.r*self.d
lf_new = likelifood(self.w_new,self.w_0,self.S_0,self.t,self.x)
lf_old = likelifood(self.w,self.w_0,self.S_0,self.t,self.x)
c = (lf_new - lf_old)/lf_old
self.w =self.w_new
self.lf=lf_new
print("End iteration of {}".format(i+1))
if np.abs(c)<self.th:
print("Convergence!! itr={}".format(i+1))
break
print("End the learning")
return self.w_newirisデータセットを準備して、2値分類の問題にします。
import pandas as pd
from sklearn import datasets
iris = datasets.load_iris()
df = pd.DataFrame(iris.data, columns=iris.feature_names)
df['target'] = iris.target_names[iris.target]
data =df.target=="setosa"
df_reg_index = df[data].index
df_reg=df.drop(index=df_reg_index)
y=df_reg.target
x=df_reg.drop(columns="target")一応特徴量は正規化してからベイズロジスティック回帰しましょう。
from sklearn.preprocessing import OrdinalEncoder
oe = OrdinalEncoder()
encoded = oe.fit_transform(np.array(y).reshape(-1,1))
t=encoded
stats = x.describe().T
def norm(x):
return (x - stats['mean']) /stats['std']
m=len(x.columns)
w_0 =np.zeros((m,1))
S=np.eye(m)
x=norm(x).T
SG=SGD(w_0=w_0, S_0=S, x=x, t=t, rate=0.0001,th=0.00001)
w_new=SG.learn(10000)事後分布の計算には\(w_{new} \)だけでなく、分散も必要でした。
$$\begin{eqnarray}
\Lambda (w) &=&-S_0 ^{-1} +\sum \sigma(w^T x_n ) (1- \sigma(w^T x_n ) ) +{\rm const}
\end{eqnarray}$$
y=sigmoid(np.dot(w_new.T , x))
S=-np.linalg.inv(S) + np.sum(y*(1-y)*np.dot(x.T, x))
\(p(w|t) \sim \mathcal{N}(w|w_{new},\Lambda (w_{new} )) \) であったので、適当にパラメーターをサンプルして混合行列を書いてみます。
from sklearn.metrics import confusion_matrix
w_s =np.random.multivariate_normal(w_new.reshape(-1,), S,5)
for i in w_s :
y= sigmoid(np.dot(i , x))
pred=[]
for i in range(max(y.shape)):
pred.append(np.argmax([y[i], 1-y[i]]))
cm = confusion_matrix(t, pred)
print(cm)
###
[[44 6]
[ 9 41]]
[[44 6]
[ 9 41]]
[[44 6]
[ 6 44]]
[[44 6]
[ 9 41]]
[[44 6]
[ 9 41]]
###大体80%~90%の正答率です。sklearn のロジスティック回帰で行うと、95%くらいの正答率が出ます。4
まとめ
- ラプラス近似の説明をした
- ロジスティック回帰の説明をした
- ベイズロジスティック回帰の説明をした
- pythonで計算した
