「日々の学び」一覧

[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くらいが入っていて、対角とあってそう。

[mathjax] 階層型クラスタリングについてわかりやすい解説を聞いたので頭の中にあるものを書き出してみます。 (せっかく聞いたシャープな解説が台無しになってしまっているかと思いますが...) 本当はもっとアクティブなアウトプット方法があれば良いのですが、今のところはブログしかないのと、 WordPressのエントリという単位が、アウトプットの量として個人的に丁度良いなと思いますので。 初学なので解釈の誤りとか検討不足があるかと思いますが、本日も書いていきます！ 試しに実行してみるためのデータを作る面倒さの洗礼を受けまして、 実際に試すのは別記事になります。 [clink url=\"https://ikuty.com/2019/07/16/hierarchical_clustering_by_method/\"] 階層型クラスタリング 小さいクラスタを大きいクラスタで内包する様が木構造(階層構造)になることから、 階層型クラスタリング、という名前が付いているんでしょうかね。 最初にk個のクラスタに分割する！ みたいな問題であればk-means法、 何個のクラスタに分割すべきかはやってみて考えるのであれば、階層型クラスタリング。 木の深さとクラスタ数を対応させるところに見事な美しさを感じます。 クラスタ間距離として以下のWard距離を使うことで空間拡張性と非単調性に優れた結果を得られやすく、 よく傾向が分かっていないデータに対して一発目に適用すべき方法の様子。 begin{eqnarray} d_{ward} (C,C\') = sum_{iin Ccup C\'} || x_i - bar{x}_{Ccup C\'}||^2 - left( sum_{i in C} ||x_i-bar{x}_C||^2 + sum_{iin C\'}||x_i - bar{x}_{C\'}||^2 right) end{eqnarray} 階層型クラスタリングのアルゴリズム 以下の通り。最初、全てのデータポイントをクラスタとみなすところからスタートする。 クラスタとクラスタを併合してより大きいクラスタを作る処理であって、 繰り返しを止めなければ最終的に一つのクラスタまで成長する。 つまり、どこかで止めることで、その時点のクラスタ分割を得ることができる。 全データポイントそれぞれをクラスタとする 全てのデータポイントが初期クラスタ以外に含まれるまで以下の計算を繰り返す 各クラスタ間の距離を計算する。 最もクラスタ間距離が小さいクラスタを併合する。 クラスタを成長させる際に、うまい具合に併合相手を選ばないとよくない。 全体としてクラスタがバラけて欲しいのだけれども、クラスタが隣のクラスタを併合した結果、 次の併合先としてその隣のクラスタが対象となってしまうかもしれない。 そうすると、単調に端っこから順番に成長することになり、 「大きなクラスタ」と「その他の小さなクラスタ」みたいになる。 何らかの距離に基づいてクラスタとクラスタを併合するのだけれども、 クラスタを併合した結果、併合後のクラスタと別のクラスタの距離が、 より短くなってしまう、というケースも起こってしまう。 クラスタ内のデータポイント間の距離、とか、クラスタの重心間の距離とかでは 上記の問題が起こりやすい。 先に距離があって次に併合を考えるのではなく、併合前後の条件を使って距離を作ることで、 上記の問題が起こりにくくなる。その一つがWard距離。 Ward距離 データポイント間の距離はユークリッド距離を使うとして、クラスタ間の距離をどう決めるのか？ データポイントをクラスタ(C)とクラスタ(C\')に分けるとする。分ける前は(Ccup C\')。 (逆? クラスタ(C)とクラスタ(C\')を併合してクラスタ(Ccup C\')を作る)。 クラスタ(C)とクラスタ(C\')の距離(d_{ward}(C,C\'))を以下のように定義する。 begin{eqnarray} d_{ward} (C,C\') = sum_{iin Ccup C\'} || x_i - bar{x}_{Ccup C\'}||^2 - left( sum_{i in C} ||x_i-bar{x}_C||^2 + sum_{iin C\'}||x_i - bar{x}_{C\'}||^2 right) end{eqnarray} 併合後のクラスタ(Ccup C\')における重心との距離の2乗和( sum_{iin Ccup C\'} || x_i - bar{x}_{Ccup C\'}||^2)から、 併合前のクラスタ(C)における重心との距離の2乗和(sum_{i in C} ||x_i-bar{x}_C||^2)、 (C\')における重心との距離の2乗和(sum_{iin C\'}||x_i - bar{x}_{C\'}||^2)を引いたもの。 当たり前だけれども、大きいクラスタよりも小さいクラスタの方が凝集度が小さい。 つまり、併合前のクラスタ(C)、(C\')を併合してクラスタ(Ccup C\')を作るとなると、 (Ccup C\')の凝集度よりも、(C)、(C\')の凝集度の方が必ず小さい。 (C)の凝集度と(C\')の凝集度の和よりも、(Ccup C\')の凝集度の方が大きくなることがポイント。 要するに(Ccup C\')を(C)、(C\')に分割することは、全体の凝集度が小さくなるような操作を行うということ。 無数にある分割の仕方の中から、凝集度の総和が小さくなるような「マシな分割」を決めたい！ (sum_{i in C} ||x_i-bar{x}_C||^2+sum_{iin C\'}||x_i - bar{x}_{C\'}||^2)よりも( sum_{iin Ccup C\'} || x_i - bar{x}_{Ccup C\'}||^2)の方が大きいから、 (d_{ward})は正の値となる。(d_{ward})が大きければそれはマシな分割。 そもそも、データ達が何をもってクラスタを形成しているのかって考えるとかなり微妙で、 人間が解釈するための合理的な根拠を作り出す手段でしかないことがわかった。 データポイント間の\"距離\"を何らかの方法で定義し、\"距離が近しい塊をクラスタと言う\" という話であるとすれば、確かに合理的なのだろうけども。 モヤモヤしているのは、k個のクラスタ集合とk+1個のクラスタ集合を比較したときに、 前者と後者で異なるクラスタに属するデータポイントが発生するところ。 クラスタリングとは単に境界を引く処理なのであって、 それにどういう意味を付けるかは目的次第なのだろうと改めて認識した次第。 そもそも、これらは単なる道具でしかないので、他の手段も同じく 出力結果にどう意味付けするのかは目的次第なんだろうなと。 最初から何個のクラスタに分割すべきか分かっているときはk-means法を使えば良さそう。 今回は、そもそも何個に分割すべきか分かっていないときに使う手法。 [clink url=\"https://ikuty.com/2019/06/29/k-means_idea/\"]

[mathjax] 教師なし学習の問題。クラスタリング。 クラスタの個数を事前に指定するタイプと、自分でクラスタ数を設定できるタイプがあります。 今回、前者のk-means法をアイデアを聞いたので、まとめなおしてみようと思います。 k-means法 (C_1,cdots,C_K)がクラスタ。各々のクラスタにおける中心は(bar{x_k})(centroid)。 各々のクラスタにおけるcentroidとの距離の2乗を全部足して、 それが最小になるようなクラスタ(C_1,cdots,C_K)を求めるという問題。 begin{eqnarray} newcommand{argmin}[1]{underset{#1}{operatorname{arg},operatorname{min}};} hat{c_1},cdots,hat{c_K} = argmin {c_1,cdots,c_K} sum_{k=1}^K sum_{i in C_k} || x_i - bar{x_k}||^2 end{eqnarray} 距離の2乗なので異常値があると異常値に引っ張られて、クラスタが大きく変わってしまう。 centroid(標本平均)も異常値に引っ張られる。外れ値を除去しないといけない。 距離の2乗を計算するので、ベクトル(v_i)の各次元のスケールが揃っていないと、 計算の過程で単位が大きい要素に引っ張られてしまう。 単位の大きさは重要度とは関係ないのに、単位が大きいことが距離に影響を与えることは 避けないといけない。その対応として単位を揃える（正規化）する。 begin{eqnarray} |x_i-bar{x}|^2 = begin{pmatrix} x_{i1} - bar{x_1} \\ x_{i2} - bar{x_2} \\ vdots \\ x_{in} - bar{x_n} \\ end{pmatrix} = sum_{i=1}^{n} sqrt{x_{i1}-bar{x_1}}^2 end{eqnarray} そのままだとNP困難->近似解法 データポイントの個数、centroidの位置がそれぞれ自由だと、 定義通りに実装しようとしても(mathcal{O}(n^3))となって使えない。 begin{eqnarray} newcommand{argmin}[1]{underset{#1}{operatorname{arg},operatorname{min}};} hat{c_1},cdots,hat{c_K} = argmin {c_1,cdots,c_K} sum_{k=1}^K sum_{i in C_k} || x_i - bar{x_k}||^2 end{eqnarray} 以下のようにランダムなクラスタを初期値として設定し、 より良いクラスタに更新することで計算量を下げる。 データ数(n)、クラスタ数(k)だと(mathcal{O}(nk))。 各データポイントにランダムなクラスタを割り当てる 以下の処理を繰り返す 各データポイントに割り当てられたクラスタのクラスタ中心を計算する 各データポイントについてクラスタ中心との距離を計算する 各データポイントについて距離が最小のクラスタを設定する k-meansで実際にクラスタリングするコードを書いてみようと思いましたが、 そろそろ適当なデータを見つけてくるのではなく、シミュレーション用のデータを作るところも やってみようと思います。次回、pythonでコード書きます。

[mathjax] 特異値分解 (M times N)の行列(A)を以下の3つに分解する。 (M times M)の直行行列(V) (M times N)の対角行列(S) (N times N)の直行行列(U) 以下みたいな感じに分解するのが特異値分解。 begin{eqnarray} A &=& VSU^T \\ &=& begin{pmatrix} v_{11} & cdots & v_{1M} \\ vdots & ddots & vdots \\ v_{M1} & cdots & v_{MM} end{pmatrix} begin{pmatrix} lambda_1 & cdots & 0 \\ vdots & ddots & vdots \\ 0 & cdots & lambda_r \\ & & & 0 & \\ & & & & 0 end{pmatrix} begin{pmatrix} u_{11} & cdots & u_{N1} \\ vdots & ddots & cdots \\ u_{1N} & cdots & u_{NN} end{pmatrix} end{eqnarray} (V)は列ベクトル(v_1),(v_2),(cdots),(v_M)を横に並べたもの。 (U^T)は列ベクトル(u_1),(u_2),(cdots),(u_N)を横に並べたものの転値。 (A)のランクは(r)で、(S)は(r)個の(lambda)からなる対角行列。 (lambda_i)は対角行列なので対角の非ゼロの項だけ残す意味で以下のようにかける。 (lambda_i)自体はスカラだから前に寄せられる。 begin{eqnarray} A = sum_{i=1}^{r} lambda_i v_i u_i^T end{eqnarray} 特異値分解のやり方 やりたいことは、(v_i)、(u_i)、(lambda_i)を求めること。 (AA^T)、(A^TA)の固有値、固有ベクトルを計算することで求まる。 主成分分析の対象とするとき、(A)が縦長の行列であることが多い、 つまり高次元のベクトルであることが多いという事実から、 (AA^T)は巨大な行列となり、直接計算することは困難。 以下みたいに(AA^T)を計算する。 begin{eqnarray} AA^T &=& VSU^T(VSU^T)^T \\ &=& VSU^T US^T V^T \\ &=& VSS^TV^T \\ &=& VS_M^2 V^T end{eqnarray} 右から(V)をかける。 begin{eqnarray} AA^T V &=& VS_M^2 \\ AA^T v_i &=& lambda_i^2 v_i end{eqnarray} これは、(AA^T)を固有値(lambda_i^2)、固有ベクトル(v_i)に分解した形。 (固有値、固有ベクトルは計算できる)。 [clink url=\"https://ikuty.com/2019/06/17/eigenvalue-eigenvector/\"] 同様にして、(A^TA)を計算する。 (A^TA)の固有値、固有ベクトルの計算により(u_i)と(lambda_i)を求められる。 begin{eqnarray} A^TA &=& (VSU^T)^T VSU^T \\ &=& US_N^TS_nU^T \\ &=& US_N^2 U^T \\ (A^TA) U &=& US_N^2 \\ (A^TA) u_i &=& lambda_i^2 u_i end{eqnarray} (AA^T)と異なり、(A^TA)は大きくなく、容易に計算が可能。 (v_i)の求め方 (AA^T v_i = lambda_i^2 v_i)は大きすぎて計算できないが、(A^TA)を使って(AA^T)を計算できる。 以下のように式変形する。(u_k)の正規直行性を使う。 begin{eqnarray} A &=& sum_{i=1}^{r} lambda_i v_i u_i^T \\ Au_k &=& sum_{i=1}^{r} lambda_i v_i u_i^T u_k \\ &=& lambda_k v_k u_k ^T u_k \\ &=& lambda_k v_k \\ v_k &=& frac{1}{lambda_k} Au_k end{eqnarray} (A^TA)の特異値分解は計算可能。以下から(hat{u_i})と(hat{lambda_i})が求まる。 begin{eqnarray} (A^TA) hat{u_i} &=& hat{lambda}_i^2 hat{u_i} end{eqnarray} 得られた(hat{u_i})と(hat{lambda_i})を (v_k = frac{1}{lambda_k} Au_k)に入れることで、(v_k)が求まる。