カテゴリー: 機械学習

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] 参考書籍を読んでいく. 今回はロジスティック回帰. 未知のサンプルのクラスの所属確率を見積もることができる, という重要なことが書かれている. また重みの学習は勾配降下法のフレームワークの中で統一されていて一連の理解が進む. 一度は理解しておいた方が良い内容だと思った. パーセプトロンと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}} ]

[mathjax] 球面集中現象を理解するために記憶にないことが多すぎるので 理解に必要そうなことを少しずつ復習していきます。 まず極座標の直行座標変換。極座標を使って球の体積を求めてみます。 次の記事で球の体積を求める時に使う微小体積がヤコビアンになることを確認します。 極座標の直行座標変換 2次元のとき、極座標((r),(phi)) を直行座標で表すと、 begin{eqnarray} x &=& r cos(phi) \\ y &=& r sin(phi) end{eqnarray} 3次元のとき、極座標((r),(phi_1),(phi_2))を直行座標で表すと、 begin{eqnarray} x &=& r cos(phi_1) \\ y &=& r sin(phi_1) cos(phi_2) \\ z &=& r sin(phi_1) sin(phi_2) end{eqnarray} ここまではイメージできるのだけれども、(N)に飛ばしたときにどうなるのか。 N次元のとき、極座標((r),(phi_1),(phi_2),(cdots),(phi_{N-1}))を直行座標で表すと 以下のようになるらしい。 begin{eqnarray} x_1 &=& r cos(phi_1) \\ x_2 &=& r sin(phi_1) cos(phi_2) \\ x_3 &=& r sin(phi_1) sin(phi_2) cos(phi_3) \\ vdots \\ x_{N-1} &=& r sin (phi_1) cdots sin(phi_{N-2}) cos(phi_{N-1}) \\ x_{N} &=& r sin (phi_1) cdots sin(phi_{N-2}) sin(phi_{N-1}) end{eqnarray} この後使わないけれど、逆変換は以下の通りになるらしい。 begin{eqnarray} r &=& sqrt{x_1^2 + cdots + x_{N-1}^2 + x_N^2} \\ phi_1 &=& arccosfrac{x_1}{sqrt{x_1^2 + cdots + x_{N-1}^2 + x_N^2}} \\ phi_2 &=& arccosfrac{x_2}{sqrt{x_2^2 + cdots + x_{N-1}^2 + x_N^2}} \\ vdots \\ phi_{N-2} &=& arccos frac{x_{N-2}}{sqrt{x_{N-2}^2 + x_{N-1}^2 + x_N^2}} \\ phi_{N-1} &=& begin{cases} arccos frac{x_{N-1}}{sqrt{x_{N-1}^2+X_N^2}} & X_N ge 0 \\ -arccos frac{x_{N-1}}{sqrt{x_{N-1}^2+X_N^2}} & X_N lt 0 end{cases} end{eqnarray} 3次元極座標系における体積計算の例 極座標((r),(phi_1),(phi_2))において(r)を微小量(Delta r)だけ増やす。 (2pi r)は半径を(r)とする円の直径だということを考慮すると、 (phi_1)から微小角(Delta phi_1)増やしたときの(r Delta phi_1)は円弧の一部。 (Delta phi_1)を限りなく小さくすると(r Delta phi_1)は直線と考えられる。 (2pi rsin(phi_1) )が半径を(rsin(phi_1))とする円の直径と考えると、 (phi_2)から微小角(Delta phi_2)増やしたときの(r sin(phi_1) Delta phi_2)は円弧の一部。 (Delta phi_2)を限りなく小さくすると(r sin(phi_1) Delta phi_2)は直線と考えられる。 (Delta r, Delta phi_1, Delta phi_2)が限りなくゼロに近くときに 上記のそれぞれの直線を辺とする直方体ができる。その体積を(dV)とすると、 begin{eqnarray} dV &=& Delta r cdot r Delta phi_1 cdot r sin(phi_1) Delta phi_2 \\ &=& r^2 sin(phi_1) Delta r Delta phi_1 Delta phi_2 end{eqnarray} 直交座標(x,y,z)の位置を(U(x,y,z))とする。 直交座標系における球の体積はざっくり以下で表せる。 begin{eqnarray} int U(x,y,z) dV end{eqnarray} 直行座標の極座標変換は前述の通り以下。この座標を(U(r,phi_1,phi_2))とする。 begin{eqnarray} x &=& r cos(phi_1) \\ y &=& r sin(phi_1) cos(phi_2) \\ z &=& r sin(phi_1) sin(phi_2) end{eqnarray} 極座標系において以下だから、 begin{eqnarray} 0 le r le 1 \\ 0 le phi_1 le pi\\ 0 le phi_2 2pi end{eqnarray} それぞれについて積分すると体積になる。 begin{eqnarray} int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{1} U(r,phi_1,phi_2)r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 end{eqnarray} つまり以下。 begin{eqnarray} int U(x,y,z) dV = int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{1} U(r,phi_1,phi_2)r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 end{eqnarray} ここで( U(r,phi_1,phi_2) =1 )として積分範囲を(r)とすると球の体積の公式になる。 begin{eqnarray} int U(x,y,z) dV &=& int_{0}^{2pi} Biggl( int_{0}^{pi} Bigl( int_{0}^{r} r^2 sin(phi_1) dr Bigr) d phi_1 Biggr) dphi_2 \\ &=& int_{0}^{2pi} int_{0}^{pi} Bigl( int_{0}^{r} frac{sin(phi_1)}{3} bigl[ r^3 bigr]_{0}^{r} Bigr) dphi_1 dphi_2 \\ &=& int_{0}^{2pi} int_{0}^{pi} frac{r^3 sin(phi_1)}{3} dphi_1 dphi_2 \\ &=& int_{0}^{2pi} -frac{r^3}{3} bigl[ cos(phi_1) bigr]_{0}^{pi} dphi_2 \\ &=& int_{0}^{2pi} frac{2r^3}{3} dphi_2 \\ &=& frac{2r^3}{3} bigl[ phi_2bigr]_{0}^{2pi} \\ &=& frac{4pi r^3}{3} end{eqnarray} 次回、微小体積とヤコビアンが等しい理由を書いてみます。

[mathjax] sihoutteって何て読むのか...と思うけども\"シルエット\"だそう。フランス語源。 ポートレート写真を背景白、顔を黒に減色した例の\"シルエット\"。\"輪郭\"みたいな。 データ達を複数のクラスタに分割したとして、各々のデータがそのクラスタに存在することの 収まりの良さを表すことができる。 本来(相対的に)他のクラスタに属しているべきデータ達の割合が分かったり、 クラスタ境界で(相対的に)どちらのクラスタに属しても良さそうなデータ達の割合が分かったりする。 sklearnは教師なしクラスタリングまでで、それをグラフ化しないといけないのだけど、 matplotlibに対応するのがなく、\"線を書く\"みたいなコマンドを並べて作っていかないといけない様子。 yellowbrickというパッケージを使うとそれもやってくれる。 Yellowbrick is a suite of visual diagnostic tools called “Visualizers” that extend the Scikit-Learn API to allow human steering of the model selection process. In a nutshell, Yellowbrick combines scikit-learn with matplotlib in the best tradition of the scikit-learn documentation, but to produce visualizations for your models! sklearnとmatplotloibだけで出力するバージョンと、yellowbrickで出力するバージョンの 両方を試してみた(前者はクラスタのラベルを出力できずsihoutteグラフとの対応関係が理解できずに 中途半端で終了)。 sihoutte係数 全データ達は(x_i)。データ(x_i)がクラスタ(A)に属し、最近傍にクラスタ(B)があるという状況。 (A)のクラスタ中心は(x_A)、(B)のクラスタ中心は(x_B)。 \"最近傍のクラスタ中心との距離の平均\"から\"自分のクラスタ中心との距離の平均\"を引いた値。 わざわざ1行に全ての変数が出てくるように書いてみる。 begin{eqnarray} s_i = frac{sum_{i=1}^{n}|x_i-x_B|/n - sum_{i=1}^{n}|x_i-x_A|/n }{max bigl{sum_{i=1}^{n}|x_i-x_A|/n,sum_{i=1}^{n}|x_i-x_B|/n bigr}} end{eqnarray} データが(n)個あるので、sihoutte係数も(n)個できる。 自分が属しているクラスタ(A)よりも、隣のクラスタ(B)の方が居心地が良いと 分子はマイナスになる。 今属しているクラスタでも隣のクラスタでも、どちらもさほど居心地が変わらないとゼロ近辺になる。 大きければ迷いなく今のクラスタで良いことを表せる。 sklearnとmatplotlibだけでsihoutte係数 乱数で作った偽データに6-meansをかけた図。 どうやってもクラスタ中心のラベルが取り出せない!!..(スルー..) import numpy as np import matplotlib.pyplot as plt import seaborn as sns mu = [[0,0], [20,20], [50,50], [40,30], [40,10], [20,40]] sigma = [ [[30,20],[20,50]], [[20,30],[10,20]], [[60,40],[20,20]], [[60,20],[20,60]] ,[[30,10],[10,30]],[[50,20],[20,50]] ] points = 100 clusters = [] for index in range(len(mu)): cluster = np.random.multivariate_normal(mu[index], sigma[index], points) dig = np.full((points,1),index+1, dtype=int) cluster = np.hstack((cluster,dig)) clusters = np.r_[clusters,cluster] if len(clusters) > 0 else cluster plt.scatter(x=clusters[:,0], y=clusters[:,1],c=clusters[:,2]) from sklearn.cluster import KMeans model = KMeans(n_clusters=6, init=\'random\',max_iter=10) y_km = model.fit_predict(clusters[:,:2]) labels = model.labels_ centers = model.cluster_centers_ centers fig = plt.figure() ax1 = fig.add_subplot(1,1,1) ax1.scatter(x=clusters[:,0], y=clusters[:,1],c=labels) ax2 = fig.add_subplot(1,1,1) ax2.scatter(x=centers[:,0], y=centers[:,1], alpha=0.5,s=600,c=\"pink\",linewidth=2,edgecolors=\"red\") sklearnとmatplotlibだけでsihoutte係数のグラフを出してみる。 \"線を書く\"みたいなコマンドを並べて規格通りの図を作るのか...。 from sklearn.metrics import silhouette_samples # cluster数 num_clusters=6 #全データのsilhouette係数を取得 silhouette_vals = silhouette_samples(clusters[:,:2],y_km,metric=\'euclidean\') cluster_labels = np.unique(y_km) min,max = 0,0 yticks = [] for i,c in enumerate(cluster_labels): c_silhouette_vals = silhouette_vals[y_km == c] c_silhouette_vals.sort() max += len(c_silhouette_vals) plt.barh( range(min,max), c_silhouette_vals, height=1.0, edgecolor=\'none\' ) yticks.append((min+max)/2.) min += len(c_silhouette_vals) avg = np.mean(silhouette_vals) plt.axvline(avg, color=\'red\', linestyle=\"--\") plt.yticks(yticks, cluster_labels + 1) plt.ylabel(\'Cluster\') plt.xlabel(\'Silhouette coefficient\') plt.show() 全体的に茶色のクラスタのsilhouette係数が低め。 残念ながら、どのクラスタと対応するのか出力できず何の考察も出来ず...。 yellowbrickでsilhouette係数を出力 無茶苦茶簡単に出せる。 from yellowbrick.cluster import SilhouetteVisualizer sv = SilhouetteVisualizer(model) sv.fit(clusters[:,:2]) sv.poof() 出てきた図。自作したものと全然合っていないように見える.. 自作したものとラベルが合っていないだけだと信じたい。 似た形の塊があるので.. 以上、完全な失敗だけれども一旦終了...

k-means法を実行する際に妥当なkを決めたいという欲求があります。 クラスタ集合の凝集度を定量化することでkと凝集度の関係を得られます。 複数のkについてkと凝集度の関係を取得し、 そこから妥当なkを発見する方法がエルボー法です。 この記事においては以下を書いてみます。 クラスタ集合の凝集度を定量化する方法 kと凝集度の関係を表すグラフの書き方 クラスタ集合の凝集度を定量化する方法 これ以上無いくらいシンプルなやり方です。 データポイントとcentroidの距離の2乗和(SSE)を凝集度として使います。 それを全クラスタに関して計算し足し合わせます。 足し合わせた値を歪み(distortion)と言っているところが多いです。 クラスタ数kが増えるにつれてこの値は小さくなるだろうと予測できます。 横軸にk、縦軸にSSEの和をとるグラフを書いてみると、 kの増加に伴ってSSEの和が小さくなること、減少幅は次第に小さくなりそうです。 (すべてのケースでそうなるのかは不明) 減少幅が緩やかになる地点のkが費用対効果が高いkですし、 ほぼ減少が止まる(サチる)地点のkを採用することは、 それ以上kを増加させても意味がないという事実を根拠として使えそうです。 擬似データへの適用 前回から使っている多次元正規分布に従う擬似クラスタデータに対してk-meansを適用します。 k=1から9まで適用し、k-distortionグラフを書いてみます。(だいぶPythonに慣れてきた..。) import numpy as np import matplotlib.pyplot as plt import seaborn as sns mu = [[0,0], [20,20], [50,50], [40,30], [40,10], [20,40]] sigma = [ [[30,20],[20,50]], [[20,30],[10,20]], [[60,40],[20,20]], [[60,20],[20,60]] ,[[30,10],[10,30]],[[50,20],[20,50]] ] points = 100 clusters = [] for index in range(len(mu)): cluster = np.random.multivariate_normal(mu[index], sigma[index], points) dig = np.full((points,1),index+1, dtype=int) cluster = np.hstack((cluster,dig)) clusters = np.r_[clusters,cluster] if len(clusters) > 0 else cluster plt.scatter(x=clusters[:,0], y=clusters[:,1],c=clusters[:,2]) from sklearn.cluster import KMeans kmeans_model = KMeans(n_clusters=6, init=\'random\',max_iter=10).fit(clusters[:,:2]) labels = kmeans_model.labels_ distortions = [] numClusters = 10 for i in range(1,numClusters): kmeans_model = KMeans(n_clusters=i, init=\'random\',max_iter=10) kmeans_model.fit(clusters[:,:2]) distortions.append(kmeans_model.inertia_) plt.plot(range(1,numClusters),distortions,marker=\'o\') plt.xlabel(\'Number of clusters\') plt.ylabel(\'Distortion\') plt.show() 複数の擬似点を中心に2次元正規分布に従う散らばりをもったデータ達です。 k-distortionグラフです。もの凄い期待通りなグラフができました。 k=5か6あたりでdistortionがサチっています。 6-meansを適用し、centroidを重ねてみました。 散布図を見ての想像になりますが、 確かに4-meansと5-meansでは、5-meansの方が凝集していそうです。 5-meansと6-meansだと、4と5程凝集度合いに変化がなさそうです。

前回適当に作ったデータを使って階層型クラスタリングを試してみる。 階層型クラスタリングの様子をdendrogramというツールを使って可視化できるのと、 クラスタ間距離の違いによってクラスタの作られ方が大きく変化することを dendrogramの違いを見ながら確認してみる。 テストデータを作る データは前回と同様に以下のように作る。 import numpy as np import matplotlib.pyplot as plt import seaborn as sns mu = [[0,0], [20,20], [50,50], [40,30], [40,10], [20,40]] sigma = [ [[30,20],[20,50]], [[20,30],[10,20]], [[60,40],[20,20]], [[60,20],[20,60]] ,[[30,10],[10,30]],[[50,20],[20,50]] ] points = 100 clusteres = [] for index in range(len(mu)): cluster = np.random.multivariate_normal(mu[index], sigma[index], points) dig = np.full((points,1),index+1, dtype=int) cluster = np.hstack((cluster,dig)) clusteres = np.r_[clusteres,cluster] if len(clusteres) > 0 else cluster plt.scatter(x=clusteres[:,0], y=clusteres[:,1],c=clusteres[:,2]) dendrogram 階層型クラスタリングの可視化を行う方法の1つとしてdendrogramがある。 葉(=各データポイント)から始まり、クラスタがまとまっていく様を木構造で表した図。 根は全てのデータポイントを1つのクラスタで表した状態。 左部分木と右部分木に1度別れるところが2つのクラスタに別れた状態。 グラフの縦軸はクラスタ間距離を表していて、dendrogramによって各クラスタの離れ度合いがわかる。 from sklearn.preprocessing import StandardScaler from scipy.cluster.hierarchy import linkage, dendrogram, fcluster SS = StandardScaler() scaled_clusters = SS.fit_transform(clusteres) result = linkage(scaled_clusters, method = \"ward\") dendrogram(result) plt.show() クラスタを作っていく際にどうクラスタ間距離を設定するかで木の構造がだいぶ変わってくる。 空間拡散性 以下はクラスタ間距離として最短距離法を使った場合。 (クラスタ間距離としてそれぞれに含まれるデータポイント間の距離の最小値） 1度併合してクラスタが大きくなると、次に、そのクラスタと隣のクラスタがくっつきやすい。 (データ数が多すぎてよくわからないけれど)下のdendrogramは 右、その１つ左、その1つ左...というように、右から順番に併合する傾向が強いのがわかる。 つまり1度クラスタが作られると、そのクラスタは次のクラスタになりやすい。 1度作られたクラスタは次にクラスタになりにくい性質、 (=あちこちでバランスよくクラスタが作られていって最後にそれぞれがまとまる性質) を空間拡散性というらしい。ほぅ。下の図は空間拡張性が無いことを表す。 クラスタ数を多くしていったとき、空間拡張性がないと\"でかいクラスタ\"と\"小さいクラスタ\"が共存する アンバランスな構造になりそう。 距離の単調性 以下は重心法を使った場合。 重心法というのは、クラスタとクラスタの距離として各クラスタの重心間の距離を使う方法。 まぁ妥当な距離な感じがするが..。 本来は、同時多発的にクラスタが出来上がり、より大きくなるにつれて、 より遠いクラスタと併合するのが\"あるべき姿\"なのだけれども、 1度クラスタがくっつくと、他のクラスタとの距離が以前よりも小さくなってしまう状態が発生する。 距離の単調性がない、というらしい..。むむ. dendrogramの縦方向はクラスタ間距離なので、 距離の単調性があるのであれば、左右の部分木は両方とも下方向に伸びるはず。 無いのであれば、左右のいずれかが上方向に伸びる。 ward法を使った場合 以下はWard法を使った場合。詳細は以下。無茶苦茶面倒臭い。 [clink url=\"https://ikuty.com/2019/07/09/hierarchical_clustering/\"] 空間拡散性、距離の単調性があって(バランス木という言葉は全然違うが..) バランス良く木ができている。 データ数が多すぎて全然見えないけど...

multivariate_normalを使って6個のデータポイント近隣にデータ達を発生させる。 薄くクラスタを見つけられそうだけれど境界は曖昧で大分被っているという状況。 そんなデータ達に6-meansをかけてみたとき、どうクラスタが出来るのかという実験。 被っている部分がどう別のクラスタに入るのか確認するため。 import numpy as np import matplotlib.pyplot as plt import seaborn as sns mu = [[0,0], [20,20], [50,50], [40,30], [40,10], [20,40]] sigma = [ [[30,20],[20,50]], [[20,30],[10,20]], [[60,40],[20,20]], [[60,20],[20,60]] ,[[30,10],[10,30]],[[50,20],[20,50]] ] points = 100 clusteres = [] for index in range(len(mu)): cluster = np.random.multivariate_normal(mu[index], sigma[index], points) dig = np.full((points,1),index+1, dtype=int) cluster = np.hstack((cluster,dig)) clusteres = np.r_[clusteres,cluster] if len(clusteres) > 0 else cluster plt.scatter(x=clusteres[:,0], y=clusteres[:,1],c=clusteres[:,2]) おもむろにsklearnのKMeansを使ってみる。 n_clustersは6、max_iterを10に設定してみた。 クラスタ中心がわかるように重ねてみた。 from sklearn.cluster import KMeans kmeans_model = KMeans(n_clusters=6, init=\'random\',max_iter=10).fit(clusteres[:,:2]) labels = kmeans_model.labels_ centers = kmeans_model.cluster_centers_ centers fig = plt.figure() ax1 = fig.add_subplot(1,1,1) ax1.scatter(x=clusteres[:,0], y=clusteres[:,1],c=labels) ax2 = fig.add_subplot(1,1,1) ax2.scatter(x=centers[:,0], y=centers[:,1], alpha=0.5,s=600,c=\"pink\",linewidth=2,edgecolors=\"red\") 当然のごとく、微妙に被っていた部分は異なるクラスタに分類された。 今回は当たり前なことを確認して終了。 2次元データじゃつまらない...。

クラスタリングについて凄まじくわかりやすい説明を聞いて感動しました。 k-meansを自力で書くことで理解が深まりそうなので、参考にしながらアプトプットしてみます。 前回のエントリで、多次元正規分布に従う標本を作れるようになったので、 もともと存在する中心から適当に散らばるように作ったデータポイントに対して、 どの程度クラスタ分割できるのか確認する例です。 見てわからないと仕方がないので2次元とします。 早速、2つの中心から適当に散らばるようなデータポイントを作ってみます。 cluster1 = np.random.multivariate_normal([-10,-10], [[20,10],[30,20]], 100) cluster2 = np.random.multivariate_normal([10,10], [[50,5],[5,50]], 100) X = np.r_[cluster1, cluster2] plt.scatter(X[:,0], X[:,1]) 合計200個の点に、確率0.5でラベル1を振って(=確率0.5でラベル0)、 それぞれ色を付けて表示してみます。これを初期クラスタとします。 0クラスタ、1クラスタそれぞれ混じっていて、このままでは適切なクラスタ分割ではありません。 y = np.random.binomial(n = 1, p = 0.5, size = 200) plt.scatter(X[:,0], X[:,1],c=y) ラベル0のデータをX1、ラベル1のデータをX2という風に名前を付けて、 pandasデータフレーム化します。 さらに振ったラベル列を付けます。列追加は.assign()で出来るみたい。 import pandas as pd center = pd.DataFrame(X, columns = [\"X1\", \"X2\"]) center.shape # (200, 2) center = center.assign(y=y) center.shape # (200, 3) 0クラスタと1クラスタでまとめます。SQLで言うところのGroupBy句。 メソッド名もまんまgroupbyなのでわかりやすい。 GroupByした後、平均を計算する。0クラスタ、1クラスタ毎のX1,X2の平均が出ているっぽい。 center = center.groupby(y).mean() # X1 X2 y # 0 -1.701594 -1.186065 0 # 1 1.278607 0.954203 1 いよいよ、各データポイントのクラスタを更新して、適切なクラスタに振っていきます。 pandasデータフレームからNumPyのArrayを取るのは結構面倒くさい。 locにより行を取得できる。その際pandasの列名を指定できる。 SQLでSELECT書いているのに近い。 クラスタ0、クラスタ1それぞれのデータポイントについて、 各クラスタにおけるクラスタ中心と距離を求めます。距離は2乗和。 center_cluster_0 = center.loc[0, [\"X1\", \"X2\"]].values center_cluster_1 = center.loc[1, [\"X1\", \"X2\"]].values center0 = center.loc[0,[\"X1\",\"X2\"]].values center1 = center.loc[0,[\"X1\",\"X2\"]].values dist0 = (X[:,0].reshape(200) - center0[0])**2 + (X[:,1].reshape(200) - center0[1])**2 dist1 = (X[:,0].reshape(200) - center1[0])**2 + (X[:,1].reshape(200) - center1[1])**2 クラスタ0のクラスタ中心との距離よりも、クラスタ1のクラスタ中心との距離の方が大きいということは クラスタ0に属するべきです。論理式の結果がTrueかFalseかのいずれかになることを使って 2つのクラスタに分けます。かなり無理やりな気がします。 真偽の扱いですが、Pythonの言語仕様上、0か0に類する値はFalse、それ以外はTrueです。 new_cluster_id = dist0 < dist1 元Webプログラマが普通にforを使って書いてみると以下の通りとなります。 new_cluster_id2 = [] for index in range(len(dist0)): new_cluster_id2.append( 0 if (dist0[index] < dist1[index]) else 1) np.reshape(new_cluster_id2,200) 生成した入力値に新しいクラスタIDを振って色分け表示してみます。 new_cluster_id = dist0 < dist1 plt.scatter(X[:,0], X[:,1], c = new_cluster_id) 綺麗にわかれた気がしますが、大元のクラスタが交差する部分において、 別のクラスタに分類されたデータがあるようです。

[mathjax] そろそろ適当なデータを見つけてきて手法を試すのとは別に、 自力でデータを作って試してみたいと思い、NumPyを使った生成法を調べてみた。 一口に乱数といっても、正規分布に従う標本の生成のこと。 多次元正規分布に従う標本をmultivariate_normalで生成して表示してみる。 1次元正規分布に従う標本 その前に普通の乱数。 平均(mu)、標準偏差(sigma)の正規分布(N(mu,sigma))に従う標本の生成。 numpy.random.normal((mu,sigma,n))。 以下、(mu=50,sigma=10)として1000個の標本を作って、 階級数100のヒストグラムとして表示する。 import numpy as np import matplotlib.pyplot as plt values = np.random.normal(50, 10,1000) plt .hist(values,bins=100) 多次元正規分布に従う標本 多次元の乱数を作るのに必要なのは、 各次元における平均(mu=begin{pmatrix}mu_1 \\ mu_2 end{pmatrix})と、分散共分散行列(sum=begin{pmatrix}sigma_{11} & sigma_{21} \\ sigma_{12} & sigma_{22}end{pmatrix})。 ちなみに、(sum=begin{pmatrix}sigma_{11} & sigma_{21} \\ sigma_{12} & sigma_{22}end{pmatrix} = begin{pmatrix} sigma_1^2 & sigma_{21} \\ sigma_{12} & sigma_2^2 end{pmatrix} = begin{pmatrix} sigma_1^2 & Cov(X_1,X_2) \\ Cov(X_2,X_1) & sigma_2^2 end{pmatrix} )。 対角に分散、それ以外にクロスする共分散が入る。 共分散は昔書いてた。 begin{eqnarray} r_{xy} &=& frac{C_{xy}}{S_x S_y} \\ &=& frac{sum_{i=1}^n(x_i-bar{x})(y_i-bar{y})/n}{sqrt{sum_{i=1}^n{(x_i-bar{x})^2}/n} sqrt{sum_{i=1}^n{(y_i-bar{y})^2}/n}} \\ &=& frac{sum_{i=1}^n(x_i-bar{x})(y_i-bar{y})}{sqrt{sum_{i=1}^n{(x_i-bar{x})^2}} sqrt{sum_{i=1}^n{(y_i-bar{y})^2}}} \\ end{eqnarray} [clink url=\"https://ikuty.com/2018/10/17/covariance_zero/\"] 平均(mu=begin{pmatrix} 0 \\ 0end{pmatrix})、分散共分散行列(sum=begin{pmatrix}30 & 20 \\ 20 & 50end{pmatrix})の2次元正規分布に従う標本を1000個作る。 import numpy as np import matplotlib.pyplot as plt import seaborn as sns mu = [0,0] sigma = [[30,20],[20,50]] values = np.random.multivariate_normal(mu, sigma, 1000) sns.jointplot(x=values[:,0],y=values[:,1]) 散布図とヒストグラムを良い感じに表示してくれるseaborn.jointplotを使って表示してみる。 第1軸、第2軸とも平均は0っぽい。 第1軸は-15から15、第2軸は-25から25くらいが入っていて、対角とあってそう。