「日々の学び」一覧

やりたいことは以下の2つ。 ansibleでaws-cliをインストールする ansibleでインストールしたaws-cliでs3コマンドを打てるようにする なお、相手には既にpipがインストールがしてあるものとします。 ansibleを実行するために最小構成でPythonをインストールしたもののpipは入れていない、 という状況であれば、先にpipをインストールする必要があります。 リージョン、S3のkey,secretは仮に以下とします。 事前にAWSのコンソールで設定,取得してください。 region: ap-northeast-1 s3.key: AHJKOKODAJOIFAJDJDIOA s3.secret: AugioaiARJOIfjop20FJIOADOiFJAODA ファイル達 構成は以下の通りです。(※)のファイルが核心です。 stagingとかになってますが、もちろん成立する範囲で修正してください。 ├──provision.yml (※) ├──ansible.cfg ├──group_vars │ └───staging.yml (※) ├──hosts | └───staging ├──host_vars | └───default.yml └──roles └─awscli (※) ├─templates | └─config.conf.j2 | └─credentials.conf.j2 └─tasks └─main.yml group_vars/staging.ymlに設定を書きます。 user: ubuntu s3: region: ap-northeast-1 # S3.region key: AHJKOKODAJOIFAJDJDIOA # S3.key secret: AugioaiARJOIfjop20FJIOADOiFJAODA # S3.secret roles/awscli/templates/config.conf.j2にaws-cliの設定を書きます。 s3.regionが評価され値が入ります。相手の~/.aws/configに配置します。 [default] output = json region = {{s3.region}} roles/awscli/templates/credentials.conf.j2にs３の設定を書きます。 s3.keyとs3.secretが評価され値が入ります。相手の~/.aws/credentialsに配置します。 [default] aws_access_key_id = {{s3.key}} aws_secret_access_key = {{s3.secret}} rokes/awscli/tasks/main.ymlに状態を定義します。 内容は以下の通りです。 1)　aws-cliがpip installされた状態 2) ~/.aws/以下に設定ファイルがコピーされた状態 --- - name: install aws-cli pip: name: awscli - name: create .aws dir file: dest: /home/{{user}}/.aws state: directory owner: \"{{user}}\" group: \"{{user}}\" - name: copy config template: src: config.conf.j2 dest: /home/{{user}}/.aws/config owner: \"{{user}}\" group: \"{{user}}\" - name: copy credential template: src: credential.conf.j2 dest: /home/{{user}}/.aws/credentials owner: \"{{user}}\" group: \"{{user}}\" ~ ~ Playbook(provision.yml)は以下の通りです。 - hosts: remote_machine remote_user: \"{{ user }}\" gather_facts: \"{{ check_env | default(True) }}\" become: yes become_method: sudo roles: - { role: awscli } 実行結果 Playbookを実行します。 $ ansible-playbook -i hosts/staging provisiong.yml 相手のユーザディレクトリに.awsというディレクトリが作られ、中にファイルが作られます。 ~/ └─.aws ├─config └─credentials 相手側でaws s3 lsコマンドを打って設定しろと言われなければ成功です。 $ aws s3 ls 2019-10-20 11:11:20 hogehoge おわり。

[mathjax] マルコフ連鎖モンテカルロ法。2変量正規分布からGibbs Samplingする方法を考えてみました。 式を流してもよくわからないので、行間ゼロで理解できるまで細切れにして書いてみます。 Gibbs Sampling 無茶苦茶一般的に書かれている書籍だと以下みたいになっている。 ステップごとに(theta=(theta_1,theta_2,cdots,theta_n))の各要素を、 その要素以外の要素の条件付き確率でランダムに発生させる。 begin{eqnarray} theta^0 &=& (theta_1^0,theta_2^0,cdots,theta_n^0) \\ theta_i^{t+1} &sim& p(theta_i|theta_1^{t+1},cdots,theta_{i-1}^{t+1},theta_{i+1}^t,cdots,theta_n^t) end{eqnarray} 2次元空間であれば、ある点を決める際に、 片方の変数(theta_1)を固定して条件付き確率(P(theta_2|theta_1))を最大にするパラメタ(hat{theta_2})を決める。 その時点で((theta_1,hat{theta_2}))が決まる。 次は変数(hat{theta_2})を固定して条件付き確率(P(theta_1|hat{theta_2}))を最大にするパラメタ(hat{theta_1})を決める。 2次元正規分布であれば、片方の確率変数を固定したときの条件付き確率が1変数の正規分布になるから、 これがすごくやりやすい。 Gibbs Samplingは、このように条件付き確率を計算できる必要がある。 2変量正規分布の条件付き確率 2変量正規分布の片方の確率変数を固定すると1変数の正規分布が出てきます。 山の輪切りにして出てきた正規分布の母数を調べたいのですが、導出が無茶苦茶大変そうです。 出てきた正規分布の母数について式変形すると出てくるようです。こちらを参考にさせて頂きました。 この関係式を利用してひたすら式変形していきます。 begin{eqnarray} f(x|y) &=& frac{f(x,y)}{f(x)} end{eqnarray} 今、確率変数(y=y\')を固定し確率変数(x)の正規分布を考えます。 2変量正規分布の分散共分散行列が(sum)であるとします。 begin{eqnarray} sum &=& begin{pmatrix} sigma_{xx} & sigma_{xy} \\ sigma_{yx} & sigma_{yy} end{pmatrix} end{eqnarray} 出てくる正規分布の平均は以下となります。 確率変数が(x)なのに固定した(y)が出てくるのがポイント。 begin{eqnarray} mu\' = bar{x} + frac{sigma_{xy}}{sigma_{yy}} (y\'-bar{y}) end{eqnarray} また、出てくる正規分布の分散は以下となります。 begin{eqnarray} sigma\'^2 &=& sigma_{xx} - frac{sigma_{xy}}{sigma_{yy}} sigma_{yx} end{eqnarray} Python実装例 実際にPythonコードにしてみた図です。2変量正規分布の条件付き確率さえわかってしまえば あとはコードに落とすだけです。 import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats n_dim = 2 def gibbs_sampling(mu, sigma, sample_size): samples = [] start = [0, 0] samples.append(start) search_dim = 0 for i in range(sample_size): search_dim = 0 if search_dim == n_dim-1 else search_dim + 1 prev_sample = samples[-1][:] s11 = sigma[search_dim][search_dim] s12 = sigma[search_dim][search_dim - 1] s21 = sigma[search_dim -1 ][search_dim] s22 = sigma[search_dim - 1][search_dim - 1] mu_x = mu[search_dim] mu_y = mu[search_dim-1] _y = prev_sample[search_dim - 1] new_mean = mu_x + s12/float(s22) * (_y - mu_y) new_sigma = s11 - s12/float(s22) * s21 sample_x = np.random.normal(loc=new_mean, scale=np.power(new_sigma, .5), size=1) prev_sample[search_dim] = sample_x[0] samples.append(prev_sample) return np.array(samples) nu = np.ones(2) covariance = np.array([[0.5, 0.5], [0.5, 3]]) samples = gibbs_sampling(nu, covariance, 1000) fig, ax1 = plt.subplots(figsize=(6, 6)) ax1.scatter(sample[:, 0], sample[:, 1], marker=\"o\", facecolor=\"none\", alpha=1., s=30., edgecolor=\"C0\", label=\"Samples\" ) 実行結果

[mathjax] 最尤推定とMAP推定とベイズの定理は繋がっていたので、 記憶が定かなうちに思いの丈を書き出してみるテスト。俯瞰してみると面白い。 あるデータ達(x)が観測されていて、それらは未知のパラメータを持つ確率分布から発生している。 観測されたデータ達(x)を使って、それらのデータを発生させたモデルのパラメータを推定したい。 確率密度関数の中に2つの変数があって。 片方を定数、片方を確率変数として扱うことで2通りの見方ができる。 例えば(n)回のコイントスで(k)回表が出る確率が(theta)だとしてベルヌイ分布の確率密度関数は (k)と(theta)のどちらが確率変数だとしても意味がある。 begin{eqnarray} f(k;theta)=theta^k (1-theta)^{n-k} end{eqnarray} 表が出る確率(theta)が定数だと思って、確率変数(x)の確率密度関数と思う。 単なる確率変数(x)の確率密度関数の中に(theta)という定数がある。尤度。 尤度は確率変数(x)の確率密度関数!。 begin{eqnarray} p(X=x|theta) end{eqnarray} 尤度(p(x|theta))を最大にする(theta)を推定するのが最尤推定。 begin{eqnarray} newcommand{argmax}{mathop{rm arg~max}limits} hat{theta} = argmax_{theta} p(X=x|theta) end{eqnarray} 事後確率と事前確率には関係があって以下のようになる。ベイズの定理。 begin{eqnarray} p(theta|x)=frac{p(x|theta) p(theta)}{p(x)} end{eqnarray} ちなみに、(p(x))は以下のようにしておくとわかりやすい。 同時確率と周辺確率の関係。表を書いて縦、横がクロスするところが同時確率だとして、 縦、横いずれかの方向に同時確率を足し合わせる操作にあたるらしい。 なにか確率変数が独立でなければならない、というのは気にしない。 begin{eqnarray} p(x) = int p(x,theta) dtheta end{eqnarray} なので以下みたいに書き直せる。最後の比例のところは...。 左辺は事後確率分布、右辺は尤度と事前確率分布の積!!。 begin{eqnarray} p(theta|x) &=& frac{p(x|theta) p(theta)}{p(x)} \\ &=& frac{p(x|theta) p(theta)}{int p(x,theta) dtheta} \\ &propto& p(x|theta) p(theta) end{eqnarray} (p(theta))は確率変数(theta)の確率分布。尤度(p(x|theta))は(theta)に対して定数。(x)に対して変数。 ということで、右辺は確率分布(p(theta))を尤度(p(x|theta))を使って変形した確率分布。 で、左辺の(p(theta|x))は右辺の(p(x|theta) p(theta))を定数倍した確率分布。 データを観測していない状態で立てた(p(theta))があって、 観測したデータを使って求めた尤度(p(x|theta))が得られたことで、 左辺の(p(theta|x))が得られた、という状況。 (p(theta|x))は確率変数(theta)の確率分布なので、 最尤推定とベイズの定理を俯瞰してみると、最尤推定が点推定である一方で、 ベイズの定理では確率分布が得られるという具合で異なる。 (観測値が極端なデータだったとき、最尤推定は極端な推定結果が得られるだけだけれども、 ベイズの定理で得られる事後確率分布は確率分布なので様子がわかる..??) 事後確率分布を最大化する(theta_{MAP})を求めるのがMAP推定。(点推定) begin{eqnarray} hat{theta_{MAP}} = argmax_{theta} p(theta|x) end{eqnarray} 尤度をかけて得られた事後分布と同じ形になる便利な分布があって、 観測データ達の分布と対応して決まっている(共役事前分布)。 ベルヌイ分布の共役事前分布はベータ分布。

以前、\"正規分布に従う確率変数の二乗和はカイ二乗分布に従うことの証明\"という記事を書いた。 記事タイトルの通り、正規分布に従う確率変数の二乗和はカイ二乗分布に従う。 [clink url=\"https://ikuty.com/2018/08/01/chi-square-distribution\"] 実際にデータを生成して確かめてみる。 まずは、scipi.stats.chi2.pdfを使って各自由度と対応する確率密度関数を書いてみる。 \"pdf\" ってPortableDocumentFormatではなく、ProbabilityDensityFunction(確率密度関数)。 import numpy as np import matplotlib.pyplot as plt from scipy import stats # 0から8まで1000個のデータを等間隔で生成する x = np.linspace(0, 8, 1000) fig, ax = plt.subplots(1,1) linestyles = [\':\', \'--\', \'-.\', \'-\'] deg_of_freedom = [1, 2, 3, 4] for k in deg_of_freedom: linestyle = linestyles[k-1] ax.plot(x, stats.chi2.pdf(x, k), linestyle=linestyle, label=r\'$k=%i$\' % k) plt.xlim(0, 8) plt.ylim(0, 1.0) plt.legend() plt.show() 定義に従って各自由度ごとに2乗和を足し合わせてヒストグラムを作ってみる。 (ループのリストが自由度) cum = 0 for i in [1,2,3,4]: # 標準正規分布に従う乱数を1000個生成 x = np.random.normal(0, 1, 1000) # 2乗 x2 = x**2 cum += x2 plt.figure(figsize=(7,5)) plt.title(\"chi2 distribution.[k=1]\") plt.hist(cum, 80) 全部k=1ってなってしまった...。左上、右上、左下、右下の順に1,2,3,4。 頻度がこうなので、各階級の頻度の相対値を考えると上記の確率密度関数の形になりそう。 k=1,2が大きくことなるのがわかるし、3,4と進むにつれて変化が少なくなる。

LaravelをAPIサーバにして同期的にsklearnのPCAを実行するアプリを作ってみました。 jQyery/bootstrap/chart.jsがフロント、APIサーバはLaravel+MySQL。 Laravel製APIがGET/POSTに対してPythonコードを実行します(Shellで...)。 exec()でPythonを起動するため無茶苦茶重いし、 ろくにエラーハンドリングできません。 結論から書けば同期的なアプリをこの構造で作るのは無理があります。 バックエンドが無茶苦茶重くてどうせバッチ実行になるのであれば、 上記の問題は結構問題なくなって、これでも良いかなと思い始めます。 MS系のInteroperabilityで、多言語が動的に結合するやつがありますが、 あんな感じでLL言語をglueできれば楽なのになと思います。 PSRの多言語拡張みたいなやつで、PHPからPythonのクラスを使うとか...

[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次元データじゃつまらない...。