目次

1. はじめに

   • ロジスティック回帰とは

     • ロジスティック回帰とは

       ベルヌーイ分布(0, 1)に従う目的変数の，クラスに属する確率を回帰する機械学習モデル．

       1

       ロジスティック回帰は二値分類モデルです．をはじめに持ってくる

       主に二値分類で使われる．

       (クラスに属する確率が50%以下→クラス0，50%以上→クラス1って感じで二値に分類する．)

       今回は二値分類のロジスティック回帰を実装する．

       <aside> 💡

       クラスに属する確率とは ロジスティック回帰で二値分類するときに，実際は，まずクラスに属する確率をまず求める．(~%の確率で1である．みたいに．)

       その後，確率50%を閾値に0, 1どちらに分類するのかを決める．

       クラスに属する確率は，二値の片方についてしか考えない．

       なぜなら一方のクラスに属する確率が決まれば，もう片方の確率は決まるから．(確率の総和は100%であることから)

       これが，クラスに属する確率である．

       </aside>

     • 前提知識 ロジスティック回帰を学ぶにあたり，前提知識となる線形回帰について，復習する．

       • 線形回帰(単回帰，重回帰)

       1つ，または複数の説明変数から，1つの連続値を予測するモデル．

       $$ \begin{align*} &n:データ数\\ &m:次元数\\ &X\in \mathbb{R}^{n\times m}:説明変数\\ &\boldsymbol{y}\in \mathbb{R}^{n}:目的変数\\ &\boldsymbol{w}\in \mathbb{R}^{m}:パラメータ\\ &b:バイアス項 \end{align*} $$

       線形回帰はこの式で表される．

       $$ \boldsymbol{\hat{y}}=X\boldsymbol{w} $$

       $\boldsymbol{\hat{y}}=X\boldsymbol{w}$は何者だろうか．

       それにあたり，まず1つの説明変数の線形回帰，単回帰について考える．

       それは，ただの一次関数．

       最初の単回帰から説明変数をベクトルにする．

       $$ y=wx+b $$

       では複数の説明変数の線形回帰，重回帰だとどうなるのか．

       それは，単回帰をひたすら足すようになる．

       wがベクトルになってる．

       $$ \boldsymbol{\hat{y}}=\boldsymbol{w_0}\boldsymbol{x_0}+\boldsymbol{w_1}\boldsymbol{x_1}+\boldsymbol{w_2}\boldsymbol{x_2}+\cdots+\boldsymbol{w_m}\boldsymbol{x_m}\\ ※\boldsymbol{b}=\boldsymbol{w_0}\boldsymbol{x_0}=1\boldsymbol{x_0}=\boldsymbol{x_0} $$

       これを，行列を使って整理すると以下のようになる．

       $$ \boldsymbol{\hat{y}}=\sum_{i=1}^{n}({\boldsymbol{w}^{T}\boldsymbol{x_i}})=X\boldsymbol{w} $$

       この資料では， $X\boldsymbol{w}$よりも $\sum_{i=1}^{n}({\boldsymbol{w}^{T}\boldsymbol{x_i}})$の方が多く登場する．

       (詳しくは線形回帰の実装をみて)

       線形回帰の回帰直線は以下のようなグラフになる．

       

       ひとまず，線形回帰の式は ，$\boldsymbol{\hat{y}}=X\boldsymbol{w}$で，図のような直線が引かれるということを理解してほしい．

     • 線形回帰との違い

       さっき述べたように，ロジスティック回帰は，クラスに属する確率を回帰する．

       確率を回帰するという点が線形回帰との大きな違いである．

       この点はミソなので，よく覚えておいてほしい．

       (逆に言えば，線形回帰は，確率を回帰していない．)

2. 気持ち 線形回帰から，ロジスティック回帰への一連の流れを，数式なしで図で理解する．

   • 線形回帰 まず，線形回帰では以下のような図が作られる．

     データは連続値で，それに対してモデルがフィットした図．

     

     二値分類 ロジスティック回帰は二値分類であるので，データが二値のケースを考える．

   まず前提として，線形回帰とは取り扱うデータがことなるのである．

   縦軸を見れば，二値分類のときは，0~1になっていることが伺える． (適当に，縦軸の値が0以下をとっていればクラス0に，0以上をとっていればクラス1に，データを加工した．)

   

   では，この図にいい感じの直線を引けるだろうか？

   →否．

   ということで，この図にいい感じの曲線を引いてみる．

   <aside> 💡

   いい感じの直線とは，データをうまくなぞるような，残差が最小である直線のことを指す．

   二値分類のこのようなデータは直線で表現できるとは言い難い．

   (実際に直線を引こうとしてみたらわかる．)

   </aside>

   • シグモイド関数の適用

   シグモイド関数に線形回帰を適用すると，図のような曲線に変換される．

   <aside> 💡

   ロジスティック回帰でのこのような変換は，ロジスティック関数(シグモイド関数)を利用しているため，ロジット変換と呼ばれたりする．

   機械学習の分野では，ロジスティック関数とシグモイド関数はほぼ同じものとして扱うことが多い．

   厳密にはロジスティック関数は．シグモイド関数の特定の形式を指す．

   </aside>

   直線よりかはうまく線を引けているはず．

   

   しかし，まだデータと曲線がかけ離れている．これは，パラメータを最適化することで解決する．

   • 最適なパラメータを適用する．

   最適化されたパラメータを適用する．

   どうやって最適化するのかは後に述べるが，パラメータを最適化することで，いい感じの曲線を引くことができる．

   

   このようなモデルをロジスティック回帰と言う．

   ロジスティック回帰は，二値のデータをうまく表現できている．

3. 理論 気持ちでやったことを，数式で表す．実装するには数式の理解が必須だから． 資料を二つ開いて，気持ちと理論を見比べるとわかりやすいかも．

   • 線形回帰

   気持ちでは，資料の書き方的に，データが既にあって，それに対し最適化された線形回帰を用意していたようにみえるが，それは誤解なので注意．

   まず，線形回帰を用意する．これがロジスティック回帰の基盤だからね．

   (tに意味はない．適当に置いた，代入するための変数．)

   $$ \begin{align*} &n:データ数\\ &m:次元数\\ &X\in \mathbb{R}^{n\times m}:説明変数\\ &\boldsymbol{y}\in \mathbb{R}^{n}:目的変数\\ &\boldsymbol{w}\in \mathbb{R}^{m}:パラメータ\\ \end{align*}

   $$

   $$ \boldsymbol{t}=\sum_{i=1}^{n}({\boldsymbol{w}^{T}\boldsymbol{x_i}})=X\boldsymbol{w} $$

   • シグモイド関数

   まずシグモイド関数を用意する．

   $$ \sigma(\boldsymbol{t}) = \frac{1}{1+e^{-\boldsymbol{t}}} $$

   このシグモイド関数に線形回帰を適用する．

   $$ \boldsymbol{\hat{y}} =\sigma(X\boldsymbol{w})=\sum_{i=1}^{n}\sigma({\boldsymbol{w}^{T}\boldsymbol{x_i}})= \frac{1}{1-e^{-X\boldsymbol{w}}} $$

   • パラメータ推定

   次に，最適なパラメータを求める．

   パラメータを最適化し，データに対して当てはまりの良い曲線にする．

   パラメータを最適化するには，パラメータについての具体的な式がわからなければいけない．(???を求めるということ．)

   $$ \boldsymbol{w}= ??? $$

   そのために，最尤法を利用し損失関数を導出する．

   その後，最適化アルゴリズムで損失が最小のパラメータを求める．

   困惑した人は，「パラメータを最適化することで，データに対し当てはまりの良い曲線にする．」だけ理解していればよし．

4. 実装

class LogisticRegression():
    def __init__(self, tol: float = 0.0001, max_iter: int = 100):
        self.tol = tol # 近似の許容誤差．
        self.max_iter = max_iter # 最大で何回更新するか．

    @staticmethod
    def _sigmoid(Xw):
        return ...

    def fit(self, x, y):
        self.w = ...
        self.r = ...
        tol_vec = ...
        tol_ = ...
        while ...:

    def predict(self, x):
        ...
        return y_pred

class LogisticRegression():
    def __init__(self, tol: float=0.0001, max_iter: int=100):
        self.tol = tol
        self.max_iter = max_iter

		# シグモイド関数を実装する．
    @staticmethod
    def _sigmoid(...):
        return ...

		# シグモイド関数に線形回帰を適用する．
    def predict_proba(self, x):
        return ...

    def _get_w(self, x, y, y_hat):
        for _ in range(X.shape[0]):
            self.R[_,_] = y_hat[_] * (1 - y_hat[_])
        xr = x.T @ self.R
        w_new = (np.linalg.pinv(xr @ x) @ xr) @ (x @ self.w - np.linalg.pinv(self.R) @ (y_hat - y))
        diff = w_new - self.w
        return w_new, diff

    def fit(self, x, y):
        self.w = np.random.randn(x.shape[1])
        self.R = np.identity(x.shape[0])
        tol_vec = np.full(x.shape[1],self.tol)
        for i in range(self.max_iter):
            y_hat = self._predict(x)
            w_new, diff = self._get_w(x, y, y_hat)
            if np.all(diff < tol_vec):
                break
            elif np.any(diff > tol_vec):
                self.w = w_new

    def predict(self, x):
        # 一方のクラスに属する確率を算出．学習済みモデルから．取り出せるようにしてほしい．
        # 説明変数のデータ数繰り返し．
        # 0.5を閾値にクラスを分ける．
        return y_pred
        
        def predict(self, x):
        self.predict_proba = _sigmoid(np.dot(x, self.w))
        y_pred = self.predict_proba.copy()
        for _ in range(x.shape[0]):
            if y_pred[_] > 0.5:
                y_pred[_] = 1
            elif y_pred[_] < 0.5:
                y_pred[_] = 0
        return y_pred

class LogisticRegression():
    def __init__(self, tol: float=0.0001, max_iter: int=100):
        self.tol = tol # パラメータの許容誤差的なもの．更新後のパラメータと，更新前のパラメータの差がこれより小さくなればOKということ．
        self.max_iter = max_iter # 最大の更新回数．

# シグモイド関数を実装する．
    @staticmethod
    def _sigmoid(...): # メソッドとして使用することはないので先頭に_をつける
        return ...

    def fit(self, x, y):
        # パラメータの初期化．ランダムでOK
        # 対角行列の初期化．データ数のshapeの単位行列作成して初期化すると良い．
        # パラメータと同じshapeのベクトルで，中の要素がtolで初期化する．これを基準にbreakするか決める．     
        # パラメータと同じshapeのベクトルで，新パラメータと旧パラメータの差を計測するベクトルを作る．絶対に0より大きくならないといけないのでnp.infとかで初期化すると良い．
        # イテレーション数を初期化．この変数に更新数を+1して，更新回数を計測する．
        # while文で，diffのどれか一つでもtol_vecより大きかったらtrue．かつ今の更新回数がmax回数より小さかったらtrueにする．
            # 一方のクラスに属する確率を算出する．いわゆるy＾を作る．
            # データ数だけfor文で繰り返し．
                # 単位行列に対角行列の各要素を代入していく．
            # XRを算出
            # 新たなパラメータを求める．
            # 新パラメータと旧パラメータの差をとり，diffに代入．
            # 更新数+1
            # パラメータを更新する．新パラメータを旧パラメータとするということ．

    def predict(self, x):
        # 一方のクラスに属する確率をpredict_probaとして算出する．後で取り出せるようにしといて．
        # ここでpredict_probaをpredとしてコピーしておく．
        # データ数だけ繰り返す．
            # predが0.5より大きいなら，そのpredはクラス1とする．
            # でなければクラス0
        return y_pred

class LogisticRegression():
    def __init__(self, tol: float=0.0001, max_iter: int=100):
        self.tol = tol # 停止基準の許容範囲
        self.max_iter = max_iter # 最大の更新回数．

    @staticmethod
    def _sigmoid(...): # メソッドとして使用することはないので先頭に_をつける
		    # シグモイド関数を実装する．
        return ...

    def fit(self, x, y):
        # 学習部分

    def predict(self, x):
        # 推論部分
        return ...

答え

class SLogisticRegression():
    def __init__(self, tol: float=0.0001, max_iter: int=100):
        self.tol = tol
        self.max_iter = max_iter

    @staticmethod
    def _sigmoid(Xw):
        return 1 / (1 + np.exp(-Xw))

    def fit(self, x, y):
        self.w = np.random.randn(x.shape[1])
        self.R = np.identity(x.shape[0])
        tol_vec = np.full(x.shape[1],self.tol)
        diff = np.full(x.shape[1], np.inf)
        tol_=0
        while np.any(diff > tol_vec) and (tol_ < self.max_iter):
            y_hat = self._sigmoid(x @ self.w)
            for _ in range(X.shape[0]):
                self.R[_,_] = y_hat[_] * (1 - y_hat[_])
            xr = x.T @ self.R
            w_new = (np.linalg.solve(xr @ x, xr)) @ (x @ self.w - np.linalg.solve(self.R, (y_hat - y)))
            diff = w_new - self.w
            tol_ += 1
            self.w = w_new

    def predict(self, x):
        self.predict_proba = _sigmoid(np.dot(x, self.w))
        y_pred = self.predict_proba.copy()
        for _ in range(x.shape[0]):
            if y_pred[_] > 0.5:
                y_pred[_] = 1
            elif y_pred[_] < 0.5:
                y_pred[_] = 0
        return y_pred

model = LogisticRegression(
    tol=0.0001,
    max_iter=100
)

start = time.time()
model.fit(X,y)
end = time.time()
time_diff = end - start
print('実行時間',time_diff)
print("coef_",model.w)

pred = model.predict(X)
print("accuracy",accuracy_score(y, pred))
print("一個目のデータのクラス1に属する確率", model.predict_proba[0])
print("更新回数", model.iter)

参考文献

https://youtu.be/1ShUUblhb4s?si=iiDtIdBovXrHvum9

付録