Unsupervised Learning: PCA （主成分分析）の流れ
PCA (Principal Component Analysis) とは
データ全体から互いに独立な成分（主成分）を抜き出して、主成分によりもとのデータを表現すること。
分析の中では、データの次元を削減（圧縮）する手段としてもちいられる。また、説明変数の多重共線性の回避の手段として重要。（でも解釈は難しくなるよ）
まずは対象となる特徴量行列として表された各アイテム間の分散共分散行列をとり、その固有値・固有ベクトルを求めることから始まる。
こう書くと大学１年の冬を思い出す。（以下読み飛ばし推奨）
大学入学当初は授業（証明しか解説しない。しかも所々省略する。）と試験（証明には全く触れず、計算問題のみ）のギャップにみごとにはまって落第寸前の日々を送っていた。長時間かけて一つ一つ証明を写経して理解し、自分で解説できるようになってもテストではボロボロ。。。友人に勉強を教えていたのに成績は断然自分のほうが低い。頭の悪さを呪いながら生きていた。
線形代数はそんな時代（１年の冬学期）に苦学して習得し、そしてすぐに証明も含めて忘れ去った分野である。
３年になって、遅ればせながら友人の一人から「授業なんてでないで、教務課でもらった過去問だけ解いてれば優は楽勝だ」と教えられて「え？過去問なんてあるんだ。しかも教務課でもらえるの？というか毎回似たような内容の問題をだすコマなんてあるの？」と目から鱗が落ちる思いだった。そういえば大学受験でも過去問を解いていなかったな。
友人たちは誰からそんなことを教えてもらっていたのか。小ヒルベリーエレジー現象である。
ダブルスクールにこもっていてサークルをおろそかにしていた自分にも原因はあるのだが。。
singapp.hatenablog.com
ちなみに、写経して証明を理解していくという手法は、社会人になってから他の数学・統計の分野を自学するにあたって大変役に立ってます。
GPAの低下という大いなる犠牲を払って学んだことは無駄ではなかった(´;ω;｀)
並行して、過去問を大事にするという学びもCPAとCFAを取得したときに助けてくれた。
余談の余談ですが、当時周りに数多くいた本当に頭の良い人たちは証明も実践も軽くこなしていました。（やっぱり自分の頭が良くないだけやん）
（回想終わり）
さて、本題にもどる
求めた固有ベクトルを使って、特徴量行列を分散共分散行列の固有値の数まで減少させることが主成分分析である。
分散共分散行列の各固有ベクトルの大きさがすなわち固有ベクトルで特徴量を表現する際の分散の大きさとなるので、大きいほうから順番にとっていけばできあがり。
Structure
sklearn.decomposition.PCA(n_components=None, *, copy=True, whiten=False, svd_solver='auto', tol=0.0, iterated_power='auto', random_state=None
主なParameters
| Parameters | Description |
| n_components | 主成分の数 |
| iterated_power | solverの反復回数 |
主なAttributes
| Attributes | Description |
| components_ | 主成分を返す |
| explained_variance_ | 各主成分の大きさ、つまり元の特徴量空間の分散 |
| explained_variance_ratio_ | 上記の分散のパーセンテージ |
| n_components | 主成分の数 |
PCA を １からやってみる
流れ
- 準備： Sample Dataの作成
- Sample Dataを可視化して傾向をみる
- PCAの実行
- 結果の確認
Sample Dataの作成
今日もUCI（いつもお世話になっております）からアワビを取り寄せ。
import requests import io import pandas as pd url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data' data = requests.get(url ,verify = False) df = pd.read_csv(io.StringIO(data.content.decode('utf-8')), header=None) #こちらからNameデータを #http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.names column_names =['Sex', 'Length', 'Diameter', 'Height', 'Whole weight', 'Shucked weight', 'Viscera weight','Shell weight', 'Rings'] df.columns=column_names #Sexをget_dummiesでBooleanに変換 df=pd.get_dummies(df['Sex']).join(df).drop('Sex',axis=1)
Sample Dataの可視化
import seaborn as sns import matplotlib.pyplot as plt g= sns.pairplot(df) g.fig.set_figheight(10) g.fig.set_figwidth(10) plt.show()
説明変数間で相関の高そうなものが散見される
PCAの実行
from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA X = df.drop('Rings', axis=1) y = df['Rings'] #説明変数の標準化 sc= StandardScaler() X_std = sc.fit_transform(X) pca = PCA() pca.fit(X_std) #主成分 X_pca = pca.transform(X_std) pd.DataFrame(X_pca, columns=["PC{}".format(x + 1) for x in range(len(X_pca[1]))]).head()
結果の確認
Explained Variance Ratio (寄与率)
まずは分散の寄与率を確認
pc_df=pd.DataFrame(pca.explained_variance_ratio_, index=["PC{}".format(x + 1) for x in range(len(X_pca[1]))]) print(pc_df) plt.bar([n for n in range(1, len(pca.explained_variance_ratio_)+1)], pca.explained_variance_ratio_) plt.show()
Cumulative Explained Variance (累積寄与率）
次に累積寄与率を図示してみる
plt.plot(np.cumsum(pca.explained_variance_ratio_),marker=".") plt.xlabel('Number of Components') plt.ylabel('Cumulative Explained Variance') plt.grid(True) plt.show()
第２主成分までで85％に届いているので、ここまでの成分を使えば十分であることがわかる。実務上は累積寄与率80％程度で切り捨てられる。
主成分の数を調整してPCAの再実行
結果の確認を受けて、第２主成分までを使う主成分分析をやり直す
pca = PCA(n_components=2) pca.fit(X_std) X_pca = pca.transform(X_std) print(X_pca.shape)
プロット
採用した第１主成分と第２主成分における特徴量の寄与度をプロットすることで、各主成分の解釈を試みる
plt.figure(figsize=(7,7)) for x, y, name in zip(pca.components_[0], pca.components_[1], X.columns[1:]): plt.text(x, y, name) plt.scatter(pca.components_[0], pca.components_[1], alpha=0.8) plt.grid() plt.xlabel("PC1") plt.ylabel("PC2") plt.show()
え、これ元のデータのうちほとんどの特徴が意味をなしてないじゃん。。。（ランクが低すぎる）
第１主成分はM、I（とLength）、F（とその他）という性別を評価した成分。
第２成分はLength、M（とその他）、Iを評価した成分。
どうやらアワビの特徴は性別でほとんど分類されてしまうらしい。あえて性別を除いて主成分分析をしてみたほうがよかったかも。
アルゴリズム思考術 問題解決の最強ツール
- 作者:ブライアン クリスチャン,トム グリフィス
- 発売日: 2017/10/31
- メディア: Kindle版
読むべき時に読めた本。ちょうど最適化アルゴリズムを仕事で使う必要があったので、本書からのインサイトはすぐに役に立った。（自分で組んだわけではなく、納品されたものを検査する側）
こういった本を読んでいて関心するのは、公私含めて普段なにげなく行っている行動の多く（例えば、本棚やメールの整理から配船計画まで）について大学院等の研究結果が公開されていることである。世間では当たり前のことなのだが、自分はほとんど意識していなかった。
これに気が付いたとき、会社で働いているときのモヤモヤ感が一つ晴れた気がした。もっと研究結果を仕事に応用しよう。
一個人に担える発明・発見の量は微々たるものか、もしくはすでに他のだれかが発見して使っているものでありうる。
日本では事務系・技術系と区分された業務と思想が一般的である。一方、発見・発明の必要性は技術系に限らず事務系にも当てはまるのである。
事務系ももっと外と交流して知見を広げることで会社の生産性をぐっと上げることができるのではないだろうか。なにも産学連携とは決して研究開発や技術系分野に限ったもんものではなく、もっと広義に捉えて活用すべきである。また、こういった機能をコンサルばかりに頼っていては会社にknowledgeを蓄積するノウハウが内なわれてしまう。
今更こんなことを考えるとは、潤沢な投資をされた大学で教育を受けた身ながら大変恐縮である。
以下、本の内容のうち覚えておきたいものを抜粋
本書では、代表的な事例に対するアルゴリズムが学術的にはどのように考えられているか、簡単なロジックとともに紹介されている。今後必要になったときに適宜文献を参照したい。大事なのは絶対的な解を見つける方法より少ない計算量による合理的な解見つける方法。
- 最適停止問題
- 最良の伴侶の見つけ方として巷では紹介されている。大事なのは何回「検討」する機会があるか。質が未知もしくは分布が均等ならば37％点以降に今までで最適なものを選ぶ
- ソート
- バブルソートは問題外。O(n)、O(n^2)という試行回数の基準をもちいて比較する。有効なのは比較数え上げ（分割して、分割された部分集合のなかでソートする）
- ファイリング（キャッシュ）
- 最長未使用時間アルゴリズムがよい（最近使ったものを近くへ）
- 汎用スケジューリング
- 各仕事の重要度を処理時間でわり、単位時間重要度で降順ソートする
- アポを取るときはこちらから候補を絞ったほうが相手のメモリ消費を低減できる。
Exploratory Data Analysis - 探索的データ解析
Exploratory Data Analysisとは
データの理解をすること。ざっくりと、
データの型、分布の形状や性質、欠損値や異常値の確認
前処理を含む
- inforによるデータ型、欠損値の確認
- 上記で発見した不具合、不都合の処理
- boxplotを使った異常値の確認
- describeによる基本統計量の確認
- pairplotで散布図に
- 相関係数をheatmapで視覚化
- 必要であればOne-Hot処理
実践編
データ準備
UCIからアワビのとりよせ
Ring（年輪）を予想するモデルを作る
import requests import io import pandas as pd url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.data' data = requests.get(url ,verify = False) df = pd.read_csv(io.StringIO(data.content.decode('utf-8')), header=None) #こちらからNameデータを #http://archive.ics.uci.edu/ml/machine-learning-databases/abalone/abalone.names column_names =['Sex', 'Length', 'Diameter', 'Height', 'Whole weight', 'Shucked weight', 'Viscera weight','Shell weight', 'Rings'] df.columns=column_names
データの観察
眺める
df.head() df.tail() df.info()
列はきれいに集められている模様。アワビのオス・メスってどうやってみわけるんだろう。
欠損値なし。データも数値として認識されている。
データの前処理
データの成型
- 型を変更したり
- to_numeric()
- to_string(),
_ to_datetime()
- データsplitしたり
- テーブルを結合したり
- pd.merge()
- .join()
いろいろやってみる
欠損値処理
- dropna()
- fillna()
Box plotで分布の広がりと異常値を観察
import matplotlib as mpl import matplotlib.pyplot as plt #Ringsは目的変数なのでdrop df.drop('Rings',axis=1).boxplot(figsize=(20,6)) plt.show()
Whole weightの広がりが大きい。異常値であると判断すれば処理。
実力不足により自分では処理できないので、Zepprixさまの記事を参考にさせていただきました。自分でも記事書きたいです。
qiita.com
基本統計量の確認
df.describe()
とてもお手軽。HeightとShuckedのminが気になるところ。
組み合わせ散布図 Pair Plot
import seaborn as sns sns.pairplot(df)
これまたとてもお手軽。Ringsにたいして関係が深そうなものを更なる探索対象にする。
全体的にまとまりにかける。Heightは除いてもよさそう。
相関行列のHeat Map
cor = df.corr() plt.figure(figsize=(10,6)) sns.heatmap(cor, square=True, cmap='plasma', vmax=1,vmin=-1, annot=True) plt.show()
やはり全体的にRingに対する相関が低め。
その後
所与のデータでうまくいきそうなら、次にはモデルの選択とチューニングへ
ダメそうなら特徴量エンジニアリングで目的変数と説明変数の関係性が何か見いだせないか分析
特徴量エンジニアリング
Unsupervised-Larningによって特徴量を操作
- 次元の削減
- 次元が大きすぎると、Machine-Learningの必要時間が増える
- 多重共線性の回避
- 統計的なあれ
k-menas
第三回：Ensemble Model - Bagging (Bootstrap Aggregating)
Baggingとは
Bootstrapping（サンプル集合から復元抽出によりsub-sample setを作成）により作成したsub-sample setそれぞれに対してモデルを作成し、モデルの結果を集約して予測すること。
結果の集約はClassifierの場合は多数決、Regressionの場合は平均値をとる。
Structure
sklearn.ensemble.BaggingClassifier(base_estimator=None, n_estimators=10, *, max_samples=1.0, max_features=1.0, bootstrap=True, bootstrap_features=False, oob_score=False, warm_start=False, n_jobs=None, random_state=None, verbose=0
主なParameters
| Parameters | Description |
| base_estimator | 一つ一つのsub-sample setに用いるモデルを選択 好みがなければDecision Treeがおすすめ（それはつまりRandom Forest..） |
| n_estimators | 何個のsub-setに対してモデルを作成するか |
| max_samples | 元のsampleの何割を一つのsub-sampleとして抽出するか 割合が大きすぎると過学習しがち |
| n_features | sub-sampleで説明変数を何種類使うか 過学習しているときには種類を減らす |
Sample
estimatorとしてK-Neighborを使う
import pandas as pd from sklearn.ensemble import BaggingClassifier from sklearn.neighbors import KNeighborsClassifier from sklearn.model_selection import train_test_split from sklearn.datasets import load_breast_cancer cancer = load_breast_cancer() X= pd.DataFrame(cancer.data, columns=cancer.feature_names) y = pd.Series(cancer.target) X_train, X_test, y_train, y_test = train_test_split(X, y,test_size=0.5 ,random_state = 0) models = {'kNN': KNeighborsClassifier(), 'bagging': BaggingClassifier(KNeighborsClassifier(), n_estimators=50, random_state=0)} scores = {} for model_name, model in models.items(): model.fit(X_train,y_train) scores[(model_name, 'train_score')] = model.score(X_train,y_train) scores[(model_name, 'test_score')] = model.score(X_test,y_test) pd.Series(scores).unstack()
K-Neighbor単体の場合とくらべて多少test_scoreが改善された
回帰モデルの評価
回帰モデルの評価
分類モデルでは分類結果に対しての評価に加えて、計算された確率に対しての評価を行う必要があった。
しかし、回帰モデルでは目的変数が数値であるために計算結果に対して直感的な指標でモデルを評価することができる
主なモデルは
- 平均二乗誤差（MSE: Mean Squared Error ）
- 平均絶対誤差（MAE: Mean Absolute Error）
- 中央値絶対誤差（MedAE: Median Absolute Error）
- 決定係数（R2）
分類モデルの評価はこちら
singapp.hatenablog.com
singapp.hatenablog.com
singapp.hatenablog.com
平均二乗誤差（MSE: Mean Squared Error ）
予測値と目的変数の差（残差）の２乗を足しあげたもの（SSE: Sum of squared errors）をサンプル数でわったもの
線形回帰分析でおなじみ
sklearn.metrics.mean_squared_error(y_true,
y_pred,
*,
sample_weight=None,
multioutput='uniform_average',
squared=True)
crross_validateの引数
- neg_mean_squared_error
平均絶対誤差（MAE: Mean Absolute Error）
残差の絶対値を足しあげて、サンプル数で割ったもの。２乗されていないのでMSEと比べて外れ値の影響を受けにくい
sklearn.metrics.mean_absolute_error(y_true,
y_pred,
*,
sample_weight=None,
multioutput='uniform_average')
crross_validateの引数
- neg_mean_absolute_error
中央値絶対誤差（MedAE: Median Absolute Error）
残差の絶対値の中央値。MAEよりのさらに外れ値に対して強い
sklearn.metrics.median_absolute_error(y_true,
y_pred,
*,
multioutput='uniform_average')
crross_validateの引数
- neg_median_absolute_error
決定係数（R2）
検証データの平均値で予測をした場合の残差平方和（SST: Sum of Squared Total）と残差平方和SSEの比率
平均値での予測をベースにモデルでどれだけ２乗誤差を削れたかを示す指標
sklearn.metrics.r2_score(y_true,
y_pred,
*,
sample_weight=None,
multioutput='uniform_average')
crross_validateの引数
- r2
実践
データの準備
from sklearn.datasets import load_boston boston = load_boston() X = pd.DataFrame(boston.data, columns=boston.feature_names) y = pd.Series(boston.target, name='MEDV')
Sample
特徴として
- Cross ValidationとStandard Scalingを同時に行うためにpiplelineを使う
- Cross Validationではk-foldに代わってRepeated random subsamplingを行う
- データの並びに偏りがあるため、単純なk-foldでは一部のスコアが極端に悪い
- 複数の評価値を同時に計算するためにcross_val_scoreに代わってcross_validateを使う
- 最終値は各検証結果の平均値にする
from sklearn.preprocessing import StandardScaler from sklearn.pipeline import Pipeline from sklearn.model_selection import cross_validate from sklearn.model_selection import RepeatedKFold from sklearn.linear_model import LinearRegression, Ridge from sklearn.tree import DecisionTreeRegressor #検証したいモデルをdict形式に格納 models = { 'LinearRgression': LinearRegression(), 'Ridge': Ridge(random_state=0,max_iter=1200000), 'DecisionTreeRegressor': DecisionTreeRegressor(random_state=0), 'LinearSVR': LinearSVR(random_state=0,max_iter=2400000) } #評価指標の引数をリスト化 scoring = ('neg_mean_squared_error', 'neg_mean_absolute_error', 'neg_median_absolute_error','r2') #PipeLineに入れるScalerオブジェクトを作成 scaler= StandardScaler() #cross_validateに入れるRepeated Random subsumplingのオブジェクトを作成 cv=RepeatedKFold(n_splits=5, n_repeats=5) scores = {} for model_name, model in models.items(): pipeline = Pipeline([('transformer', scaler), ('estimator', model)]) scores[model_name] = pd.DataFrame(cross_validate(pipeline, X, y, cv = cv,scoring=scoring)).mean() scores = pd.DataFrame(scores).unstack().unstack()
Hot Encoder
Hot Encoderとは
特徴量（行列X、多分類のときにはyも)をOne Hot Vectorつまり1つの成分が1で残りの成分が0となるような列ベクトルの集合に変換する機能（classだったり、methodだったり）
pythonとそのLibraryではいろいろ用意されている。
- OneHotEncoder
- get_dummies
- label_binirize
列ベクトルにて表現されている特徴量がカテゴリーの時（数値でない場合）に、特徴を0 or 1で列ごとに１つずつの特徴を表すように変換。そうしないと意味のない数値差（例えばasciiコードの数値差等）が学習データとして採用されてしまう。
カテゴリーの数だけ列を増やすので計算負荷がかかる。
Dataの準備
配管工のおじさんが大好きなキノコのデータを例として準備する
（画像を貼りたいが危険なのでガマン）
Column Nameの格納の仕方が汚かったために思わぬところで前処理の実践をすることに。。。
#UCIのデータベースからキノコデータをダウンロード
url='https://archive.ics.uci.edu/ml/machine-learning-databases/mushroom/agaricus-lepiota.data'
data = requests.get(url, verify=False)
df = pd.read_csv(io.StringIO(data.content.decode('utf-8')),header=None)
#Column Nameを取得
url2='https://archive.ics.uci.edu/ml/machine-learning-databases/mushroom/agaricus-lepiota.names'
name_data = requests.get(url2, verify=False)
name_data= pd.read_csv(io.StringIO(name_data.content.decode('utf-8')),sep='\n',engine='python',error_bad_lines=False)
name_data = name_data.iloc[79:112,0].map(lambda x: str(x).strip())
name_data = name_data[name_data.map(lambda x: x[0].isnumeric())==True]
name_data = name_data.map(lambda x: x.split(':')[0].split('.')[1].strip())
name_data[0] = 'classes'
name_data=name_data.sort_index().reset_index(drop=True)
#データに付加
df.columns=name_data
OneHotEncoder
Class
注意点として、nanをうまく扱えないので事前に処理が必要
変数をいい具合に全体から判断して勝手に詰めてくれる
Structure
sklearn.preprocessing.OneHotEncoder(*,
categories='auto',
drop=None,
sparse=True,
dtype=<class 'numpy.float64'>,
handle_unknown='error')
主なParameters
| Parameters | Description |
| categories | defaultはautoで列の中からカテゴリを自動で作ってくれる。自分で指定もできる。 |
| sparse | True: array, メモリを食う False: Sparse matrix（疎行列） |
Sample
from sklearn.preprocessing import OneHotEncoder henc= OneHotEncoder(sparse=False) hotMushroom=henc.fit_transform(df[['gill-color', 'gill-attachment', 'odor', 'cap-color']]) print(hotMushroom)
get_dummies
nanを勝手に処理してくれる
Structure
pandas.get_dummies(data,
prefix=None,
prefix_sep='_',
dummy_na=False,
columns=None,
sparse=False,
drop_first=False,
dtype=None)
主なParameters
| Parameters | Description |
| data | One-Hotにしたいデータを入れる 統計ではdummy |
| dummy_na | nanをdummyデータとして処理するか |
| sparse | True: array, メモリを食う False: Sparse matrix（疎行列） |
| drop_first | カテゴリー数がk個あった時、k-1個のダミー変数を作成する 統計の自由度の観点から |
Sample
import pandas as pd hotMushroom=pd.get_dummies(df[['gill-color', 'gill-attachment', 'odor', 'cap-color']],sparse=False) hotMushroom.head()
label_binarize
Label(y)用
１列に含まれる多種のLabelを複数列に置き換える。
Structure
sklearn.preprocessing.label_binarize(y,
*,
classes,
neg_label=0,
pos_label=1,
sparse_output=False)
主なParameters
| Parameters | Description |
| y | array 次のclass指定を考慮すると現実的に１列なので目的変数が対象かな |
| classes | カテゴリーを指定 対象列に対してuniqueの結果を流用しできる |
| sparse_output | True: array, メモリを食う False: Sparse matrix（疎行列） |
Sample
from sklearn.preprocessing import label_binarize hotMushroom = label_binarize(df['classes'],classes=df['classes'].unique(),sparse_output=False) hotMushroom
