カテゴリー: 教養

SageMakerはローカルで使うことができるので、それを試してみた。 この記事を書くにあたって以下の公式の記事を参考にしています。 オンプレミス環境から Amazon SageMaker を利用する 機械学習のHelloWorld アヤメデータをschikit-learnの決定木分類器で学習して種類を予測する. 結構いろいろなところで機械学習のHelloWorldとして使われている例題を題材にしていく. sagemaker-python-sdk/scikit_learn_iris 既に SageMaker用のサンプルコードがあるので、 これをローカルで学習・推論できるように修正していく。 構成は以下の通り. 公式のブログの通り、SageMaker Notebook用に書かれた.ipynb を Local用に微修正するだけで動く。 SageMaker Notebookで動かすための.ipynb .ipynbから呼ぶschikit-learnコード SageMakerのサンプルはSageMakerのJupyterNotebookで動くように書かれているがが、 ちょっと修正するだけでローカルで動くようになる様子。(1つしか試してないけど) 前準備 学習と推論をローカルで行うが、そのために裏でDockerのコンテナが走る。 ローカルコンピュータ用にDockerをインストールしておく必要がある。 CredentialsとIAM 以下が必要。 AmazonSageMakerFullAccess 権限をもった IAM ユーザの Credential AmazonSageMakerFullAccess の IAM ロール ローカルコードからAWSリソースにアクセスするために aws configure を使って設定する。 Credentialsが書かれたcsvをダウンロードし aws configure の応答に答えていく。 $ pip install awscli --upgrade --user $ aws configure AWS Access Key ID [None]: ****************** AWS Secret Access Key [None]: ************************************ Default region name [None]: ap-northeast-1 Default output format [None]: json SageMaker PythonSDKインストール SageMaker PythonSDKをインストールする。 実行するコードに応じてSDKのバージョンを指定することができる。 $ pip install -U sagemaker >=2.15 バージョンを指定しない場合は以下の通り。 $ pip install sagemaker ローカルのJupyter Notebookでファイルを修正 scikit_learn_estimator_example_with_batch_transform.ipynb を ローカルのJupyter Notebookで修正していく。 SageMaker ローカルSessionを開始 SageMakerを想定したコードは以下。 # S3 prefix prefix = \"Scikit-iris\" import sagemaker from sagemaker import get_execution_role sagemaker_session = sagemaker.Session() # Get a SageMaker-compatible role used by this Notebook Instance. role = get_execution_role() それをローカルで動かすために以下のように修正する localSession()というセッションが用意されているのでそれを使用する。 ローカルでは get_execution_role()では ロールを取得できないので直接ロールのARNを指定する。 # S3 prefix prefix = \"Scikit-iris\" import sagemaker from sagemaker import get_execution_role # LocalSession()を使用する sagemaker_session = sagemaker.local.LocalSession() # sagemaker.Session()から変更 # Get a SageMaker-compatible role used by this Notebook Instance. # ローカルでは get_execution_role()は使えない。直接ロールのARNを指定する。 # role = get_execution_role() role = \'arn:aws:iam::(12桁のAWSアカウントID):role/(ロール名)\' 学習用データの準備 (変更なし) 学習用データが巨大であればS3にデータを準備する(と書かれている). アヤメデータは軽量なので、ローカルファイルに保存する。 import numpy as np import os from sklearn import datasets # Load Iris dataset, then join labels and features iris = datasets.load_iris() joined_iris = np.insert(iris.data, 0, iris.target, axis=1) # Create directory and write csv os.makedirs(\"./data\", exist_ok=True) np.savetxt(\"./data/iris.csv\", joined_iris, delimiter=\",\", fmt=\"%1.1f, %1.3f, %1.3f, %1.3f, %1.3f\") その後、用意したローカルデータをSageMaker Python SDKに食わせる。 WORK_DIRECTORY = \"data\" train_input = sagemaker_session.upload_data( WORK_DIRECTORY, key_prefix=\"{}/{}\".format(prefix, WORK_DIRECTORY) ) Scikit learn Estimator scikit-learnの機械学習は以下の3段構成になっている。 Estimator: 与えられたデータから学習(fit)する Transformer: 与えられたデータを変換(transform)する Predictor: 与えられたデータから結果を予測(Predict)する SageMakerは機械学習プラットフォームであって、かなり多くのライブラリや手法がサポートされている。 その中で、scikit-learnもサポートされていて、SKLearn Estimatorとして使用できる。 要は、schikit-learn のI/Fに準じたコードを SageMaker に内包することができる。 SKLearn Estimatorに scikit-learn コードを食わせると SageMakerから SKLearn インスタンスとして操作できる. 例えば、.ipynbで以下のように書く. from sagemaker.sklearn.estimator import SKLearn FRAMEWORK_VERSION = \"0.23-1\" script_path = \"scikit_learn_iris.py\" sklearn = SKLearn( entry_point=script_path, framework_version=FRAMEWORK_VERSION, instance_type=\"local\", role=role, sagemaker_session=sagemaker_session, hyperparameters={\"max_leaf_nodes\": 30}, ) SKLearnにentry_pointとして渡しているのがscikit-learnのコード本体。 内容は以下。普通の決定木分類のコードにSageMakerとのIFに関わるコードが追加されている。 実行時引数として、SM_MODEL_DIR、SM_OUTPUT_DATA_DIR、SM_CHANNEL_TRAINが渡される。 fitで学習した結果(つまり係数)をシリアライズしSM_MODEL_DIRに保存する。 model_fnでは、SM_MODEL_DIRにシリアライズされた係数をデシリアライズし、 scikit-learnの決定木分類木オブジェクトを返す。 # Copyright 2018 Amazon.com, Inc. or its affiliates. All Rights Reserved. # # Licensed under the Apache License, Version 2.0 (the \"License\"). # You may not use this file except in compliance with the License. # A copy of the License is located at # # http://www.apache.org/licenses/LICENSE-2.0 # # or in the \"license\" file accompanying this file. This file is distributed # on an \"AS IS\" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either # express or implied. See the License for the specific language governing # permissions and limitations under the License. from __future__ import print_function import argparse import os import joblib import pandas as pd from sklearn import tree if __name__ == \"__main__\": parser = argparse.ArgumentParser() # Hyperparameters are described here. In this simple example we are just including one hyperparameter. parser.add_argument(\"--max_leaf_nodes\", type=int, default=-1) # Sagemaker specific arguments. Defaults are set in the environment variables. parser.add_argument(\"--output-data-dir\", type=str, default=os.environ[\"SM_OUTPUT_DATA_DIR\"]) parser.add_argument(\"--model-dir\", type=str, default=os.environ[\"SM_MODEL_DIR\"]) parser.add_argument(\"--train\", type=str, default=os.environ[\"SM_CHANNEL_TRAIN\"]) args = parser.parse_args() # Take the set of files and read them all into a single pandas dataframe input_files = [os.path.join(args.train, file) for file in os.listdir(args.train)] if len(input_files) == 0: raise ValueError( ( \"There are no files in {}.n\" + \"This usually indicates that the channel ({}) was incorrectly specified,n\" + \"the data specification in S3 was incorrectly specified or the role specifiedn\" + \"does not have permission to access the data.\" ).format(args.train, \"train\") ) raw_data = [pd.read_csv(file, header=None, engine=\"python\") for file in input_files] train_data = pd.concat(raw_data) # labels are in the first column train_y = train_data.iloc[:, 0] train_X = train_data.iloc[:, 1:] # Here we support a single hyperparameter, \'max_leaf_nodes\'. Note that you can add as many # as your training my require in the ArgumentParser above. max_leaf_nodes = args.max_leaf_nodes # Now use scikit-learn\'s decision tree classifier to train the model. clf = tree.DecisionTreeClassifier(max_leaf_nodes=max_leaf_nodes) clf = clf.fit(train_X, train_y) # Print the coefficients of the trained classifier, and save the coefficients joblib.dump(clf, os.path.join(args.model_dir, \"model.joblib\")) def model_fn(model_dir): \"\"\"Deserialized and return fitted model Note that this should have the same name as the serialized model in the main method \"\"\" clf = joblib.load(os.path.join(model_dir, \"model.joblib\")) return clf 学習 SageMaker(またはローカル)の.ipynb は、SKLearnインスタンスに対して fit() を実行するだけで良い。 sklearn.fit({\"train\": train_input}) 推論 なんと、推論はWebインターフェースになっている。 推論コンテナ内でnginxが動作し、PythonWebAppが wsgi(gunicorn) を介してnginxから入力/応答する。 SageMaker(またはローカル)の.ipynbからは、SKLearnインスタンスに対してdeploy()を実行する。 推論コンテナへのインターフェースとなるインスタンスが生成され、後はこのインスタンスに対してpredict()を呼ぶ。 推論用にデータを集めて predict()を実行する例。 テストデータと推論の結果を並べて表示している。 うまくいっていれば同じになるはず。 (こんな風に訓練データとテストデータを拾って良いのかはさておき...) predictor = sklearn.deploy(initial_instance_count=1, instance_type=\"local\") import itertools import pandas as pd shape = pd.read_csv(\"data/iris.csv\", header=None) a = [50 * i for i in range(3)] b = [40 + i for i in range(10)] indices = [i + j for i, j in itertools.product(a, b)] test_data = shape.iloc[indices[:-1]] test_X = test_data.iloc[:, 1:] test_y = test_data.iloc[:, 0] print(predictor.predict(test_X.values)) print(test_y.values) /invocationsというURLに対してPOSTリクエストが発行されている。 応答は以下の通り、 テストデータの説明変数と、predict()の結果得られた値が一致していそう。 hqy7i6eoyi-algo-1-vq8df | 2021-05-22 16:26:50,617 INFO - sagemaker-containers - No GPUs detected (normal if no gpus installed) [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. 2. 2. 2. 2. 2. 2. 2. 2.]hqy7i6eoyi-algo-1-vq8df | 172.23.0.1 - - [22/May/2021:16:26:51 +0000] \"POST /invocations HTTP/1.1\" 200 360 \"-\" \"python-urllib3/1.26.4\" [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. 2. 2. 2. 2. 2. 2. 2. 2.] まとめ SageMaker用の機械学習のHelloWorldをローカルで動かしてみた。 （かなり雑だけども), SageMakerのサンプルコードをちょっと修正するだけでローカルで動くことがわかった。

[mathjax] 分散と標準偏差を計算しやすく変形できる。 いちいち偏差(x_i-bar{x})を計算しておかなくても、2乗和(x_i^2)と平均(bar{x})がわかっていればOK。 begin{eqnarray} s^2 &=& frac{1}{n} sum_{i=1}^n (x_i - bar{x})^2 \\ &=& frac{1}{n} sum_{i=1}^n ( x_i^2 -2 x_i bar{x} + bar{x}^2 ) \\ &=& frac{1}{n} ( sum_{i=1}^n x_i^2 -2 bar{x} sum_{i=1}^n x_i + bar{x}^2 sum_{i=1}^n 1) \\ &=& frac{1}{n} ( sum_{i=1}^n x_i^2 -2 n bar{x}^2 + nbar{x}^2 ) \\ &=& frac{1}{n} ( sum_{i=1}^n x_i^2 - nbar{x}^2 ) \\ &=& frac{1}{n} sum_{i=1}^n x_i^2 - bar{x}^2 \\ s &=& sqrt{frac{1}{n} sum_{i=1}^n x_i^2 - bar{x}^2 } end{eqnarray} 以下みたい使える。平均と標準偏差と2乗和の関係。 begin{eqnarray} sum_{i=1}^n (x_i - bar{x})^2 &=& sum_{i=1}^n x_i^2 - nbar{x}^2 \\ ns^2 &=& sum_{i=1}^n x_i^2 - nbar{x}^2 \\ sum_{i=1}^n x_i^2 &=& n(s^2 + bar{x} ) end{eqnarray}

[mathjax] 参考書籍を読んでいく. 今回はロジスティック回帰. 未知のサンプルのクラスの所属確率を見積もることができる, という重要なことが書かれている. また重みの学習は勾配降下法のフレームワークの中で統一されていて一連の理解が進む. 一度は理解しておいた方が良い内容だと思った. パーセプトロンとADALINE 線形和(z=w^T)を考えたとき 「サンプル(x)が決定境界の上にある」 というのは(w^Tx=0). (w^T x > 0), (w^T x < 0) を決定境界の 「こちら側」 か 「向こう側」 か対応させる. パーセプトロンのアイデアは以下の通り. 線形和 (z=w^T x)を決定境界とし, (phi(z)in (-1,+1))を予測ラベルとする 予測が誤った回数, つまり真のラベルと予測ラベルの差の合計がゼロになるように(w)を学習する (phi(z))はステップ関数で連続ではないので解析的な式の変形ができない 決定境界により完全に分離できなければ学習が収束することはない 決定境界から離れている度合いを連続で凸な関数として定義し, その関数が極小となる(w)を求めていこう, というアイデアに拡張したのがADALINE. (phi(z) = z). コスト関数(J(w) = frac{1}{2} sum_i left( hat{y^{(i)} - phi ( z^{(i)})} right)^2 ) トレーニングセット全体から(J(w))の極小化を考えることが可能. (勾配降下法) トレーニングセットの一部を使って(J(w))の極小化を考え, 他のトレーニングセットで反復することも可能. (確率的勾配降下法) ロジスティック回帰 決定境界を表す線形和(z=w^T x)について(phi(z) in (mathbb{R}| 0 leq z leq +1))を考える. 線形和(z)がクラス1に分類される確率を表す(phi(z))を考える. (phi(z))をどう作るのか 線形和を(0)から(1)の実数全体に変換する関数を考える. オッズ比の対数をとったロジット関数(logit(p))を使う. ここで (p)は(x)が与えられたときにサンプルデータがクラス1に属する確率. begin{eqnarray} logit(p(y=1 | x)) &=& log frac{p}{1-p} \\ &=& w_0 x_0 + w_w x_w + cdots + w_m x_m \\ &=& sum_{i=0}^m w_i x_i \\ &=& w^T x end{eqnarray} つまり, 線形和(z = w^T x in (mathbb{R} | 0 leq z leq 1) ) (p)は未知であるということが重要. この式から(p)を予測していく. (f(p) = log frac{p}{1-p}) の逆関数が分かれば良い. begin{eqnarray} f(p) &=& log frac{p}{1-p} \\ e^{f(p)} &=& frac{p}{1-p} \\ (1-p)e^{f(p)} &=& p \\ e^{f(p)} &=& p(1+e^{f(p)}) \\ p &=& frac{e^{f(p)}}{1+e^{f(p)}} \\ p &=& frac{1}{1+e^{-f(p)}} end{eqnarray} (f(p)=z)としていたので以下のようになる. (zとpが混在していてたぶん数式としての表現は間違ってる.. ) begin{eqnarray} p &=& frac{1}{1+e^{-z}} end{eqnarray} (p)は確率なので(p in (mathbb{R}|0 leq p leq 1)). (phi(z) in (mathbb{R}|0 leq p leq 1) )とすると, begin{eqnarray} phi(z) = frac{1}{1+e^{-z}} end{eqnarray} (z-phi(z))をGeoGebraでプロットすると以下のような感じ. (z)を大きくしていくと(phi(z))は限りなく(1)に近づく. (z)を小さくしていくと(phi(z))は限りなく(0)に近づく. (phi(z)=)の出力は正事象が起こる確率であると解釈される. (phi(z)=0.8)である場合, 80%の確率で正事象が起こり,20%の確率で正事象が起こらないことを表す. トレーニングデータからパラメタを学習し, ある(phi(z))が得られたとして, 未知のサンプルのクラスの所属関係を確率で見積もることができる. 天気予報で晴れの確率70%、晴れでない確率30%のように. 予測の際に,クラスの所属関係の確率を見積もることができるのが神がかるレベルで重要!! もし(phi(z))を2値分類に当てはめるのであれば, (phi(z)=0.5)を境に上半分を(1)、下半分を(0)と予測するものとする. begin{eqnarray} hat{y} = begin{cases} 1 & phi(z) ge 0.5 \\ 0 & phi(z) lt 0.5 end{cases} end{eqnarray} それって, (z=0)を境に左半分, 右半分に割り付けるのと同じなので, begin{eqnarray} hat{y} = begin{cases} 1 & z ge 0 \\ 0 & z lt 0 end{cases} end{eqnarray} ロジスティック回帰の重みの学習 では(w)をどうやって学習するか. 式をイジるだけで面倒なだけで,これ以上は洞察得られないなと思うので, 流れだけ書いてみる. 学生時代に形態素解析っぽい何かのアルゴリズムを考えたことがあって, 形態素どうしのベターな接続を求め際に, 接続コストをコーパスから作った確率で表して, 尤度を最大化する問題として解いた記憶があって, 今更合点が行くという謎. 勉強不足で尤度最大化のアイデアがあってそうした訳ではないので, もっと勉強していたらよかったな..と. ADALINEの説明の中で, 真の値と予測値の差の2乗誤差をコスト関数として, 連続で凸な関数の極小化の問題として(w)を求めていった. ロジスティック回帰の場合, 予測は(phi(z))も連続で凸な関数なので, 同じロジックで(w)を求めることもできる. コスト関数(J(w))は同様に以下. begin{eqnarray} J(w) &=& sum_i frac{1}{2} left( phi(z^{(i)}) - y^{(i)} right)^2 end{eqnarray} (phi(z)=frac{1}{1+e^{-z}})だと, これを解析的に微分するのは不可能. トレーニングデータが(m)個あった場合に(i)番目のトレーニングデータにおける確率は以下. (x)も(w)もパラメタですよ,って書き方. begin{eqnarray} P(y|x;w)=P(y^{(i)}|x^{(i)};w) end{eqnarray} 求めたいのは, 全てのトレーニングデータについて(P(y))の積が最大になる(w). 積,つまり尤度を書き下すと以下のようになる. begin{eqnarray} L(w) &=& prod_{i=1}^m left( P(y^{(i)}) | x^{(i)};w right) \\ &=& prod_{i=1}^m left( phi(z^{(i)}) right)^{y^{(i)}} + left( 1-y^{(i)}right) log left( 1-phi(z^{(i)}) right) end{eqnarray} 尤度そのものの最大化するよりも対数をとったものを最大化する方が楽になる. (対数をとると指数の肩が定数倍になり積が和になる. 相当簡単になる. ) begin{eqnarray} l(w) &=& log L(w) \\ &=& sum_{i=1}^{n} left( y^{(i)} log( phi(z^{(i)}) ) + (1-y^{(i)}) log (1-phi(z^{(i)})) right) end{eqnarray} これを最大化するということは, これに(-1)を掛けたものを最小化するということ. そもそもコスト関数(J(w))を立てて最小化しようとしていたので, これを(J(w))と置いてしまう. begin{eqnarray} J(w) = sum_{i=1}^{n} left( -y^{(i)} log( phi(z^{(i)}) ) - (1-y^{(i)}) log (1-phi(z^{(i)})) right) end{eqnarray} ADALINEの勾配降下法で, コスト関数を最小にするように(w)を更新していった. つまり, (w := w - Delta w)という更新をしていった. ロジスティック回帰において, コスト関数を対数尤度で書き直していて, 対数尤度を最大化するように(w)を更新する. つまり, (w := w + Delta w)という更新をおこなう. 対数尤度を最大化する(w)を更新する話は, コスト関数を最小化する(w)を更新する話と同じ, なので, 結局のところ勾配降下法と同じ式となる. begin{eqnarray} Delta w_j &:=& - eta frac{partial J}{partial w_j} end{eqnarray} 尤度関数の偏導関数はシグモイド関数の偏導関数を使って以下のようになるらしい. (省略) シグモイド関数の偏導関数は以下 begin{eqnarray} frac{partial phi(z)}{partial z} &=& frac{1}{partial z} frac{1}{1-e^{-z}} \\ &=& frac{1}{1+e^{-1}} left( 1- frac{1}{1+e^{-z}} right) \\ &=& phi(z) (1-phi(z)) end{eqnarray} 対数尤度の偏導関数はシグモイド関数の偏導関数を使って以下. begin{eqnarray} frac{partial l(w)}{partial w_j} &=& left( y frac{1}{phi(z)} - (1-y) frac{1}{1-phi(z)} right) frac{partial}{partial w_j} phi(z) \\ &=& left( yfrac{1}{phi(z)} - (1-y) frac{1}{1-phi(z)} right) phi(z)(1-phi(z)) frac{partial}{partial w_j} z \\ &=& left( y(1-phi(z))-(1-y)phi(z) right) x_j \\ &=& left( y-phi(z) right) x_i end{eqnarray} コスト関数の最小化と対数尤度関数の最大化が同じ,というロジックが効いて, ADALINEの更新と全く同じになる. begin{eqnarray} w_j &:=& - eta frac{partial J}{partial w_j} \\ &=& eta sum_{i=1}^n left( y^{(i)} -phi(z^{(i)}) right) x_j^{(i)} end{eqnarray} 奇跡的..

[mathjax] 参考にした書籍でこの順序で誘導されていて理解しやすかったです. パーセプトロンによる一番簡単な教師あり学習を理解する ADALINEにより学習を凸で連続なコスト関数の最小化問題として捉える パーセプトロンの学習 2値のラベル((+1),(-1))が付与されているサンプルデータが与えられたとする. ラベル値が(-1)であるデータと, ラベル値が(+1)であるデータに分類できる. それらの境界が線形和で表される問題である場合, その境界を求めることができれば, 未知のデータに対してその境界の(+1)側か(-1)側かを判定できる. (m)次のサンプルデータと重みの線形和を考える. [ boldsymbol{z} = boldsymbol{w}^T boldsymbol{x} = w_1 x_1 + w_2 x_2 + cdots x_m x_m ] (z)に関するステップ関数を用意する. こうすることで (z)-(phi(z)) 空間において(z=theta)を境に線形和(z)が(-1),(+1)に分かれることを表現できる. [ varphi(z) = begin{cases} 1 & (z gt theta) \\ -1 & (z leq theta) end{cases} ] (theta)を左辺に移動し, (x_0 = 1), (w_0 = -theta) とすると, 線形和に(theta)の項を組み入れることができる. 新しい線形和(z\')が(0)になる箇所を境にステップ関数の応答値が切り替わる. [ varphi(z) = begin{cases} 1 & (z-theta gt 0) \\ -1 & (z-theta leq 0) end{cases} ] [ z\' = z - theta = w_0 x_0 + w_1 x_1 + w_2 x_2 + cdots x_m x_m ] あるパラメータ(boldsymbol{w})が存在するとして, データ(boldsymbol{x})について, (boldsymbol{w}^T boldsymbol{x}=0)を境界にステップ関数の応答値が切り替わる. 未知のデータ(x)に対して, (boldsymbol{w}^T boldsymbol{x})はラベルの予測値(+1), (-1)を応答する. 予測値(varphi(w^Tx))が正解であるサンプルデータに含まれる(y^{i})と同じとなる(w)を探したい. それを探していく. 最初,(w)に初期値を設定し, 以下の手続きによって(w)を更新していく. 上付きの添字はサンプルデータ内の順序数を表している. (y^{(i)})は(i)番目のサンプルデータ. 右下の添字は該当サンプルデータの各次元を表している. (y^{(i)_m})は(y^{(i)})の(m)番目の要素. 現在の(w)を使って予測値を計算する. (hat{y}=w^T x) (w)を更新する. (Delta w_j := eta(y^{i}-hat{y})x^{i}_j) 1個のサンプルデータで1回(w)が更新される. そもそもデータセットを綺麗に線形和が表す決定境界で分離できないような状態だと、 何度やっても予測したクラスラベルと真のクラスラベルの差が埋まらない。 あっちを立てればこっちが立たず、みたいになる。 予測したクラスラベルと真のクラスラベルを比較するこの方法だと、 (Delta w_j := eta(y^{i}-hat{y})x^{i}_j) が良い値なのか悪い値なのか、繰り返さないとわからない。 予測と真の値の比較を凸で連続な関数で表すことができれば、その関数の凸の底に向かうやり方で、 「今より良い次の値」に更新する方法を解析的に求めることができて都合が良い。 真の値と予測した値の差が「損失」とか「コスト」とか呼んで、 「コスト関数」を最小化するパラメータが良いパラメータなんだろう、という話が進んでいく。 ADALINEの学習 線形和をステップ関数に変換する部分があった。 こうして(w^T x varphi(w^T x))空間で (w^T x)が(-1),(+1)に分かれる(w)を求めようとした. [ varphi(z) = begin{cases} 1 & (z-theta gt 0) \\ -1 & (z-theta leq 0) end{cases} ] そうではなく以下の恒等式(左辺と右辺が同じ)を使い,真の値と予測した値の差の求め方を変えると, 差が凸で連続な関数で表せるようになり (解析的に微分できるようになり), 一番小さいところは? という話がしやすくなる. [ varphi(z)=z ] 最小二乗法のときやったやつで、差の2乗を足し合わると符号がキャンセルされてよくて、差を式で表せる。 [ J(w) = sum_i left( y^{(i)} - varphi(w^T x^{(i)}) right)^2 ] こういうのを誤差平方和と言うようで(frac{1}{2})倍するらしい. [ J(w) =frac{1}{2} sum_i left( y^{(i)} - varphi(w^T x^{(i)}) right)^2 ] (J(w))がなんで連続で凸なのかは知ったことではないが、 連続で凸であれば(J\'(w)=0)を解くと極小となる(w)が求まるのは高校生のときにやった話. (J(w))の接線(多次元だと接っする面)の傾きが(J\'(w))。 ただ(w)には次数があって傾きは以下(それぞれの次数毎の極限). [ nabla J(w) = frac{partial J(w)}{partial w_j} ] (w J(w))空間で(nabla J(w))は(w)における傾きなので, (w)に(-nabla J(w) cdot 1)を足すと, (J(w))が極小になる箇所に近づく. どれだけ行けば良いかわからないので(-nabla J(w) cdot eta)を足すことにする. (w)を以下の通り更新する. begin{eqnarray} w &:=& w + nabla J(w) \\ &:=& w -nabla J(w) cdot eta end{eqnarray} ちなみに(nabla J(w))は式をこねくり回すと計算できる. 最終的に更新は以下の通りとなる. [ w := w - eta sum_i left( y^{(i)} -varphi left( z^{(i)} right) right) x_j^{(i)} ] パーセプトロンの更新とそっくりな形が出てきて感動する... あっちは(varphi(z))が不連続なステップ関数でこっちは連続な関数. 式をこねくり回す際に, (varphi(z) = z)の恒等式の関係を使っていない. 何か連続な関数であればこれが成り立つ. 形式上,例えば(varphi(z))を以下(シグモイド)とすることでロジスティック回帰 (本当か? 次回確認.). [ varphi (z) = frac{1}{1+e^{-z}} ]

ある値の近傍を取る状態から別の値の近傍を取る状態へ遷移する時系列データについて、 どちらの状態でもない過渡状態を経て状態が遷移するみたい。 今回知りたいのは理論ではなく、定常状態に遷移したことをいかにして検知するかという物理。 まぁいかようにでも検知できなくもないけれども、非定常状態の継続時間は不明であるという。 なるべく遅れなしに定常状態への遷移を検知したい。 ARモデルとかMAモデルとかの出番ではない。 DeepLearningとかやってる暇はないw 過去の微小時間内に含まれるデータだけから次が定常状態なのかを検知したい。 前状態と後状態のベースラインがどんな値であっても、 それを考慮することなしに過渡応答の後にくる定常状態の先頭を検知できる、という特徴がほしい。 移動平均だとベースラインがケースごとに変わるので扱いづらい。 微分であればベースラインがゼロに揃うのが良いな、と思った。 結局、今度はゼロへの収束を検知する問題に差し代わるだけだけども、 前ラインと後ラインの2つも見ないでも、確実に楽にはなっている。 微分するとノイズが増えるが...。 微分値が上下の閾値を超えたラインから先、 ゼロ収束の判定条件が一定時間継続した最初のデータポイントがソコ。 これは過渡状態の終了検知なだけだから、それとRawデータの変化をみる。

[mathjax] 各地点において関数の値を最大にするベクトル((frac{partial f}{partial x_0},frac{partial f}{partial x_1}))を全地点に対して計算したものを勾配とかいう。 ある地点において、このベクトルの方向に向かうことにより最も関数の値を大きくする。 で、今後のために正負を反転して関数の値を最小にするベクトルを考えることにした。 関数の値を小さくする操作を繰り返していけば、いずれ\"最小値\"が見つかるはず。 というモチベを続けるのが勾配降下法。学習率(eta)を使って以下みたいに書ける。。 begin{eqnarray} x_0 = x_0 - eta frac{partial f}{partial x_0} \\ x_1 = x_1 - eta frac{partial f}{partial x_1} end{eqnarray} ということで(f(x_0,x_1)=x_0^2+x_1^2)の最小値を初期値((3.0,4.0))、 学習率(eta=0.1)に設定して計算してみる。 import numpy as np def numerical_gradient(f, x): h = 1e-4 grad = np.zeros_like(x) for idx in range(x.size): tmp_val = x[idx] x[idx] = tmp_val + h fxh1 = f(x) x[idx] = tmp_val - h fxh2 = f(x) grad[idx] = (fxh1 - fxh2) / (2*h) x[idx] = tmp_val return grad def gradient_descent(f, init_x, lr=0.01, step_num=100): x = init_x for i in range(step_num): grad = numerical_gradient(f,x) x -= lr * grad return x def function2(x): return x[0]**2 + x[1]**2 init_x = np.array([-3.0, 4.0]) v = gradient_descent(function2, init_x=init_x, lr=0.1, step_num=100) v # array([-6.11110793e-10, 8.14814391e-10]) ((0,0))に収束した。 ニューラルネットワークの勾配 損失関数を重みパラメータで微分する。以下みたいな感じ。 損失関数の大小を見るとして、例えば(w_{11})以外の重みを固定したとして(w_{11})をわずかに 増やしたときに損失関数の値がどれだけ大きくなるか。 損失関数の値はパラメータ(W)と入力(x)から決まるベクトルだけれども、それぞれ乱数と入力値が設定されている。 begin{eqnarray} W= begin{pmatrix} w_{11} & w_{12} & w_{13} \\ w_{21} & w_{22} & w_{23} end{pmatrix}, frac{partial L}{partial W}= begin{pmatrix} frac{partial L}{partial w_{11}} & frac{partial L}{partial w_{12}} & frac{partial L}{partial w_{13}} \\ frac{partial L}{partial w_{21}} & frac{partial L}{partial w_{22}} & frac{partial L}{partial w_{23}} end{pmatrix} end{eqnarray} 重み(W)が乱数で決まるネットワークがあるとする。このネットワークは入力と重みの積を出力 として返す。出力はSoftmaxを経由するとする。 ネットワークの出力と教師データのクロスエントロピー誤差を誤差として使う。 その前に、数値微分関数を多次元対応する。 普通、配列の次元が(n)個になると(n)重ループが必要になるけれども、 Numpy.nditer()を使うと(n)乗ループを1回のループにまとめることができる。 下のmulti_indexが((0,0),(0,1),(0,2),(1,0),(1,1),(1,2))みたいに イテレータが(n)次のタプルを返す。反復回数はタプルの要素数の直積。 Numpy配列にそのタプルでアクセスすることで晴れて全ての要素にアクセスできる。 def numerical_gradient_md(f, x): h = 1e-4 grad = np.zeros_like(x) it = np.nditer(x, flags=[\'multi_index\'], op_flags=[\'readwrite\']) while not it.finished: idx = it.multi_index tmp_val = x[idx] x[idx] = tmp_val + h fxh1 = f(x) # f(x+h) x[idx] = tmp_val - h fxh2 = f(x) # f(x-h) grad[idx] = (fxh1 - fxh2) / (2*h) x[idx] = tmp_val # 値を元に戻す it.iternext() return grad 初期値(x=(0.6,0.9))、教師データ(t=(0,0,1))をネットワークに入力する。 predict()は(1 times 3)を返す。 それをSoftmax()を通して、(t)とのクロスエントロピー誤差を求めたものが以下。 import numpy as np def cross_entropy_error(y, t): if y.ndim == 1: t = t.reshape(1, t.size) y = y.reshape(1,y.size) batch_size = y.shape[0] delta = 1e-7 return -np.sum( t * np.log( y + delta)) / batch_size def softmax(x): c = np.max(x) return np.exp(x-c) / np.sum(np.exp(x-c)) import sys, os sys.path.append(os.pardir) import numpy as np class simpleNet: def __init__(self): self.W = np.random.randn(2,3) def predict(self, x): return np.dot(x, self.W) def loss(self, x, t): z = self.predict(x) y = softmax(z) loss = cross_entropy_error(y, t) return loss net = simpleNet() x = np.array([0.6, 0.9]) p = net.predict(x) t = np.array([0, 0, 1]) net.loss(x, t) # 0.9463818740797788 このlossを(W)で微分したのが以下。 あえてパラメータ(W)を引数にとり損失関数の値を計算する(f(W))を定義することで、 数値微分が何と何の演算なのかをわかりやすくしている。 実際は(f(W))は(W)とは関係なく(x)と(t)だけから結果を返すけれども、 損失関数(f(W))を(W)で微分するという操作が自明になるようにコードを合わせている。 def f(W): return net.loss(x, t) dW = numerical_gradient_md(f, net.W) dW # array([[ 0.07627371, 0.49923236, -0.57550607], # [ 0.11441057, 0.74884853, -0.8632591 ]]) 結果の解釈 上記の(w),(W),(t)から(frac{partial L}{partial W})が求まった。 損失関数が何か複雑な形をしているという状況で、 (frac{partial L}{partial w_{11}})は(w_{11})がわずかに動いたときに損失関数の値が変化する量を表している。 それが(w_{11})から(w_{32})まで6個分存在する。 begin{eqnarray} frac{partial L}{partial W} = begin{pmatrix} frac{partial L}{partial w_{11}} & frac{partial L}{partial w_{21}} & frac{partial L}{partial w_{31}} \\ frac{partial L}{partial w_{12}} & frac{partial L}{partial w_{22}} & frac{partial L}{partial w_{32}} end{pmatrix} = begin{pmatrix} 0.07627371 & 0.49923236 & -0.57550607 \\ 0.11441057 & 0.74884853 & -0.8632591 end{pmatrix} end{eqnarray}

[mathjax] 2変数関数(f(x_0,x_1))を各変数で偏微分する。 地点((i,j))におけるベクトル((frac{partial f(x_0,j)}{partial x_0},frac{partial f(i,x_1)}{partial x_1}))を全地点で記録していき、ベクトル場を得る。 このベクトル場が勾配(gradient)。 (f(x_0,x_1)=x_0^2+x_1^2)について、(-4.0 le x_0 le 4.0)、(-4.0 le x_1 le 4.0)の範囲で、 勾配を求めてみる。また、勾配を可視化してみる。 まず、2変数関数(f(x_0,x_1))の偏微分係数を求める関数の定義。 ((3.0,3.0))の偏微分係数は((6.00..,6.00..))。 def numerical_gradient(f, x): h = 10e-4 grad = np.zeros_like(x) for idx in range(x.size): tmp_val = x[idx] x[idx] = tmp_val + h fxh1 = f(x) x[idx] = tmp_val - h fxh2 = f(x) grad[idx] = (fxh1 - fxh2) / 2*h x[idx] = tmp_val return grad def function2(x): return x[0]**2 + x[1]**2 p = np.array([3.0,3.0]) v = numerical_gradient(function2, p) v # array([6.e-06, 6.e-06]) (-4.0 le x_0 le 4.0)、(-4.0 le x_1 le 4.0)の範囲((0.5)刻み)で偏微分係数を求めて、 ベクトル場っぽく表示してみる。matplotlibのquiver()は便利。 各地点において関数の値を最も増やす方向が表示されている。 w_range = 4 dw = 0.5 w0 = np.arange(-w_range, w_range, dw) w1 = np.arange(-w_range, w_range, dw) wn = w0.shape[0] diff_w0 = np.zeros((len(w0), len(w1))) diff_w1 = np.zeros((len(w0), len(w1))) for i0 in range(wn): for i1 in range(wn): d = numerical_gradient(function2, np.array([ w0[i0], w1[i1] ])) diff_w0[i1, i0], diff_w1[i1, i0] = (d[0], d[1]) plt.xlabel(\'$x_0$\',fontsize=14) #x軸のラベル plt.ylabel(\'$x_1$\',fontsize=14) #y軸のラベル plt.xticks(range(-w_range,w_range+1,1)) #x軸に表示する値 plt.yticks(range(-w_range,w_range+1,1)) #y軸に表示する値 plt.quiver(w0, w1, diff_w0, diff_w1) plt.show() 値が大きい方向に矢印が向いている。例えば((-3.0,3.0))における偏微分係数は((-6.0,6.0))。 左上方向へのベクトル。 参考にしている本にはことわりが書いてあり、勾配にマイナスをつけたものを図にしている。 その場合、関数の値を最も減らす方向が表示されることになる。 各地点において、この勾配を参照することで、どちらに移動すれば関数の値を最も小さくできるかがわかる。

引き続き、ゼロDの写経を続ける。今回は、学生の頃が懐かしい懐ワード、数値微分。 全然Deepに入れないけれどおっさんの復習。解析的な微分と数値微分が一致するところを確認してみる。 昔と違うのは、PythonとJupyterNotebookで超絶簡単に実験できるし、 こうやってWordPressでLaTeXで記事を書いたりできる点。 [mathjax] まず、微分の基本的な考え方は以下の通り。高校数学の数3の範囲。 begin{eqnarray} frac{df(x)}{fx} = lim_{hrightarrow infty} frac{f(x+h)-f(x)}{h} end{eqnarray} 情報系学科に入って最初の方でEuler法とRunge-Kutta法を教わってコードを 書いたりレポート書いたりする。懐すぎる..。 または、基本情報の試験かなんかで、小さい値と小さい値どうしの計算で発生する問題が現れる。 ゼロDにはこの定義を少し改良した方法が載っている。へぇ。 begin{eqnarray} frac{df(x)}{fx} = lim_{hrightarrow infty} frac{f(x+h)-f(x-h)}{2h} end{eqnarray} 写経なので、がんばって数値微分を書いて動かしてみる。 簡単な2次関数(f(x))。 begin{eqnarray} f(x) &=& x^2 - 5x +3 \\ f\'(x) &=& 2x - 5 end{eqnarray} def numerical_diff(f, x): h = 10e-4 return (f(x+h) - f(x-h)) / (2*h) (f(x))と、(x=2.5)のところの(f\'(x))をmatplotlibで書いてみる。懐い... import matplotlib.pyplot as plt import numpy as np def f(x): return x**2 - 5*x + 3 x = np.arange(-10, 10, 0.1) y = f(x) dy = numerical_diff(f,x) plt.plot(x, y) x1 = -2.5 dy1 = numerical_diff(f, x1) y1 = f(x1) # y-y1 = dy1(x-x1) -> y = dy1(x-x1) + y1 j = lambda x: dy1 * (x-x1) + y1 plt.plot(x,j(x)) plt.xlabel(\'x\') plt.ylabel(\'y\') plt.grid() plt.show() 偏微分 2変数以上の関数の数値微分は以下の通り。片方を止める。 数値微分の方法は上記と同じものを使った。 begin{eqnarray} frac{partial f(x_0,x_1)}{partial x_0} &=& lim_{hrightarrow infty} frac{f(x_0 +h,x_1)-f(x_0-h,x_1)}{2h} \\ frac{partial f(x_0,x_1)}{partial x_1} &=& lim_{hrightarrow infty} frac{f(x_0,x_1+h)-f(x_0,x_1-h)}{2h} end{eqnarray} ((x_0,x_1)=(1,1))における(x_0)に対する偏微分(frac{partial f(x_0,x_1)}{x_0})、(x_1)に対する偏微分(frac{partial f(x_0,x_1)}{x_1})を求めてみる。 ちゃんと(frac{partial f(x_0,1.0)}{x_0}=2.00..)、(frac{partial f(1.0,x_1)}{x_1}=2.00..)になった。 import matplotlib.pyplot as plt import numpy as np from mpl_toolkits.mplot3d import Axes3D def f(x): return x[0]**2 + x[1]**2 X = np.meshgrid(np.arange(-5., 5., 0.2),np.arange(-5., 5., 0.2)) Z = f(X) fig = plt.figure(figsize=(6, 6)) axes = fig.add_subplot(111, projection=\'3d\') axes.plot_surface(X[0],X[1], Z) f0 = lambda x: x**2 + 1.0**2 f1 = lambda x: 1.0**2 + x**2 df0 = numerical_diff(f0, 1.0) df1 = numerical_diff(f1, 1.0) print(df0) # 2.0000000000000018 print(df1) # 2.0000000000000018 plt.show()

台風で自宅に篭れるから勉強時間をとれるな..、と見積もってたのだけれども、 近所の多摩川がマジで溢れそうでそれどころではなく...。 時間が空いてしまったがゼロから作るDeepLearningを読んで実際に実装する作業を再開する。 今後、パラメータを更新していくのだが、どういう方針でパラメータを更新するか決めておく必要がある。 教師ありデータを使った学習を扱っている訳で、訓練データと対応する教師データが与えられている前提。 何かの学習をした結果のモデルの出力と教師データの差を「誤差」として、「誤差」が小さくなるように パラメータを決めていこうという方針。 例えば手書き文字認識で言うところの「認識精度」を指標に使ってしまうと、 モデルの出力が微小に変化したところで「認識精度」は微小に変化しない状況が発生する。 「認識精度」が変化するときは一気に変化する。これではパラメータをどの方向にずらして良いかわからない。 ※SVMの解説で非線形分離を行う決定境界を0/1損失で決めることの問題点に通じる。 非線形分離を行う決定境界も損失関数により微小な変化に追従して決めていく。 2乗和誤差,クロスエントロピー誤差 ということで誤差関数を導入する。 (y_k)はモデルの出力で、最終段でSoftMax関数を通してある。 まずは2乗和誤差。まぁ簡単で、正解と出力の差を2乗した値を足す。 [mathjax] begin{eqnarray} E = frac{1}{2} sum_{k=1}^N (y_k - t_k)^2 end{eqnarray} 次にクロスエントロピー誤差。 begin{eqnarray} E = - sum_{k=1}^N t_k log y_k end{eqnarray} どっちでも良いんじゃ..、と思う訳だけれども、非線形分離問題で決定境界を決めるときに、 正解をより正解として、誤りをより誤りとして表現できる誤差がより優秀なので、 クロスエントロピー誤差の方が適切ではある。 クロスエントロピー誤差の方は(t_k)がゼロの項はゼロになるので、(t_k)が1の項だけ計算すれば良い。 つまり、正解が1のケースについてのみ誤差値が発生する。 (-log y_k)は(y_k)がゼロに近いと急激に値が大きくなる。 これにより、(t_k=1)なのにゼロに近い(y_k)が出力されたときに大きなペナルティを与えられる。 クロスエントロピー誤差関数の実装 バッチ(並列実行)対応のクロスエントロピー誤差関数。 def cross_entropy_error(y, t): if y.ndim == 1: t = t.reshape(1, t.size) y = y.reshape(1,y.size) batch_size = y.shape[0] delta = 10e-7 return -np.sum( t * np.log( y + delta)) / batch_size 1次元のデータを与えてみる。 t1 = [0, 0, 1, 0, 0, 0, 0, 0, 0, 0] y1 = [0.1, 0.05, 0.6, 0.0, 0.05, 0.1, 0.0, 0.1, 0.0, 0.0] cross_entropy_error(np.array(y1), np.array(t1)) # 0.5108239571007129 2次元のデータを与えてみる。 t2 = [ [0, 0, 1, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0, 0]] y2 = [[0.1, 0.05, 0.6, 0.0, 0.05, 0.1, 0.0, 0.1, 0.0, 0.0], [0.1, 0.05, 0.6, 0.0, 0.05, 0.1, 0.0, 0.1, 0.0, 0.0]] cross_entropy_error(np.array(y2), np.array(t2)) # 7.163167257532494 バッチ対応が簡単に書けるところがかなり美しい。

MNISTから手書き文字認識用のデータセットをロードするクラスを作ってみた。 ロードしたデータセットをpickleでシリアライズ、デシリアライズする機能付き。 後々改造する予定でここに貼ったのはメンテしない予定。 たったこれだけ書くのに40分も要してしまった...。 Pythonの冗長感が半端ないけども慣れるしかない。 マジックコードだらけだけども、MNISTの手書き文字認識用データ取得専用だから仕方ない。 get_image()関数により訓練データ、訓練ラベル、テストデータ、テストラベルを取れる。 データは1次元になって入る。つまり1ファイルごとに28*28のサイズがある。 import os import urllib.request import numpy as np import gzip import pickle class MnistLoader: def __init__(self, mnistdir): self.url_base = \'http://yann.lecun.com/exdb/mnist/\' self.dataset_dir = mnistdir self.save_path = self.dataset_dir + \'mnist.pkl\' self.key_file = { \'train_img\':\'train-images-idx3-ubyte.gz\', \'train_label\':\'train-labels-idx1-ubyte.gz\', \'test_img\':\'t10k-images-idx3-ubyte.gz\', \'test_label\':\'t10k-labels-idx1-ubyte.gz\' } self.dataset = {} if os.path.exists(self.save_path): with open(self.save_path,\'rb\') as f: self.dataset = pickle.load(f) else: self.__load_mnist() self.dataset[\'train_img\'] = self.__load_img(self.key_file[\'train_img\']) self.dataset[\'train_label\'] = self.__load_label(self.key_file[\'train_label\']) self.dataset[\'test_img\'] = self.__load_img(self.key_file[\'test_img\']) self.dataset[\'test_label\'] = self.__load_label(self.key_file[\'test_label\']) with open(self.save_path, \'wb\') as f: pickle.dump(dataset, f, -1) def __load_mnist(self): \'\'\' load mnist data and store to file \'\'\' for v in self.key_file.values(): file_path = self.dataset_dir + v urllib.request.urlretrieve(url_base + v, self.dataset_dir + v) def __load_img(self, file_name): file_path = self.dataset_dir + file_name with gzip.open(file_path, \'rb\') as f: data = np.frombuffer(f.read(), np.uint8, offset=16) data = data.reshape(-1, 784) return data def __load_label(self,file_name): file_path = dataset_dir + \'/\' + file_name with gzip.open(file_path, \'rb\') as f: labels = np.frombuffer(f.read(), np.uint8, offset=8) return labels def get_image(self): return (self.dataset[\'train_img\'] , self.dataset[\'train_label\']), (self.dataset[\'test_img\'], self.dataset[\'test_label\']) mnist_loader = MnistLoader(mnistdir=\'/Users/ikuty/Documents/mnist/\') PILを使って訓練データの1枚目を表示するテスト。 import sys, os sys.path.append(os.pardir) import numpy as np from PIL import Image def img_show(img): pil_img = Image.fromarray(np.uint8(img)) pil_img.show() (x_train, t_train), (x_test, t_test) = mnist_loader.get_image() img = x_train[0] label = t_train[0] bimg = img.reshape(28, 28) img_show(bimg) こんなのが出る。(28x28しかなくて小さすぎなので256x256に引き伸ばして表示。)