ヤマカサのプログラミング勉強日記

プログラミングに関する日記とどうでもよい雑記からなるブログです。

Python機械学習プログラミング 達人データサイエンティストによる理論と実践 part. 3

Python機械学習プログラミング 達人データサイエンティストによる理論と実践

第3章

scikit-learn

ライブラリを使った機械学習を行います.テストデータを標準化して学習を行います.

from sklearn import datasets
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Perceptron
from sklearn.metrics import accuracy_score

iris = datasets.load_iris()

X = iris.data[:, [2, 3]]
y = iris.target

print('Class labels:', np.unique(y))

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1, stratify=y)

print('Label counts in y', np.bincount(y))
print('Label counts in y_train', np.bincount(y_train))
print('Label counts in y_test', np.bincount(y_test))

sc = StandardScaler()

sc.fit(X_train)
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)

ppn = Perceptron(eta0=0.1, random_state=1)
ppn.fit(X_train_std, y_train)

y_pred = ppn.predict(X_test_std)
print('Misclassified examples: %d' % (y_test != y_pred).sum())
print('Accuracy: %.3f' % accuracy_score(y_test, y_pred))
print('Accuracy: %.3f' % ppn.score(X_test_std, y_test))

ロジスティック関数の重みの学習

自称の起こりやすさを表した指標であるオッズ比は,

\begin{align} \dfrac{p}{1 - p} \end{align}

と表されます.シグモイド関数 $\phi (z)$ を

\begin{align} \phi (z) = \dfrac{1}{1 + e^{-z}} \end{align}

とします.コスト関数は,

\begin{align} J(\boldsymbol{w}) = \dfrac{1}{2} \sum_{i}( \phi(z) - \hat{y}) ^{2} \end{align}

であり,ロジスティックモデルの構築時に最大化したい尤度を

\begin{align} L(\boldsymbol{w}) = P(\boldsymbol{y} | \boldsymbol{x}; \boldsymbol{w}) = \prod_{n}^{i=1} (\phi (z^{(i)}))^{y(i)} (1 - \phi (z^{(i)}))^{(1 - y{(i)})} \end{align}

とします.計算のために対数尤度

\begin{align} I(\boldsymbol{w}) = \log L(\boldsymbol{w}) = \sum_{i = 1}^{n} [y^{(i)} \log(\phi (z^{(i)})) + (1-y^{(i)}) \log (1 - \phi (z^{(i)}))] \end{align}

を最大化するようにします.

勾配降下法による重みパラメータの更新

評価関数の最小化を行うために対数尤度の符号を反転したものを考えます.

\begin{align} I(\boldsymbol{w}) = - \sum_{i = 1}^{n} [y^{(i)} \log(\phi (z^{(i)})) + (1-y^{(i)}) \log (1 - \phi (z^{(i)}))] \end{align}

偏微分である

\begin{align} \dfrac{\partial I(\boldsymbol{w}) }{\partial \boldsymbol{w}}= \left [ \dfrac{\partial I(\boldsymbol{w})}{\partial w_{0}}, \cdots , \dfrac{\partial I(\boldsymbol{w})}{\partial w_{m}} \right ] ^{\mathrm{T}} \end{align}

について考えます.

\begin{align} \dfrac{\partial \phi (z)}{\partial z} = \phi (z) (1 - \phi (z)) \end{align}

\begin{align} \dfrac{\partial I(\boldsymbol{w}) }{\partial w_{j}} &=- \left ( \dfrac{y}{\phi(z)} - \dfrac{1 - y}{1 - \phi(z)}\right) \dfrac{\partial \phi (z)}{\partial w_{j}} \\ &= - \left ( \dfrac{y}{\phi(z)} - \dfrac{1 - y}{1 - \phi(z)}\right) \phi(z) (1 - \phi (z)) \dfrac{\partial z}{\partial w_{j}} \\ & = - (y - \phi(z)) x_{j} \end{align}

勾配降下法に基づくロジスティック回帰分析

class LogisticRegressionGD(object):
    """
    勾配降下法に基づくロジスティック回帰分析

    """

    def __init__(self, eta=0.01, n_iter=50, random_state=1):
        self.eta = eta
        self.n_iter = n_iter
        self.random_state = random_state

    def fit(self, X, y):
        rgen = np.random.RandomState(self.random_state)
        self.w_ = rgen.normal(loc=0.0, scale=0.01, size=1 + X.shape[1])
        self.cost_ = []

        for i in range(self.n_iter):
            net_input = self.net_input(X)
            output = self.activation(net_input)
            errors = (y - output)
            self.w_[1:] += self.eta * X.T.dot(errors)
            self.w_[0] += self.eta * errors.sum()
            cost = -y.dot(np.log(output)) - ((1 - y).dot(np.log(1 - output)))
            self.cost_.append(cost)
        return self

    def net_input(self, X):
        return np.dot(X, self.w_[1:]) + self.w_[0]

    def activation(self, z):
        """ロジスティックシグモイド活性化関数"""
        return 1.0 / (1.0 + np.exp(-np.clip(z, -250, 250)))

    def predict(self, X):
        return np.where(self.net_input(X) >= 0.0, 1, 0)

X_train_01_subset = X_train_std[(y_train == 0) | (y_train == 1)]
y_train_01_subset = y_train[(y_train == 0) | (y_train == 1)]

lrgd = logistic_regression_gd.LogisticRegressionGD(eta=0.05, n_iter=1000, random_state=1)
lrgd.fit(X_train_01_subset,
         y_train_01_subset)

plot_decision_regions(X=X_train_01_subset,
                      y=y_train_01_subset,
                      classifier=lrgd)

plt.xlabel('petal length [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')

plt.tight_layout()
plt.show()

ロジスティック回帰分析

感想

久しぶりに微分の計算をやったので商の微分とか忘れてました.