記事・メモ一覧

[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)が求まる。

[mathjax] 線形代数もやりなおします。流石に時間ないので微妙に分かるところまでですが。 第1弾は固有値、固有ベクトル。 具体例を使って追ってみるやつ 行列(A)とベクトル(vec{v_1})をデタラメに選んで内積を取ってみる。 (A=begin{pmatrix} 2 & 5 \\ 3 & 7 end{pmatrix})。(vec{v_1}=begin{pmatrix} 1 \\ 4end{pmatrix})。 begin{eqnarray} vec{v_2} = Avec{v} = begin{pmatrix} 2 & 5 \\ 3 & 7 end{pmatrix} begin{pmatrix} 1 \\ 4end{pmatrix} = begin{pmatrix} 2 times 1 + 5 times 4 \\ 3 times 1 + 7 times 4 end{pmatrix} = begin{pmatrix} 22 \\ 31 end{pmatrix} end{eqnarray} (vec{v_1})と(vec{v_2})は向きが違うし大きさも違う。 (vec{v_1})に(A)をかけることで、回転して伸ばしてる。 (vec{v_1})と(vec{v_2})の向きが同じになるように(vec{v_1})を選べないか...。 下みたいな組み合わせだと、左と右が同じ向きになる。 begin{eqnarray} begin{pmatrix} 2 & 5 \\ 3 & 7 end{pmatrix} begin{pmatrix}0.575249 \\ 0.817978end{pmatrix} = begin{pmatrix}5.2403 \\ 7.4515end{pmatrix} end{eqnarray} 伸ばす分を外に出すと以下みたいになる。 begin{eqnarray} begin{pmatrix} 2 & 5 \\ 3 & 7 end{pmatrix} begin{pmatrix}0.575249 \\ 0.817978end{pmatrix} = 9.10977 begin{pmatrix}0.575249 \\ 0.817978end{pmatrix} end{eqnarray} 長さを伸ばす分の係数(9.10977)が行列(A)の固有値。 begin{pmatrix}0.575249 \\ 0.817978end{pmatrix}が行列(A)の固有ベクトル。 次数が多いやつ (A)は(n times n)の正方行列。固有ベクトル(vec{x})は(1 times n)。固有値(lambda)はスカラー。 begin{eqnarray} A vec{x} &=& lambda vec{x} end{eqnarray} 以下みたいに変形して(vec{x})と(lambda)を求める。 begin{eqnarray} (A-lambda E) vec{x} &=& vec{0} end{eqnarray} ((A-lambda E))に逆行列があると(Evec{x} = vec{0})とかになってゼロベクトルじゃない(vec{x})を求められないから、 逆行列が存在しない、という式を立てるらしい。。 \"逆行列が存在しない\"だけだと、(vec{x})も(lambda)も複数ありそうだけど絶対に1つしか存在しないらしい。 逆行列が存在しないのは、行列式がゼロ、とするらしい。 begin{eqnarray} det(A-lambda E) = 0 end{eqnarray} 手計算が合わないのでそれ以降は省略... まとめ 行列(A)の固有値、固有ベクトルは、 (A vec{x} = lambda vec{x})となる(vec{x})と(lambda)のセットを求めることでした。 向き、大きさの変換を表す(A)を、向きを表すベクトル(vec{x})と大きさ(lambda)で表しなおす..ような。

[mathjax] タイトルの通り。Losso回帰と違って損失関数を偏微分するだけで出来そうなのでやってみる。 Ridge回帰は線形回帰の1種だけれども、損失関数として最小二乗法をそのまま使わず、 (L_2)ノルムの制約を付けたものを使う((L_2)正則化)。 データとモデル 教師データ(boldsymbol{y})、訓練データ(boldsymbol{x})があるとする。 (または目的変数(boldsymbol{y})、説明変数(boldsymbol{x})があるとする。) 例えば(p)次の属性データが(n)個あり、それらと結果の対応が分かっている状況。 begin{eqnarray} boldsymbol{y} &=& begin{pmatrix} y_1 \\ y_2 \\ vdots \\ y_p end{pmatrix} , boldsymbol{x} &=& begin{pmatrix} x_{11} & x_{21} & cdots & x_{n1} \\ x_{12} & x_{22} & cdots & x_{n2} \\ vdots & vdots & ddots & vdots \\ x_{1p} & x_{2p} & cdots & x_{np} end{pmatrix} end{eqnarray} モデルは以下。特徴ベクトル(boldsymbol{w})は訓練データの重み。 特徴空間において損失を最小化する特徴ベクトルを求める問題。 begin{eqnarray} boldsymbol{y} &=& boldsymbol{w} boldsymbol{x} + k \\ boldsymbol{w} &=& begin{pmatrix} w_1 & w_2& cdots &w_p end{pmatrix} end{eqnarray} 損失関数 普通の2乗損失に正則化項((L_2)ノルムを定数倍した値)を付けたものを損失関数として利用する。 正則化項の係数はハイパーパラメータとして調整する値。逆数なのはsklearnに従う。 begin{eqnarray} L(boldsymbol{w}) = |boldsymbol{y} - boldsymbol{w} boldsymbol{x}|^2 +C |boldsymbol{w}|^2 end{eqnarray} 特徴ベクトルは以下。(mathjaxでargminが出せない...) begin{eqnarray} newcommand{argmin}[1]{underset{#1}{operatorname{arg},operatorname{min}};} boldsymbol{w} = argmin w L(boldsymbol{w}) = argmin w |boldsymbol{y} - boldsymbol{w} boldsymbol{x}|^2 + C |boldsymbol{w}|^2 end{eqnarray} 特徴ベクトルを求める 勾配=0と置けば上の式の解を得られる。 損失関数が微分可能だからできる技。 begin{eqnarray} frac{partial L(boldsymbol{w})}{partial boldsymbol{w}} &=& 2 boldsymbol{w}^T (boldsymbol{y} - boldsymbol{w} boldsymbol{x}) + C boldsymbol{w} \\ &=& 0 end{eqnarray} 変形する。 begin{eqnarray} 2 boldsymbol{x}^T (boldsymbol{x}boldsymbol{w}-boldsymbol{y}) + C boldsymbol{w} &=& 0 \\ boldsymbol{x}^T (boldsymbol{x}boldsymbol{w}-boldsymbol{y}) + C boldsymbol{w} &=& 0 \\ boldsymbol{x}^T boldsymbol{x} boldsymbol{w} -boldsymbol{x}^T boldsymbol{y} + Cboldsymbol{w} &=& 0 \\ (boldsymbol{x}^T boldsymbol{x} +C E) boldsymbol{w} &=& boldsymbol{x}^T boldsymbol{y} \\ boldsymbol{w} &=& (boldsymbol{x}^T boldsymbol{x} + C E)^{-1} boldsymbol{x}^T boldsymbol{y} end{eqnarray} テストデータを作る 練習用にsklearnのbostonデータを使ってみる。 ボストンの住宅価格が目的変数、属性データが説明変数として入ってる。 import pandas as pd import numpy as np from pandas import Series,DataFrame import matplotlib.pyplot as plt from sklearn.datasets import load_boston boston = load_boston() boston_df = DataFrame(boston.data) boston_df.columns = boston.feature_names print(boston_df.head()) boston_df[\"PRICE\"] = DataFrame(boston.target) # CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT PRICE # 0 0.00632 18.0 2.31 0.0 0.538 6.575 65.2 4.0900 1.0 296.0 15.3 396.90 4.98 24.0 # 1 0.02731 0.0 7.07 0.0 0.469 6.421 78.9 4.9671 2.0 242.0 17.8 396.90 9.14 21.6 # 2 0.02729 0.0 7.07 0.0 0.469 7.185 61.1 4.9671 2.0 242.0 17.8 392.83 4.03 34.7 # 3 0.03237 0.0 2.18 0.0 0.458 6.998 45.8 6.0622 3.0 222.0 18.7 394.63 2.94 33.4 # 4 0.06905 0.0 2.18 0.0 0.458 7.147 54.2 6.0622 3.0 222.0 18.7 396.90 5.33 36.2 散布図行列を表示してみる。 PRICEと関係がありそうなZN,RM,AGE,DIS,LSTATの5個を使ってみる。 pg = sns.pairplot(boston_df) plt.show() pg.savefig(\'boston_fig.png\') 特徴ベクトルを自力で計算する これを自力で計算してみる。(C=0.01)、(C=0)、(C=100)としてみた。 begin{eqnarray} boldsymbol{w} &=& (boldsymbol{x}^T boldsymbol{x} + C E)^{-1} boldsymbol{x}^T boldsymbol{y} end{eqnarray} X_df = boston_df.drop(columns=[\'CRIM\',\'INDUS\',\'CHAS\',\'NOX\',\'RAD\',\'TAX\',\'PTRATIO\',\'B\',\'PRICE\']) X = X_df.values y = boston.target.T C1 = 0.01 C2 = 0 C3 = 100 e = np.identity(5) w1 = np.dot( np.linalg.inv(np.dot(X.T , X) + C1 * e), np.dot(X.T,y)) w2 = np.dot( np.linalg.inv(np.dot(X.T , X) + C2 * e), np.dot(X.T,y)) w3 = np.dot( np.linalg.inv(np.dot(X.T , X) + C3 * e), np.dot(X.T,y)) print(w1) # [ 0.05338557 5.40396159 -0.01209002 -0.83723303 -0.63725397] print(w2) # [ 0.05338539 5.40403743 -0.01209427 -0.83728837 -0.63725093] print(w3) # [ 0.05612977 4.76664789 0.02374402 -0.38576708 -0.66137596] (C=0)のとき、つまり最小二乗法のとき。 sklearnを使う sklearnのridge回帰モデルを使うと以下みたいになる。 from sklearn.linear_model import Ridge from sklearn.model_selection import train_test_split Xf_train,Xf_test,yf_train,yf_test = train_test_split(X,y,random_state=0) ridge = Ridge().fit(Xf_train,yf_train) print(f\"accuracy for training data:{ridge.score(Xf_train,yf_train):.2}\") print(f\"accuracy for test data:{ridge.score(Xf_test,yf_test):.2f}\") # accuracy for training data:0.68 # accuracy for test data:0.58 print(ridge.coef_) # [ 0.06350701 4.3073956 -0.02283312 -1.06820241 -0.73188192] 出てきた特徴ベクトルを並べてみる 自力で計算したものとsklearnに計算してもらったものを並べてみる。 似てるのか似ていないのかよくわからない .. けど、RMの寄与度が高いというのは似ている。 # 自力で計算 (C=100) # [ 0.05612977 4.76664789 0.02374402 -0.38576708 -0.66137596] # sklearnで計算 # [ 0.06350701 4.3073956 -0.02283312 -1.06820241 -0.73188192] 自力で計算したモデルの正答率を求めてみないとなんとも... そして、正規化項の係数の大小がどう影響するのか、あまり良くわからなかった..。 (L_2)ノルムの制約を付けると、パラメタの大小が滑らかになると言いたかったのだけども。 あと、訓練データに対して68%、テストデータに対して58%という感じで、 大して成績が良くない...。　

集合関数 集合関数。ndarrayから重複を取り除きsortした結果を返す。 2dであってもその中から要素を抜き出して1dにする。 hoges = np.array([\"hoge\",\"fuga\",\"hoge\",\"fuga\"]) print(np.unique(hoges)) # [\'fuga\' \'hoge\'] fugas = np.array([[\"hoge\",\"fuga\",\"hoge\",\"fuga\"],[\"hoge2\",\"fuga2\",\"hoge2\",\"fuga2\"]]) print(np.unique(fugas)) # [\'fuga\' \'fuga2\' \'hoge\' \'hoge2\'] ファイルI/O pandasを使わずとも、NumPyだけでファイルI/Oができる。 以下でhoges.npyという無圧縮バイナリファイルが作られる。 それを読み込んで出力する。 hoges = np.array([\"hoge\",\"fuga\",\"hoge\",\"fuga\"]) np.save(\'hoges.npy\',hoges) fugas = np.load(\'hoges.npy\') print(fugas) # [\'hoge\' \'fuga\' \'hoge\' \'fuga\'] 複数の配列を同時に書き込むこともできる。 キーワードを指定して書き込む。キーワードを指定して1つずつ読み込む。 読み込む時はキーワードを指定して参照したときに遅延ロードされる。 hoges = np.array([\"hoge1\",\"fuga1\",\"hoge1\",\"fuga1\"]) fugas = np.array([\"hoge2\",\"fuga2\",\"hoge2\",\"fuga2\"]) np.savez(\'hogefuga.npz\', hoges=hoges, fugas=fugas) hogefugas = np.load(\'hogefuga.npz\') hoges_l = hogefugas[\'hoges\'] fugas_l = hogefugas[\'fugas\'] print(hgoes_l) # [\'hoge1\' \'fuga1\' \'hoge1\' \'fuga1\'] print(fugas_l) # [\'hoge2\' \'fuga2\' \'hoge2\' \'fuga2\']

線形サポートベクトル分類器で画像認識する流れを理解したので、 定着させるために記事にしてみます。 当然、モデルの数学的な理解がないとモデルを解釈することは不可能だし、 正しいハイパーパラメータを設定することも不可能なので、数学的な理解は不可欠。 NumPy、pandas、matplotlibに慣れないと、そこまで行くのに時間がかかります。 こちらはPythonプログラミングの領域なので、数こなして慣れる他ないです。 機械学習用のサンプル画像で有名なMNISTを使ってNumPy、pandasの練習。 手書き文字認識用の画像データを読み込んでみる。サイズは28x28。各々1byte。 MNISTの手書き文字認識画像の読み込み まず読み込んでみて、データの形を出力してみる。 X_trainは、要素が3個のTupleが返る。3次。 1番外が60000。28x28の2次のndarrayが60000個入っていると読む。 1枚目の画像データはX_train[0]によりアクセスできる。 import tensorflow as tf minst = tf.keras.datasets.mnist (X_train,y_train),(X_test,y_test) = mnist.load_data() print(X_train.shape) # (60000, 28, 28) y_trainは要素が1個のTupleが返る。1次。 1枚目から60000枚目までの画像が0から9のいずれに分類されたかが入っている。 y_train[0]が4なら、1枚目の画像が4に分類された、という意味。 print(y_train.shape) # (60000,) データセットの選択 X_train,y_train、X_test,y_testから、値が5または8のものだけのViewを取得する。 そのために、まず値が5または8のものだけのインデックスを取得する。 NumPyのwhereはndarrayのうち条件を満たす要素のインデックスを返す。 X_trainに入っている60000件の2d arrayのうち、 値が5または8のインデックス(0-59999)を取得するのは以下。 index_train = np.where((X_train==5)|(X_train==8)) print(index_train) # (array([ 0, 11, 17, ..., 59995, 59997, 59999]),) index_test = np.where((X_test==5)|(X_test==8)) print(index_test) # (array([ 8, 15, 23, ..., 9988, 9991, 9998]),) インデックスを使って絞り込む。 X_train,y_train = X_train[index_train],y_train[index_train] X_test,y_test = X_test[index_test],y_test[index_test] print(X_train.shape) # (11272, 28, 28) print(X_test.shape) # (1866, 28, 28) 前処理 0-255の間の値を0-1の間の値に変換する（正規化）。 28x28の画像(2darray)を1x784(1darray)に整形する（平坦化)。 X_train,X_test = X_train / 255.0, X_test / 255.0 X_train = X_train.reshape(X_train.shape[0], X_train.shape[1] * X_train.shape[2]) X_test = X_test.reshape(X_test.shape[0], X_test.shape[1] * X_test.shape[2]) ベストなハイパーパラメータの選択 線形サポートベクトル分類器を作成する。 from sklearn.svm import LinearSVC linsvc = LinearSVC(loss=\"squared_hinge\",penalty=\"l1\",dual=False) 線形サポートベクトル分類器のハイパーパラメータCの選択 逆正則化パラメータCをGridSearchCVで探す。MBP2013Laterで学習(fit)に5分くらいかかった。 GridSearchCVからはC=0.2がbestと返ってくる。 from sklearn.model_selection import GridSearchCV param_grid = {\"C\":[0.025,0.05,0.1,0.2,0.4]} model = GridSearchCV(estimator=linsvc, param_grid=param_grid,cv=5,scoring=\"accuracy\",return_train_score=True) model.fit(X_train,y_train) print(model.cv_results_[\"mean_train_score\"]) # array([0.96291693, 0.96775192, 0.97059085, 0.97340754, 0.97626859]) print(model.cv.results_[\"mean_test_score\"]) # array([0.95626331, 0.95990064, 0.96158623, 0.9625621 , 0.96105394]) print(model.best_params_) # {\'C\': 0.2} 学習、精度評価 C=0.2を使って新しく学習させる。 linsvc = LinearSVC(loss=\"squared_hinge\",penalty=\"l1\",dual=False,C=0.2) linsvc.fit(X_train,y_train) 訓練データ、テストデータに対して正答率を求める。 訓練データについて97.2%、テストデータについて96.2%。 過学習すると訓練データが高くテストデータが低くなる。 from sklearn.metrics import accuracy_score pred_train = linsvc_best.predict(X_train) acc = accuracy_score(y_true = y_train,y_pred = pred_train) print(acc) # 0.9723207948899929 pred_test = linsvc_best.predict(X_test) acc = accuracy_score(y_true = y_test,y_pred = pred_test) print(acc) # 0.9619506966773848 モデルの解釈可能性 [mathjax] 線形SVMの決定境界(f(x))の係数をヒートマップっぽく表示して、どの係数を重要視しているかを確認する。 基本的に真ん中に画像が集まっているので、28x28の隅は使わないのが正しそう。 正則化パラメータによって係数の大きさを制御しているため、正則化パラメータを変えると係数が変わる。 今回のは(L_1)正則化なので、係数が0のものが増える..らしい（..別途調べる..)。 (f(x) = w_0 + w_1 x_1 + w_2 x_2 + cdots w_{784} x_{784}) import matplotlib.pyplot as plt weights = linsvc_best.coef_ plt.imshow(weights.reshape(28,28)) plt.colorbar() plt.show()

universal functions ndarrayの全ての要素に対して基本的な計算を実行する。 以下オペランドが1つの単項universal functions。 abs,sqrt,square,exp,log,sign,ceil,floor,rint,modf,isnan,sin,cos,arcsin,arccosなどがある。 array = np.arange(10) print(array) # [0 1 2 3 4 5 6 7 8 9] sqrt = np.sqrt(array) print(sqrt) # [0. 1. 1.41421356 1.73205081 2. 2.23606798 # 2.44948974 2.64575131 2.82842712 3. ] exp = np.exp(array) print(exp) # [1.00000000e+00 2.71828183e+00 7.38905610e+00 2.00855369e+01 # 5.45981500e+01 1.48413159e+02 4.03428793e+02 1.09663316e+03 # 2.98095799e+03 8.10308393e+03] 以下、オペランドが2つの2項universal functions。 いずれかのうち最大の値を残すmaximum()。 add,subtract,divide,power,maximum,minimum,copysign,greater,lessなどがある。 x = np.random.randn(10) y = np.random.randn(10) print(x) # [ 1.3213258 0.12423666 -1.45665939 -1.49766467 -0.6129116 2.00056744 # -0.00816571 0.63247747 0.29497652 0.80000291] print(y) # [-0.76739214 0.95151629 0.03208859 0.40641677 0.82635027 1.01773826 # 0.75601178 0.25200147 1.59929321 0.6251983 ] z = np.maximum(x,y) print(z) # [1.3213258 0.95151629 0.03208859 0.40641677 0.82635027 2.00056744 # 0.75601178 0.63247747 1.59929321 0.80000291] [mathjax] matplotlibにndarrayを引数として渡せば簡単にプロットできる。 (z=sqrt{x^2+y^2})をプロットしてみる。 import numpy as np import matplotlib.pyplot as plt points = np.arange(-5,5,0.01) xs,ys = np.meshgrid(points,points) z = np.sqrt(xs**2 +ys**2) plt.imshow(z, cmap=plt.cm.gray) plt.colorbar() plt.title(\"Image plot\") plt.show() 3項演算子 where マスクの論理値に従って2つのndarrayのうちいずれかの値を選択してリストに書く。 3項演算子を使ってPythonのlistに入れる方法は以下。 xa,xbはndarrayだが最終的なr1はPythonオブジェクト。 import numpy as np xa = np.array([1,2,3,4,5]) xb = np.array([6,7,8,9,10]) cnd = np.array([True,True,False,False,False]) r1 = [(x if c else y) for x,y,c in zip(xa,xb,cnd)] print(r1) 対して、ndarrayに対して直に3項演算子を実行するwhereがある。 import numpy as np xa = np.array([1,2,3,4,5]) xb = np.array([6,7,8,9,10]) cnd = np.array([True,True,False,False,False]) r2 = np.where(cnd,xa,xb) print(r2) 数学関数,統計関数,次元削減 (n)次のndarrayをある軸について集計して(n-1)次のndarrayにする。 集計方法としていくつかの数学関数、統計関数が用意されている。 以下5x4(2次)のndarrayについて、それぞれの列について平均を取り4列(1次)のndarrayにしている。 さらに列の平均を取りスカラーにしている。 import numpy as np ary = np.random.randn(5,4) print(ary) # [[-1.84573174 1.84169514 1.43012623 -0.5416877 ] # [-1.03660701 0.63504086 -0.12239017 -0.77822113] # [ 0.1711323 -0.16660851 -0.7928288 1.17582814] # [-0.29302267 -0.23316282 1.70611457 0.53870384] # [-0.46513289 -1.12207588 0.01930695 0.49635739]] print(ary.mean(axis=0)) # [-0.6938724 0.19097776 0.44806576 0.17819611] print(ary.mean(axis=1)) # [ 0.22110048 -0.32554436 0.09688078 0.42965823 -0.26788611] print(ary.mean()) # 0.030841804893752683

大規模高速計算を前提にC言語との接続を前提にしていて、配列処理に寄せることになる。 ndarrayで確保するメモリはPythonとは別(プロセス?)で確保される。 一通り流してみる。 shape()で配列の形を応答する。2行3列。 import numpy as np data = np.random.randn(2,3) shape = data.shape print(shape) print(data) # (2, 3) # [[ 0.79004157 0.45749364 0.90854549] # [-1.91791968 2.80050094 -0.60338724]] ndarrayを作る ndarrayを作る方法は以下。 data1 = [1,2,3,4,5] data2 = [6,7,8,9,10] data = np.array([data1,data2]) print(data) # [[ 1 2 3 4 5] # [ 6 7 8 9 10]] rng = np.arange(5) print(rng) # [0 1 2 3 4] ones = np.ones((5,5)) print(ones) # [[1. 1. 1. 1. 1.] # [1. 1. 1. 1. 1.] # [1. 1. 1. 1. 1.] # [1. 1. 1. 1. 1.] # [1. 1. 1. 1. 1.]] # 零行列 zeros = np.zeros((3,5)) print(zeros) # [[0. 0. 0. 0. 0.] # [0. 0. 0. 0. 0.] # [0. 0. 0. 0. 0.]] # 未初期化の配列確保 empties = np.empty((5,3)) print(empties) # [[-1.72723371e-077 -1.72723371e-077 2.24419447e-314] # [ 2.24421423e-314 2.24421423e-314 2.24563072e-314] # [ 2.24421559e-314 2.24563072e-314 2.24421570e-314] # [ 2.24563072e-314 2.24421558e-314 2.24563072e-314] # [ 2.24421562e-314 2.24563072e-314 2.24421577e-314]] # 指定値で埋める fulls = np.full((2,3),5) print(full) # [[5 5 5] # [5 5 5]] # 単位行列 identities = np.identity(5) print(identities) # [[1. 0. 0. 0. 0.] # [0. 1. 0. 0. 0.] # [0. 0. 1. 0. 0.] # [0. 0. 0. 1. 0.] # [0. 0. 0. 0. 1.]] ndarrayのデータ型 ndarrayで確保されるメモリのデータ型。 実際に型に従ってメモリが確保されているため、簡単にCに渡せる。 ary = np.array((1,2,3),dtype=np.float64) print(ary) # [1. 2. 3.] # float64をint32でキャスト ary_int = ary.astype(np.int32) print(ary_int) # [1 2 3] # キャストできないとコケる ary_str = np.array([\'hoge\',\'fuga\']) ary_str_int = ary_str.astype(np.int32) # ValueError: invalid literal for int() with base 10: \'hoge\' ベクトル演算 配列に寄せる醍醐味。Pythonに数値計算用のオペランドが用意されていることがあって、 割と自然に書ける。 ary = np.array([[1,2,3],[4,5,6]]) print(ary * ary) # [[ 1 4 9] # [16 25 36]] print(ary - ary) # [[0 0 0] # [0 0 0]] print(ary * 2) # [[ 2 4 6] # [ 8 10 12]] print(ary ** 2) # [[ 1 4 9] # [16 25 36]] スライスとView 巨大なメモリへのアクセス高速化のために、np.arrayに対するスライスによるアクセスは、 同じメモリを指すViewを返す。Viewに対する操作は元のメモリを変更する。 Copyする場合は明示的にCopyをする必要がある。 ary = np.arange(10) print(ary) # [0 1 2 3 4 5 6 7 8 9] ary[5] = 500 print(ary) # [ 0 1 2 3 4 500 6 7 8 9] ary[3:5] = 999 print(ary) # [ 0 1 2 900 900 500 6 7 8 9] copied = ary.copy() print(copied) # [ 0 1 2 900 900 500 6 7 8 9] 2次元のnp.array。要素へのアクセスの仕方は2通り。 ary2d = np.array([[1,2,3],[10,20,30]]) print(ary2d) # [[ 1 2 3] # [10 20 30]] print(ary2d[1]) # [10 20 30] print(ary2d[1][0]) # 10 print(ary2d[1,0]) # 10 n次元arrayへスカラーでインデックス参照するとn-1次元が戻る。 スライス参照はn次元が戻る。 ary2d = np.array([[1,2,3],[10,20,30],[100,200,300]]) print(ary2d[1]) # [10 20 30] print(ary2d[:2]) # [[ 1 2 3] # [10 20 30]] print(ary2d[1,:2]) # [10,20] Viewの選択 ndarrayから欲しいViewを選択するために色々と条件をつけられる。 例えば、bool index参照。 data = np.random.randn(7,4) print(data) # [[-0.69179761 -1.30790477 1.7224557 -0.67436315] # [ 0.45457462 0.24713663 -0.84619583 -0.31182853] # [-1.36397651 0.51770088 -1.8459593 -1.75146057] # [ 2.38626251 -0.4747874 -0.49951212 0.61803437] # [ 1.00048197 1.21838773 -0.4828001 0.9952139 ] # [ 0.17838262 1.687342 0.81501139 -1.12800811] # [ 0.65216988 -2.57185067 0.29802975 0.28870091]] recs = np.array([\'apple\',\'orange\',\'banana\',\'mountain\',\'river\',\'moon\',\'snow\']) print(recs==\'mountain\') # [False False False True False False False] print(data[recs==\'mountain\']) # [[ 2.38626251 -0.4747874 -0.49951212 0.61803437]] reshape reshape()を使って行列の形を変える。例えば1x15のndarrayを3x5のndarrayに変換。 もちろんCopyではなくView。これは頻出っぽい。 ちなみに、転値は専用のメソッド(T)が用意されている。 data1 = np.arange(15) print(data1) # [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14] data2 = data1.reshape(3,5) print(data2) # [[ 0 1 2 3 4] # [ 5 6 7 8 9] # [10 11 12 13 14]] data3 = data2.T print(data3) # [[ 0 5 10] # [ 1 6 11] # [ 2 7 12] # [ 3 8 13] # [ 4 9 14]]