カテゴリー: 教養

[mathjax] Record Linkage、Entity Recognitionなど、いわゆる「名寄せ」問題には、割とさまざまな解決策が 提案されている。その1つに確率論的な数学的背景を持つFellegi-Sunterモデルを実装したSplinkがある。 数学的背景について咀嚼して説明してくれているサイトはあまりなく、Splinkの公式ドキュメントが唯一 に近いソースだろうと思う。 今回は5万人の人物データを使って重複を除去するSplinkの主要な機能のサンプルを追ってみた。 数式を読んだところで実りは薄い、ということもあり、なるべく数式を読まずに読み進めたが、 モデルのパラメタ、予測結果を解釈する際にどうしても数式を読まないとわからない箇所があった。 特にベイズ統計について理解がかなり浅くて、正直なところ理解しきれてはいない。 Splinkの公式ドキュメントを読み進めていった際に気付いたメモを書いてみることにした。 [arst_toc tag=\"h4\"] 環境構築・準備 環境構築の実装部分は主な学習要素ではないのでClaude Codeでサクッと作ってみた。 まずはPlan modeで計画を立ててもらう。これぐらいだとほぼ言い値で使える。 DuckDB、Splinkを使う最小構成の環境を作ります。環境はDocker composeで構築します。計画を立ててください。 ⏺ DuckDB、Splinkを使うDocker Compose環境の計画を完了しました。 計画概要 1. 設計文書: .tmp/design.mdで要件、アーキテクチャ、技術仕様を定義 2. タスク管理: .tmp/task.mdで詳細な実装手順を3フェーズに分けて整理 主要コンポーネント - Jupyter Lab環境: DuckDB、Splink、必要なライブラリがプリインストール - データボリューム: サンプルデータ、notebook、出力用のディレクトリ - 最小構成: 開発・学習用途に特化した軽量環境 次のステップ 実装を開始する場合は、Phase 1から順番に進めます： 1. Docker設定ファイル作成 2. ディレクトリ構造構築 3. サンプル実装 4. テストと検証 計画の詳細は.tmp/design.mdと.tmp/task.mdをご確認ください。実装を開始しますか？ 実際にファイルを作ってもらうと以下のような構成となる。 得られる中間成果物である docker-compose.yml、Dockerfile は本記事では省略する。 . ├── data ├── docker-compose.yml ├── Dockerfile ├── notebooks ├── output └── requirements.txt 普通に http://localhost:8888 で JupyterLab が開く。 使用するサンプルデータ 5万人の人物データを使って名寄せを行うサンプル。おそらくSplinkの用途として最初に思いつくやつ。 Splinkにデータをロードする前に必要なデータクリーニング手順について説明がある。 公式によると、まずは行に一意のIDを割り当てる必要がある。 データセット内で一意となるIDであって、重複除去した後のエンティティを識別するIDのことではない。 [clink implicit=\"false\" url=\"https://moj-analytical-services.github.io/splink/demos/tutorials/01_Prerequisites.html\" imgurl=\"https://user-images.githubusercontent.com/7570107/85285114-3969ac00-b488-11ea-88ff-5fca1b34af1f.png\" title=\"Data Prerequisites\" excerpt=\"Splink では、リンクする前にデータをクリーンアップし、行に一意の ID を割り当てる必要があります。このセクションでは、Splink にデータをロードする前に必要な追加のデータクリーニング手順について説明します。\"] 使用するサンプルデータは以下の通り。 from splink import splink_datasets df = splink_datasets.historical_50k df.head() データの分布を可視化 splink.exploratoryのprofile_columnsを使って分布を可視化してみる。 from splink import DuckDBAPI from splink.exploratory import profile_columns db_api = DuckDBAPI() profile_columns(df, db_api, column_expressions=[\"first_name\", \"substr(surname,1,2)\"]) 同じ姓・名の人が大量にいることがわかる。 ブロッキングとブロッキングルールの評価 テーブル内のレコードが他のレコードと「同一かどうか」を調べるためには、 基本的には、他のすべてのレコードとの何らかの比較操作を行うこととなる。 全てのレコードについて全てのカラム同士を比較したいのなら、 対象のテーブルをCROSS JOINした結果、各カラム同士を比較することとなる。 SELECT ... FROM input_tables as l CROSS JOIN input_tables as r あるカラムが条件に合わなければ、もうその先は見ても意味がない、 というケースは多い。例えば、まず first_name 、surname が同じでなければ、 その先の比較を行わない、というのはあり得る。 SELECT ... FROM input_tables as l INNER JOIN input_tables as r ON l.first_name = r.first_name AND l.surname = r.surname このような考え方をブロッキング、ON句の条件をブロッキングルールと言う。 ただ、これだと性と名が完全一致していないレコードが残らない。 そこで、ブロッキングルールを複数定義し、いずれかが真であれば残すことができる。 ここでポイントなのが、ブロッキングルールを複数定義したとき、 それぞれのブロッキングルールで重複して選ばれるレコードが発生した場合、 Splinkが自動的に排除してくれる。 このため、ブロッキングルールを重ねがけすると、最終的に残るレコード数は一致する。 ただ、順番により、同じルールで残るレコード数は変化する。 逆に言うと、ブロッキングルールを足すことで、重複除去後のOR条件が増えていく。 積算グラフにして、ブロッキングルールとその順番の効果を見ることができる。 from splink import DuckDBAPI, block_on from splink.blocking_analysis import ( cumulative_comparisons_to_be_scored_from_blocking_rules_chart, ) blocking_rules = [ block_on(\"substr(first_name,1,3)\", \"substr(surname,1,4)\"), block_on(\"surname\", \"dob\"), block_on(\"first_name\", \"dob\"), block_on(\"postcode_fake\", \"first_name\"), block_on(\"postcode_fake\", \"surname\"), block_on(\"dob\", \"birth_place\"), block_on(\"substr(postcode_fake,1,3)\", \"dob\"), block_on(\"substr(postcode_fake,1,3)\", \"first_name\"), block_on(\"substr(postcode_fake,1,3)\", \"surname\"), block_on(\"substr(first_name,1,2)\", \"substr(surname,1,2)\", \"substr(dob,1,4)\"), ] db_api = DuckDBAPI() cumulative_comparisons_to_be_scored_from_blocking_rules_chart( table_or_tables=df, blocking_rules=blocking_rules, db_api=db_api, link_type=\"dedupe_only\", ) 積算グラフは以下の通り。積み上がっている数値は「比較の数」。 要は、論理和で条件を足していって、次第に緩和されている様子がわかる。 DuckDBでは比較の数を2,000万件以内、Athena,Sparkでは1億件以内を目安にせよとのこと。 比較の定義 Splinkは Fellegi-Sunter model モデル (というかフレームワーク) に基づいている。 https://moj-analytical-services.github.io/splink/topic_guides/theory/fellegi_sunter.html 各カラムの同士をカラムの特性に応じた距離を使って比較し、重みを計算していく。 各カラムの比較に使うためのメソッドが予め用意されているので、特性に応じて選んでいく。 以下では、first_name, sur_name に ForenameSurnameComparison が使われている。 dobにDateOfBirthComparison、birth_place、ocupationにExactMatchが使われている。 import splink.comparison_library as cl from splink import Linker, SettingsCreator settings = SettingsCreator( link_type=\"dedupe_only\", blocking_rules_to_generate_predictions=blocking_rules, comparisons=[ cl.ForenameSurnameComparison( \"first_name\", \"surname\", forename_surname_concat_col_name=\"first_name_surname_concat\", ), cl.DateOfBirthComparison( \"dob\", input_is_string=True ), cl.PostcodeComparison(\"postcode_fake\"), cl.ExactMatch(\"birth_place\").configure(term_frequency_adjustments=True), cl.ExactMatch(\"occupation\").configure(term_frequency_adjustments=True), ], retain_intermediate_calculation_columns=True, ) # Needed to apply term frequencies to first+surname comparison df[\"first_name_surname_concat\"] = df[\"first_name\"] + \" \" + df[\"surname\"] linker = Linker(df, settings, db_api=db_api) ComparisonとComparison Level ここでSplinkツール内の比較の概念の説明。以下の通り概念に名前がついている。 Data Linking Model ├─-- Comparison: Date of birth │ ├─-- ComparisonLevel: Exact match │ ├─-- ComparisonLevel: One character difference │ ├─-- ComparisonLevel: All other ├─-- Comparison: First name │ ├─-- ComparisonLevel: Exact match on first_name │ ├─-- ComparisonLevel: first_names have JaroWinklerSimilarity > 0.95 │ ├─-- ComparisonLevel: first_names have JaroWinklerSimilarity > 0.8 │ ├─-- ComparisonLevel: All other モデルのパラメタ推定 モデルの実行に必要なパラメタは以下の3つ。Splinkを用いてパラメタを得る。 ちなみに u は \"\'U\'nmatch\"、m は \"\'M\'atch\"。背後の数式の説明で現れる。 No パラメタ 説明 1 無作為に選んだレコードが一致する確率 入力データからランダムに取得した2つのレコードが一致する確率 (通常は非常に小さい数値) 2 u値(u確率) 実際には一致しないレコードの中で各 ComparisonLevel に該当するレコードの割合。具体的には、レコード同士が同じエンティティを表すにも関わらず値が異なる確率。例えば、同じ人なのにレコードによって生年月日が違う確率。これは端的には「データ品質」を表す。名前であればタイプミス、別名、ニックネーム、ミドルネーム、結婚後の姓など。 3 m値(m確率) 実際に一致するレコードの中で各 ComparisonLevel に該当するレコードの割合。具体的には、レコード同士が異なるエンティティを表すにも関わらず値が同じである確率。例えば別人なのにレコードによって性・名が同じ確率 (同姓同名)。性別は男か女かしかないので別人でも50%の確率で一致してしまう。 無作為に選んだレコードが一致する確率 入力データからランダムに抽出した2つのレコードが一致する確率を求める。 値は0.000136。すべての可能なレコードのペア比較のうち7,362.31組に1組が一致すると予想される。 合計1,279,041,753組の比較が可能なため、一致するペアは合計で約173,728.33組になると予想される、 とのこと。 linker.training.estimate_probability_two_random_records_match( [ block_on(\"first_name\", \"surname\", \"dob\"), block_on(\"substr(first_name,1,2)\", \"surname\", \"substr(postcode_fake,1,2)\"), block_on(\"dob\", \"postcode_fake\"), ], recall=0.6, ) > Probability two random records match is estimated to be 0.000136. > This means that amongst all possible pairwise record comparisons, > one in 7,362.31 are expected to match. > With 1,279,041,753 total possible comparisons, > we expect a total of around 173,728.33 matching pairs u確率の推定 実際には一致しないレコードの中でComparisonの評価結果がPositiveである確率。 基本、無作為に抽出したレコードは一致しないため、「無作為に抽出したレコード」を 「実際には一致しないレコード」として扱える、という点がミソ。 probability_two_random_records_match によって得られた値を使ってu確率を求める。 estimate_u_using_random_sampling によって、ラベルなし、つまり教師なしでu確率を得られる。 レコードのペアをランダムでサンプルして上で定義したComparisonを評価する。 ランダムサンプルなので大量の不一致が発生するが、各Comparisonにおける不一致の分布を得ている。 これは、例えば性別について、50%が一致、50%が不一致である、という分布を得ている。 一方、例えば生年月日について、一致する確率は 1%、1 文字の違いがある確率は 3%、 その他はすべて 96% の確率で発生する、という分布を得ている。 linker.training.estimate_u_using_random_sampling(max_pairs=5e6) > ----- Estimating u probabilities using random sampling ----- > > Estimated u probabilities using random sampling > > Your model is not yet fully trained. Missing estimates for: > - first_name_surname (no m values are trained). > - dob (no m values are trained). > - postcode_fake (no m values are trained). > - birth_place (no m values are trained). > - occupation (no m values are trained). m確率の推定 「実際に一致するレコード」の中で、Comparisonの評価がNegativeになる確率。 そもそも、このモデルを使って名寄せ、つまり「一致するレコード」を見つけたいのだから、 モデルを作るために「実際に一致するレコード」を計算しなければならないのは矛盾では..となる。 無作為抽出結果から求められるu確率とは異なり、m確率を求めるのは難しい。 もしラベル付けされた「一致するレコード」、つまり教師データセットがあるのであれば、 そのデータセットを使ってm確率を求められる。 例えば、日本人全員にマイナンバーが振られて、全てのレコードにマイナンバーが振られている、 というアナザーワールドがあるのであれば、マイナンバーを使ってm確率を推定する。(どういう状況??) ラベル付けされたデータがないのであれば、EMアルゴリズムでm確率を求めることになっている。 EMアルゴリズムは反復的な手法で、メモリや収束速度の点でペア数を減らす必要があり、 例ではブロッキングルールを設定している。 以下のケースでは、first_nameとsurnameをブロッキングルールとしている。 つまり、first_name, surnameが完全に一致するレコードについてペア比較を行う。 この仮定を設定したため、first_name, surname (first_name_surname) のパラメタを推定できない。 training_blocking_rule = block_on(\"first_name\", \"surname\") training_session_names = ( linker.training.estimate_parameters_using_expectation_maximisation( training_blocking_rule, estimate_without_term_frequencies=True ) ) > ----- Starting EM training session ----- > > Estimating the m probabilities of the model by blocking on: > (l.\"first_name\" = r.\"first_name\") AND (l.\"surname\" = r.\"surname\") > > Parameter estimates will be made for the following comparison(s): > - dob > - postcode_fake > - birth_place > - occupation > > Parameter estimates cannot be made for the following comparison(s) since they are used in the blocking rules: > - first_name_surname > > Iteration 1: Largest change in params was 0.248 in probability_two_random_records_match > Iteration 2: Largest change in params was 0.0929 in probability_two_random_records_match > Iteration 3: Largest change in params was -0.0237 in the m_probability of birth_place, level `Exact match on > birth_place` > Iteration 4: Largest change in params was 0.00961 in the m_probability of birth_place, level `All other >comparisons` > Iteration 5: Largest change in params was -0.00457 in the m_probability of birth_place, level `Exact match on birth_place` > Iteration 6: Largest change in params was -0.00256 in the m_probability of birth_place, level `Exact match on birth_place` > Iteration 7: Largest change in params was 0.00171 in the m_probability of dob, level `Abs date difference Iteration 8: Largest change in params was 0.00115 in the m_probability of dob, level `Abs date difference Iteration 9: Largest change in params was 0.000759 in the m_probability of dob, level `Abs date difference Iteration 10: Largest change in params was 0.000498 in the m_probability of dob, level `Abs date difference Iteration 11: Largest change in params was 0.000326 in the m_probability of dob, level `Abs date difference Iteration 12: Largest change in params was 0.000213 in the m_probability of dob, level `Abs date difference Iteration 13: Largest change in params was 0.000139 in the m_probability of dob, level `Abs date difference Iteration 14: Largest change in params was 9.04e-05 in the m_probability of dob, level `Abs date difference <= 10 year` 同様にdobをブロッキングルールに設定して実行すると、dob以外の列についてパラメタを推定できる。 training_blocking_rule = block_on(\"dob\") training_session_dob = ( linker.training.estimate_parameters_using_expectation_maximisation( training_blocking_rule, estimate_without_term_frequencies=True ) ) > ----- Starting EM training session ----- > > Estimating the m probabilities of the model by blocking on: > l.\"dob\" = r.\"dob\" > > Parameter estimates will be made for the following comparison(s): > - first_name_surname > - postcode_fake > - birth_place > - occupation > > Parameter estimates cannot be made for the following comparison(s) since they are used in the blocking rules: > - dob > > Iteration 1: Largest change in params was -0.474 in the m_probability of first_name_surname, level `Exact match on first_name_surname_concat` > Iteration 2: Largest change in params was 0.052 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 3: Largest change in params was 0.0174 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 4: Largest change in params was 0.00532 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 5: Largest change in params was 0.00165 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 6: Largest change in params was 0.00052 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 7: Largest change in params was 0.000165 in the m_probability of first_name_surname, level `All other comparisons` > Iteration 8: Largest change in params was 5.29e-05 in the m_probability of first_name_surname, level `All other comparisons` > > EM converged after 8 iterations > > Your model is not yet fully trained. Missing estimates for: > - first_name_surname (some u values are not trained). モデルパラメタの可視化 m確率、u確率の可視化。 マッチウェイトの可視化。マッチウェイトは (log_2 (m / u))で計算される。 linker.visualisations.match_weights_chart() モデルの保存と読み込み 以下でモデルを保存できる。 settings = linker.misc.save_model_to_json( \"./saved_model_from_demo.json\", overwrite=True ) 以下で保存したモデルを読み込める。 import json settings = json.load( open(\'./saved_model_from_demo.json\', \'r\') ) リンクするのに十分な情報が含まれていないレコード 「John Smith」のみを含み、他のすべてのフィールドがnullであるレコードは、 他のレコードにリンクされている可能性もあるが、潜在的なリンクを明確にするには十分な情報がない。 以下により可視化できる。 linker.evaluation.unlinkables_chart() 横軸は「マッチウェイトの閾値」。縦軸は「リンクするのに十分な情報が含まれないレコード」の割合。 マッチウェイト閾値=6.11ぐらいのところを見ると、入力データセットのレコードの約1.3%が リンクできないことが示唆される。 訓練済みモデルを使って未知データのマッチウェイトを予測 上で構築した推定モデルを使用し、どのペア比較が一致するかを予測する。 内部的には以下を行うとのこと。 blocking_rules_to_generate_predictionsの少なくとも1つと一致するペア比較を生成 Comparisonで指定されたルールを使用して、入力データの類似性を評価 推定された一致重みを使用し、要求に応じて用語頻度調整を適用して、最終的な一致重みと一致確率スコアを生成 df_predictions = linker.inference.predict(threshold_match_probability=0.2) df_predictions.as_pandas_dataframe(limit=1) > Blocking time: 0.88 seconds > Predict time: 1.91 seconds > > -- WARNING -- > You have called predict(), but there are some parameter estimates which have neither been estimated or > specified in your settings dictionary. To produce predictions the following untrained trained parameters will > use default values. > Comparison: \'first_name_surname\': > u values not fully trained records_to_plot = df_e.to_dict(orient=\"records\") linker.visualisations.waterfall_chart(records_to_plot, filter_nulls=False) predictしたマッチウェイトの可視化、数式との照合 predictしたマッチウェイトは、ウォーターフォール図で可視化できる。 マッチウェイトは、モデル内の各特徴量によって一致の証拠がどの程度提供されるかを示す中心的な指標。 (lambda)は無作為抽出した2つのレコードが一致する確率。(K=m/u)はベイズ因子。 begin{align} M &= log_2 ( frac{lambda}{1-lambda} ) + log_2 K \\ &= log_2 ( frac{lambda}{1-lambda} ) + log_2 m - log_2 u end{align} 異なる列の比較が互いに独立しているという仮定を置いていて、 2つのレコードのベイズ係数が各列比較のベイズ係数の積として扱う。 begin{eqnarray} K_{feature} = K_{first_name_surname} + K_{dob} + K_{postcode_fake} + K_{birth_place} + K_{occupation} + cdots end{eqnarray} マッチウェイトは以下の和。 begin{eqnarray} M_{observe} = M_{prior} + M_{feature} end{eqnarray} ここで begin{align} M_{prior} &= log_2 (frac{lambda}{1-lambda}) \\ M_{feature} &= M_{first_name_surname} + M_{dob} + M_{postcode_fake} + M_{birth_place} + M_{occupation} + cdots end{align} 以下のように書き換える。 begin{align} M_{observe} &= log_2 (frac{lambda}{1-lambda}) + sum_i^{feature} log_2 (frac{m_i}{u_i}) \\ &= log_2 (frac{lambda}{1-lambda}) + log_2 (prod_i^{feature} (frac{m_i}{u_i}) ) end{align} ウォーターフォール図の一番左、赤いバーは(M_{prior} = log_2 (frac{lambda}{1-lambda}))。 特徴に関する追加の知識が考慮されていない場合のマッチウェイト。 横に並んでいる薄い緑のバーは (M_{first_name_surname} + M_{dob} + M_{postcode_fake} + M_{birth_place} + M_{occupation} + cdots)。 各特徴量のマッチウェイト。 一番右の濃い緑のバーは2つのレコードの合計マッチウェイト。 begin{align} M_{feature} &= M_{first_name_surname} + M_{dob} + M_{postcode_fake} + M_{birth_place} + M_{occupation} + cdots \\ &= 8.50w end{align} まとめ 長くなったのでいったん終了。この記事では教師なし確率的名寄せパッケージSplinkを使用してモデルを作ってみた。 次の記事では、作ったモデルを使用して実際に名寄せをしてみる。 途中、DuckDBが楽しいことに気づいたので、DuckDBだけで何個か記事にしてみようと思う。

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 バッチ対応が簡単に書けるところがかなり美しい。