PythonやRを用いたアルゴリズム取引・株価分析まとめ

はじめに

まわりでシステムトレードや株価の機械学習による予測などに関心が高まってきたので、私も少し調べてみようと思いPythonやRで行われた分析・実装の事例を集めてみました。自分の資産を突っ込む気にはなれないですが、事例を知っておくだけ知っておきたいですね。

調査法

Google検索において以下のクエリで上位に表示されたサイトを中心にまとめました。
「python 機械学習 株価」
「python 機械学習 為替」
「python アルゴリズム取引」
「python machine learning stock price」
(同様にRも調べました。)

アルゴリズムトレードの理論から学ぶ

システムトレードで億万長者になるぞ! coursera で Computational Investing Part I を受けた
こちらはCourseraで開かれていたシステムトレード向けの講座を受けた方のブログです。Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Riskはもっとも重要なアルゴリズムトレードの本とされているようです。出版年が1999年と結構古いですが。

Rを用いた株価データのモニタリング

Stock Analysis using R
こちらはアルゴリズム取引とかではないですが、Rで株価を簡単に取得できるquantmodライブラリについて紹介されています。

以上のコードを実行しただけで時系列データがボックスプロットで描写されました。非常に便利そうです。複数のグラフを並べてディスプレーで見てみたいですね。図ではテクニカル分析で使われるボリンジャーバンドを一発で出してくれています。傾向を掴むための分析に役立つのではないでしょうか。

MSFT_stockPrice

RSIで株価の連動性を見る

過去のデータからビッグデータ分析で株価を予測する
ここでは、主にpandasを用いて株価の分析を行っています。RSI(Relative Strength Index)という「一定の期間変動幅の中でどれ位株価が上昇しているのか、下落しているのかをはかるもの」を計算して日経平均との相関を見ています。データさえあれば、個別株の市場連動性を見るぶんにはpandasで十分に分析できそうです。

アルゴリズムトレードのシステム構成

自動トレードボット
こちらでは、仮想通貨取引のための自動トレードボットの作成のための手順などが書かれています。
工程としては
・仮想通貨価格データ取得
・バックテストの実施
・明日の価格予想
・学習パラメータの最適化
・結果のメール送信
・APIを使った成行トレード
・ジョブのスケジューリング
などとなっているようです。
完全自動化しようと思うと、作るのは大変そうですね。

WEB屋の自分が機械学習株価予想プログラムを開発した結果
こちらはよりシンプルで
「チャートを見る限り、誰がどうみても今日値上がりする銘柄」を検索して、毎朝Slackに通知してくれるシステム
をPythonで作成されています。最終的に資産を溶かす形になっているようです。

アルゴリズムトレード向けのPythonライブラリ

pythonのアルゴリズムトレードライブラリ
こちらでは
zipline( http://www.zipline.io/
PyAlgoTrade( http://gbeced.github.io/pyalgotrade/ )
pybacktest( https://github.com/ematvey/pybacktest )
backtrader( https://github.com/mementum/backtrader )
というライブラリが紹介されています。

遺伝的アルゴリズムで為替の自動売買

pythonと遺伝的アルゴリズムで作るFX自動売買システム その1
遺伝的アルゴリズムでFX自動売買 その2 進化する売買AIの実装
遺伝的アルゴリズムでFX自動売買 その3 OandaAPIで実際に取引
こちらは、外国為替の取引をPythonで自動化させた試みです。最終的に負けてしまってはいますがシステム構成についても詳しく書かれているので、勉強になりました。データはOANDAというサービスが提供しているAPIを用いて取得し、GA(遺伝的アルゴリズム)を使って為替の売買タイミングを決めているようです。実際の売買にもOandaAPIというものを利用して完全に自動化させています。

決定木による株価リターンの予測

機械学習で未来を予測する – scikit-learn の決定木で未来の株価を予測
こちらはPythonでscikit-learnを用いて、決定木による株価の予測をされています。目的変数としてはリターンインデックスをおいています。前処理においてpandasが使われています。

ランダムフォレストによる株価の予測

Pythonで機械学習を使った株価予測のコードを書こう
こちらでは、ランダムフォレストを用いた機械学習で、ETFなどのデータを特徴量にして個別の株価予測を行っています。

アカデミックサイドでの研究事例

Stock Price Forecasting Using Information from Yahoo Finance and Google Trend
UC Berkeleyの経済学部生の研究です。こちらはR言語で、ヤフーファイナンスとGoogleトレンドの情報を用いて株価を予測している研究です。これまでの時系列手法よりもパフォーマンスが良いそうです。

最後に

事例を集めてみて、どこまでのレベルのものに手を出すべきか悩ましいと思いました。完全自動化は維持するコストがかかりそうですし、ロスが増大しないか心配です。自分としては情報を自動で取得し、リターンが発生する確率の高そうな銘柄をサジェストしてくれるレベルで十分な気がします。

OpenCV&Pythonで画像の類似度を計算させる〜イケメンの顔比較

・動機
・やりたいこと
・準備
・類似度の計算
・実行コード
・実行結果
・おまけ

動機

画像系の技術にあまり関心が無かったのですが、とある知人が福士蒼汰のような雰囲気の男性が好みであると発言されたことを発端に、福士蒼汰に最も顔の近い知人を見つけるというプライベートなミッションを仰せつかりました。
そこで、まずは最も楽だろうと思われる、画像間の類似度を計算する方法について調べてみました。顔のパーツを検知して、パーツ同士で比較するなどのレベルでは無いことをご了承下さい。ちなみに、比較画像は国民的アイドルである嵐のメンバーの画像としました。今回の実践で、与えられた画像の中で、福士蒼汰に最も近そうな嵐のメンバーの写真がわかることになります。

やりたいこと

画像間の類似度の計算

準備

・Python2.7
・OpenCV

OSX環境における準備にあたっては以下の情報が参考になりました。
PythonでOpenCVを使う
Mac OS X で OpenCV 3 + Python 2/3 の開発環境を整備する方法

類似度の計算

こちらのブログにある計算手法とコードを使いました。
How-To: Python Compare Two Images

紹介されているコードで以下の評価指標が計算できます。
・Mean Squared Error (MSE)
・Structural Similarity Index (SSIM)・・・0〜1の値を取ります。

指標について、詳しくはこちらの論文に書かれています。
Image Quality Assessment: From Error Visibility to Structural Similarity

実行コード

紹介されていたコードは画像サイズが同じでないと計算ができなかったので、まずは画像のサイズを整えるためのコードが以下のようになります。コードはこちらのものを使いました。
Python – アンチエイリアスで写真をキレイに縮小

これだけではサイズが同一にならなかったので、追加で以下のコードを実行しました。

以下は、類似度計算の実行用コードとなります。

実行結果

福士蒼汰との比較をしており、画像上部にMSEとSSIMが出ています。嵐の大野くんが一番近いようです。

original

aiba

oono

sakurai

ninomiya

matsumoto

おまけ

イケメンではないですが、出川哲朗でも計算してみました。

degawa

今回最も高い値が出てしまいました。やはり、顔のパーツを識別して、そのパーツ間の類似度の計算ができないといけないような気がします。悔しいので今後も画像認識系の技術について向き合ってみようと思います。

scikit-learnのモジュールのGitHubでの利用頻度を調べてみた

『Python機械学習プログラミング』を読んで、scikit-learnのモジュールは充実しているなぁと感じたんですが、実際にWebサイトでUser Guide( http://scikit-learn.org/stable/user_guide.html )を見た所、この本に載り切らないような数多くの機械学習手法に応じたモジュールが用意されていました。そこで、世のデータサイエンティストはどのモジュールを良く使っているのだろうと気になったので、GitHubのSearchでヒットしたCodeの数を各モジュール単位で集計してみました。検索クエリは「scikit-learn + モジュール名」なので、正確なものではないのですが、相対的な利用頻度を見るぶんには使えるのではないかと思われます。

データ集計方法

・User Guideに登場するscikit-learnのモジュール名を集めています。
・教師付き学習か教師無し学習かどうかの判断は、User Guideで紹介されているモジュールかどうかで判断しています。
・GitHubのSearchで「scikit-learn + モジュール名」でヒットした件数をそのまま使っています。(2016年9月22日時点)

可視化コード

Jupyterで実行しています。

教師付き学習編

圧倒的に多いのがSVM(Support Vector Machine)を扱っているSVCモジュールで、続いて定番のロジスティック回帰やRandom Forestが使われているようです。統計解析ではメジャーなはずの線形回帰が5位なのは、初歩的なのであまりコードがアップされていないのかもしれません。GBDTのモジュールももう少し上位にくるかと思ったんですが15位でした。DMLCのXGBoostモジュールを使っているのかもしれませんね。私も実際のところXGBoostを使ってますし。

supervised_module_barplot

教師無し学習編

主成分分析やK-mean法など因子分解などのモジュールが上位を占めています。LDA(Latent Dirichlet Allocation)がもっと上位に来ると思ったんですが、思えばGensimの方が充実しているなぁと思うので、このランキングは妥当なのかもしれません。私もLDAなどはGensimを使っていますし。

unsupervised_module_barplot

収集を終えて

・社内だとデータサイエンティストの方がいないので、scikit-learnのモジュールの利用状況を知れてマニアックな共感をすることができた。
・SVMは実践例が豊富そうなので分析事例を探せば良い発見があるかもしれない。
・scikit-learnのUser Guideは充実していたので、時間を作って向き合ってみたいと思った。

XGBoostのパラメータチューニング実践 with Python

以前の投稿で紹介したXGBoostのパラメータチューニング方法ですが、実際のデータセットに対して実行するためのプログラムを実践してみようと思います。プログラム自体はAnalytics_Vidhya/Articles/Parameter_Tuning_XGBoost_with_Example/XGBoost models.ipynbに載っているのですが、データセットがついていません。そこで、前回の投稿(不均衡なデータの分類問題について with Python)で赤ワインのデータセットを手に入れているので、こちらのデータセットを用います。誤植なのかところどころ、うまく回らなかったところがあったので、手直しをしています。

以下の工程に従って進みます。結構長いですが、辛抱強く実践してみて下さい。
・ライブラリの読み込み
・データの読み込み
・前処理
・学習用データとテスト用データの作成
・XGBoostの予測結果をもとに、AUCの数値を返すための関数の定義
・モデルの実行
・チューニング

ライブラリの読み込み

データの読み込み

前処理

スクリーンショット 2016-05-15 16.54.58

学習用データとテスト用データの作成

XGBoostの予測結果をもとに、AUCの数値を返すための関数の定義

XGBoostの予測結果から、AUCの数値を返し、特徴量に応じた重要度を出力するためのプログラムです。

モデルの実行

feature_importance1

チューニング

max_depthとmin_child_weightの数値をチューニングするためのプログラムです。

より細かい数値で再度最適なパラメータを探します。

max_depthを8、min_child_weightを1として、他のパラメータチューニングに移ります。
続いて、gammaのチューニングを行います。

gammaを0.4と置きます。
ここで、いままでにチューニングしたパラメータを用いて再度推定を行います。先ほどの0.875よりも高くなっています。

feature_importance2

続いて、subsampleとcolsample_bytreeのチューニングを行います。

より細かい範囲で再度パラメータをチューニングします。

続いて、reg_alphaをチューニングします。

範囲が粗かったので、より細かくパラメータをチューニングします。

これまでにチューニングしてきたパラメータを用いて再度推定を行います。

feature_importance3

ブログであるように試行回数を1,000回から5,000回まで増やしてみます。

88.8%まで向上しました。色々と数値いじっても、1%高めるだけにとどまってしまうのですね。

feature_importance4

とにかく、XGBoostをPythonで実行してパラメータチューニングするという一連の試行がこのコードでできそうなので、今後も使いまわしてみようと思います。

不均衡なデータの分類問題について with Python

データマイニング界隈で人気のKDnuggetsで紹介されていた、”Dealing with Unbalanced Classes, SVMs, Random Forests, and Decision Trees in Python“のプログラムが残念なことに画像だったので、写経しました。せっかくなので、紹介させていただきます。内容としては不均衡データに対する処方の紹介で、プログラムはPythonで書かれています。ライブラリさえインストールできれば皆さんもすぐに実行できるので、是非チャレンジしてみて下さい。

まずはもろもろライブラリを呼び出します。

CSV形式のデータセットをWebサイトから取得します。ワインの評価と、ワインに関した特徴量からなるデータセットです。

スクリーンショット 2016-05-09 00.24.50

分類のための目的変数を作成します。

スクリーンショット 2016-05-09 00.26.24

ランダムフォレストを実行します。

どんな結果が返ってくるのか、試しに一つだけツリーの数を2にして実行してみます。10回分のクロスバリデーションを行った推定結果が出力されています。これは、いわゆる正解率のことを指します。

ツリーの数に応じた正解率を可視化します。

正解率はツリーの数を増やすことで増すようです。

Classification score for number of trees

しかしながら、正解率は誤解されやすい指標です。不均衡データでは偏りのある方ばかりを当てていても、正解率は増してしまいます。当たりだけを予測できて、ハズレを予測できないというのは分類器として使いみちが限られると思います。そこで、悪いワインの割合を直線で引いてみます。

悪いワインの割合がそもそも多いので、悪いワインと判定しまくっていても、正解率は高いわけです。

Classification score for number of trees2

そこで、機械学習における予測精度の評価指標とされているF値を使います。
ツリーの数を増やしても、F値は良くなっていないようです。

F1 scores as a function of the number of trees

ここでは、0.5よりも大きいとする予測になる特徴量のデータを切り捨てます。その切り捨てる割合がどこが望ましいのかを以下で探していきます。

どうやら、階級値が2~4、つまり割合にして0.3~0.5のカットオフ値が望ましい水準のようです。

custom F scores

以下では、決定境界の可視化を行います。しかしながら、二次元の可視化となると、複数あるデータの中から特徴量を二つだけ選ばなければなりません。その変数を決めるに際して、変数の重要度を用います。変数の重要度はランダムフォレストで計算可能です。

Importance of each feature

こちらで、重要度が上位のものに絞って決定境界を可視化します。ここでは、ランダムフォレストのみならず、SVMや決定木も実行されています。

Decision Tree Classifier

Random Forest Classifier

Support Vector Maachine

sklearnのSVMはデフォルトではクラスごとの重み付けを行わないが、自動で重み付けを行うことが出来る。以下の例では、C=1、gamma=1でクラスごとの重み付けを行う・行わないでの決定境界を描いている。重み付けを行うことで、赤色の少ない方のデータの識別が比較的できていることが伺えるが、他方で、多くの青を誤判定している。さらなる改善にはパラメータチューニングが必要となります。

Svm without class weight

Svm with class weight

不均衡データに対するアプローチや、Pythonによる機械学習を学ぶ良い機会になりました。KDnuggetsは非常に勉強になりますね。

XGBoostやパラメータチューニングの仕方に関する調査

【目次】
・XGBoostとは
・XGBoostで用いるパラメータ一覧
・XGBoostのパラメータチューニング
・参考文献

XGBoostとは

XGBoost (eXtreme Gradient Boosting) は勾配ブースティングアルゴリズムの先進的な実装例で、データサイエンスのコンペであるKaggleで話題となっていた手法です。

ブースティングアルゴリズムとは、弱識別器(weak learners)の集団を直列的に結合することで、予測における正確性を高めようとするアルゴリズムです。任意のt時点において、モデルの予測は以前のt-1時点での結果に基づき重み付けがなされます。正しく予測されたデータに対しては、重みを小さくし、誤って予測されたデータに対しては重みを大きくします。後で学習する識別器ほど、誤ったデータに集中して学習を進めることになります。

以下はブースティングのイメージ図です。
スクリーンショット 2016-04-24 18.40.38

STEP1では全ての学習データに対して、等しい重み付けで学習を行い、決定境界を引きます。これを弱学習器による学習と言います。このケースでは毒キノコを2つ当てており、キノコを5つ当てています。
STEP2ではSTEP1で正しく識別されたデータの重みが下げられ、誤って識別されたデータの重みが上げられています。高く重み付けがなされたデータは決定境界で正しく識別されていますが、他のデータは誤って分類されています。
STEP3においてもSTEP2と同様の傾向があります。このような弱学習器による処理を繰り返すことで識別性能を高めていきます。

最終的にはこのような決定境界を引くことができるような識別器を求めていきます。
スクリーンショット 2016-04-24 18.55.53

勾配ブースティングの勾配とは、ブースティングアルゴリズムにおけるパラメータ推定での最適化手法が勾配降下法に従っているという意味での勾配です。以上が勾配ブースティングモデルの簡素な説明です。

XGBoostで用いるパラメータ一覧

XGBoostで用いるパラメータに関して、大きく分けて3つあります。

    1.全体パラメータ・・・XGBoost全体を司る。
    2.ブースターパラメータ・・・各ステップでツリーなどのブースティングを司る。
    3.学習タスクパラメータ・・・最適化タスクを司る。

以下、3つのパラメータについて、「パラメータ名」・「デフォルトの値」・「役割」・「引数」を表にしています。

1.全体パラメータ

パラメータ名 デフォルトの値 役割 引数
booster gbtree 実行するモデルのタイプをツリーモデルか線形モデルのどちらかを指定できる。 gbtree: ツリーモデル
gblinear: 線形モデル
silent 0 モデルの実行結果を出力するかどうかを決めることができる。モデルを理解する上で、0のままにしておく方が良いとされている。 0:出力結果を表示する。
1:出力結果を表示しない。
nthread not set 並列処理のためのコア数などを指定できる。フルコアで実行したい場合は何も指定しなければ自動的にフルコアになる。

2.ブースターパラメータ

パラメータ名 デフォルトの値 役割 引数
eta 0.3 学習率を調整できる。
小さくすることで、モデルの頑健性を高めることができる。
0.01〜0.2の値になることが多いらしい。
min_child_weigh 1 子ノードにおいて観察されるデータの重み付けの合計値の最小の値で、過学習を避けるために用いられる。
高い値にすることで特定のサンプルに観察されるような傾向を学習することを避けられる。ただし、高くし過ぎるとフィッティングが悪くなる。
max_depth 6 木の深さの最大値
過学習を制御するために用いられる。
高いと過学習しやすくなる。
3〜10の値になることが多いらしい。
max_leaf_nodes 木の終端ノードの最大値
max_depthの代わりに用いる
n本を指定したら、n^2個の枝を生み出す。これが指定された場合は、max_depthは無効化される。
gamma 0 分割が、損失関数の減少に繋がる場合にのみノードの分割を行う。
モデルをより保守的にする。
値は損失関数に応じて大きく変わり、チューニングが必要である。
max_delta_step 0 各木のウェイトの推定に制約をかけることができる。
0の場合は制約なしで、正数値を取るとモデルがより保守的になる。
通常は必要とされないが、不均衡データの分類の際に用いる。
subsample 1 各木においてランダムに抽出される標本の割合
小さくすることで、過学習を避けることができるが保守的なモデルとなる。
0.5〜1の値になることが多いらしい。
colsample_bytree 1 各木においてランダムに抽出される列の割合 0.5〜1の値になることが多いらしい。
colsample_bylevel 1 各レベル単位での、分割における列のsubsample比率
subsampleとcolsample_bytreeで十分なので、あまり使わないが、探索してみるのも良いかもしれない。
lambda 1 重みに関するL2正則化項
多くのデータサイエンティストは使わないが、過学習を避けるためには用いられるべき。
alpha 0 重みに関するL1正則化項
高次元の場合に用いるらしい。
scale_pos_weight 1 不均衡データの際に、0以上の値を取ることで、収束を早めることができる。

3.学習タスクパラメータ

パラメータ名 デフォルトの値 役割 引数
objective reg:linear 最小化させるべき損失関数を指定する。 binary:logistic→2項分類で確率を返す。
multi:softmax→多項分類でクラスの値を返す。
(num_classでクラス数の指定が必要)
multi:softprob→softmaxと同じだが、確率を返す。
eval_metric according to objective 検証を行うためのデータの評価指標 rmse – root mean square error
mae – mean absolute error
logloss – negative log-likelihood
error – Binary classification error rate (0.5 threshold)
merror – Multiclass classification error rate
mlogloss – Multiclass logloss
auc: Area under the curve
seed 0 ランダムなシード番号。
再現可能なデータを生み出すために、あるいはパラメータチューニングの際に用いる。

 

 

XGBoostのパラメータチューニング

複数のパラメータからなるXGBoostのチューニングは非常に複雑で、理想的なパラメータについてはケースバイケースで何とも言えないそうです。

参考文献のブログにパラメータチューニングの一般的アプローチについて触れられていたので、紹介します。

    1.相対的に高い学習率を選択する。
    一般的な学習率0.1はうまくいくが、ケースによっては0.05から0.3の間でもうまくいく。学習率を求めるために、クロスバリデーションによって最適な木の数を決定する。
    2.木に関するパラメータをチューニングする。
    固定された学習率や木の数のもとで、max_depth, min_child_weight, gamma, subsample, colsample_bytreeをチューニングする。
    3.正則化パラメータをチューニングする。
    lambda, alphaなどのパラメータをチューニングすることで、モデルの複雑さを減らし、パフォーマンスを高める。
    4.学習率を下げ、最適なパラメータを決定する。

具体的な実行に関するPythonスクリプトはこちらのGithubで紹介されています。(iPython)
Analytics_Vidhya/Articles/Parameter_Tuning_XGBoost_with_Example/XGBoost models.ipynb
この方法に従って、自社で抱えているモデルのチューニングにチャレンジしてみようと思います。

参考文献

Quick Introduction to Boosting Algorithms in Machine Learning
Complete Guide to Parameter Tuning in Gradient Boosting (GBM) in Python
Complete Guide to Parameter Tuning in XGBoost (with codes in Python)
勾配ブースティングについてざっくりと説明する
xgboost のパラメータ
OS X で XGBoost & xgboost4j をビルドする手順 2016-03-07 版

GensimのHDP(Hierarchical Dirichlet Process)をクラシック音楽情報に対して試してみる

HDP(Hierarchical Dirichlet Process)いわゆる階層ディリクレ過程を実行できるモデルがPythonのGensimライブラリにあるという情報から、あまり実行例も見当たらないので、チャレンジしてみました。

HDP(Hierarchical Dirichlet Process)

HDP(Hierarchical Dirichlet Process)は文書集合全体のトピック数と文書ごとのトピック数の推定を行うことができる手法で、中華料理店フランチャイズという仕組みを用いています。通常のLDAなどでは、分析者が任意のトピック数を決める必要がありましたが、与えられたデータからその数を推定するため、その必要がないというのがHDPを使うことの利点であると思われます。

実行までの流れ

ざっくりですが、
・コーパスの準備・文書の分かち書き(名詞のみ)
・HDPの実行
という流れです。

ちなみに実行環境は
MacBook Pro
OS X Yosemite 10.10.5
2.6 GHz Intel Core i5
メモリ8GBです。

コーパスの準備

今回は、以前手に入れた某辞典サイトのクラシック音楽情報1800件のテキストデータ(1行に1件分の文字列が入っているデータで16MBくらい)があるので、それをコーパスとして使います。参考情報として挙げているブログの助けを借りて、文書単位でMeCabにより形態素解析で分かち書きした結果から、意味を持ちやすい品詞として、「名詞」に該当するもののみを結果として返す以下のPythonスクリプトを用いました。結果はtmep.txtとして出力されます。もっと良いやり方があると思いますが、目的は達成できると思います。ちなみに、MeCab Neologd(ネオログディー)という、固有名詞などに強いシステム辞書を活用してみたかったので、その利用を前提として書いています。MeCab Neologd(ネオログディー)のインストール関連の情報は参考情報にありますので、チャレンジしてみてください。(OSXかUbuntuの方が進めやすいと思います。)

こちらのスクリプトをターミナルで実行します。(解析するディレクトリ下で実施しています。)

HDPの実行

以下のPythonスクリプトで実行しています。

HDPの結果について

topic_detail.csvの結果を見たところ、トピックの数が150個もあって、「本当にトピックの数を自動で決めれているのかなぁ」と不安に思ったのですが、実際に各文書に割り当てられているトピックの数は、先ほど出力したtopic_for_corpus.csvで見ると60個でした。そのため、今回、HDPに従って決まったトピック数は60ということになります。さらに不安に思ったので、Stack Over Flowで調べていたんですが、トピックは150個出るけど確率が割り振られていないはずと回答されていました。( Hierarchical Dirichlet Process Gensim topic number independent of corpus size

出現頻度の高い上位10のトピックは以下の通りです。

加えて、トピックごとに文書に割り当てられた数を集計してみましたが、topic0が圧倒的に多く、コーパスの特性上、含まれやすい情報がここに集まっているのではないかと思います。幅広いテーマを抽出できるかと期待していたのですが、やたらと個別具体的な「トゥーランドット」や「ワーグナー」や「カルメン」などがトピックの上位単語に上がってきています。実行方法を間違えているかもしれないし、パラメータチューニングなどをもっと頑張れば、幅広いトピックを得ることができるかもしれないので、今後の課題としたいです。

スクリーンショット 2016-04-10 18.23.37

参考情報

・トピックモデルについて
machine_learning_python/topic.md at master · poiuiop/machine_learning_python · GitHub

・HDP関連
 models.hdpmodel – Hierarchical Dirichlet Process
Online Variational Inference for the Hierarchical Dirichlet Process

・MeCab関連
mecab-ipadic-NEologd : Neologism dictionary for MeCab
形態素解析器 MeCab の新語・固有表現辞書 mecab-ipadic-NEologd のご紹介
テキストマイニングの前処理。名詞抽出、ストップワード除去、珍しい単語の除去

・Python関連
データ分析をやりたいエンジニアにおすすめ!Pythonの入門スライド13選

Word2Vecでクラシックの楽曲情報をコーパスとして類似度を出してみる

あの手この手を使って手に入れた、およそ1800曲に及ぶクラシック音楽の楽曲情報(テキスト)をもとに
、PythonのGensimライブラリーのWord2Vecを使って、任意の単語に関する類似単語を出力してみたいと思います。

まずは、手に入れたコーパスを作業フォルダに置いて、MeCabによる分かち書きを行います。
(最後の引数-bは、処理する文書のサイズが大きい際に調整します。)

後はGensimパッケージを読み込んで、

Word2Vecを計算させるだけです。

(引数のsizeは特徴ベクトルの次元数です。)

早速、トランペットについて、所与のコーパスにおける類似単語を見てみたいと思います。(類似度が最も高い単語の上位10位の結果を返しています。)

惜しいですね。
願わくば、コルネットが一番目に来てほしかったです。オケの編成上、どうしてもトロンボーンが一緒の文書で出やすいのだと思います。

続いては、ピアノです。こちらもオルガン・チェンバロは非常に近い楽器だと思うのですが、一番目がヴァイオリンというのはデータ上仕方がないのかもしれません。

続いて、ヴァイオリンですが、ヴィオラ・チェロは良いと思うのですが、ピアノやチェンバロなどが上位に来ています。

続いて、クレッシェンドですが、似たような意味はあまり観察されていません。ただし、「クライマックス」・「アルペッジョ」などと似たようなシチュエーションで登場しそうな表現な気がします。

最後に、アレグロですが、こちらは速さの序列に関しては守られていないようです。やはりコーパス次第ですかね。
プレスト > アレグレット > モデラート > アンダンテ > アダージョ > ラルゴ
この序列が守られるようなWord2Vecの実践例などがあると面白いですが。

仕事でWord2Vecを使うシーンがあるとしたら、広告文のアイデアを助けたり、語彙力の弱い人の補助的なツールとして使えるかもしれませんが、実用レベルはまだまだ遠い気がします。

おまけ

左手に関しては、最も類似した単語が「右手」という結果になっています。

参考文献

models.word2vec – Deep learning with word2vec

岩波データサイエンス Vol.2 岩波データサイエンス刊行委員会 編

iPython notebookチャレンジ(カイ二乗検定)

IPython データサイエンスクックブック』の第7章のレシピ7.4を写経して、jupyterのdownload asでhtmlをはき出し、そのままブログに貼り付けてみました。色とかおかしいところがあるので、出力方法に関しては今後色々と調べなくてはなりません。ただ、実際に触ってみて、ユーザーインターフェースに関してRとの著しい乖離はない印象です。ここでは扱いませんが、ランダムフォレストを用いた重要な回帰特徴量などもRに比べてコードの長さも遜色ありませんでした。iPython jupyterはRからPythonへ手を広げるハードルを下げてくれる良い環境だと思いました。
以下では、Datasets used in the IPython Cookbook, the definitive guide to high-performance scientific computing and data science in Pythonのtennisという、テニス大会と選手に関するデータをもとに、ロジャー・フェデラーのエースと得点との関係を扱っています。手法としては相関係数の計算と、カイ二乗検定を行っています。








In [1]:
#NumPy,pandas,SciPy.stats,matplotlibをインポートする。
import numpy as np
import pandas as pd
import scipy.stats as st
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
#ロジャー・フェデラーのデータを読み込む。カレントディレクトリにデータを置いておきましょう。
player = 'Roger Federer'
filename = "data/{name}.csv".format(
                       name=player.replace(' ', '-'))
df = pd.read_csv(filename)
In [31]:
#データを確認する。
print("Number of columns: " + str(len(df.columns)))
df[df.columns[:8]].tail()
Number of columns: 70
Out[31]:
year tournament start date type surface draw atp points atp ranking
1174 2012 Australian Open, Australia 16.01.2012 GS Outdoor: Hard Draw: 128 720 3
1175 2012 Doha, Qatar 02.01.2012 250 Outdoor: Hard Draw: 32 90 3
1176 2012 Doha, Qatar 02.01.2012 250 Outdoor: Hard Draw: 32 90 3
1177 2012 Doha, Qatar 02.01.2012 250 Outdoor: Hard Draw: 32 90 3
1178 2012 Doha, Qatar 02.01.2012 250 Outdoor: Hard Draw: 32 90 3
In [4]:
#全ポイントに占める特典の割合とエースの割合だけを表示する。
npoints = df['player1 total points total']
points = df['player1 total points won'] / npoints
aces = df['player1 aces'] / npoints
In [5]:
plt.plot(points, aces, '.')
plt.xlabel('% of points won')
plt.ylabel('% of aces')
plt.xlim(0., 1.)
plt.ylim(0.)
Out[5]:
(0.0, 0.16)
In [6]:
#二つの列だけを持つDataFrameオブジェクトを新しく作る。
df_bis = pd.DataFrame({'points': points,
                      'aces': aces}).dropna()
df_bis.tail()
Out[6]:
aces points
1173 0.024390 0.585366
1174 0.039855 0.471014
1175 0.046512 0.639535
1176 0.020202 0.606061
1177 0.069364 0.531792
In [7]:
#ピアソン相関係数を計算する。
df_bis.corr()
Out[7]:
aces points
aces 1.000000 0.255457
points 0.255457 1.000000
In [9]:
#変数の二値化
df_bis['result'] = df_bis['points'] > df_bis['points'].median()
df_bis['manyaces'] = df_bis['aces'] > df_bis['aces'].median()
In [10]:
#それぞれの可能性の頻度からなる分割表を作る。
pd.crosstab(df_bis['result'], df_bis['manyaces'])
Out[10]:
manyaces False True
result
False 300 214
True 214 299
In [11]:
#カイ二乗検定の統計値とp値の計算(2番目の数値がp値)
st.chi2_contingency(_)
Out[11]:
(27.809858855369555,
 1.3384233799633629e-07,
 1,
 array([[ 257.25024343,  256.74975657],
        [ 256.74975657,  256.25024343]]))
In [ ]:
 

Sparkの仕組み・導入メリット・インストール方法など

【目次】
・Sparkについて
・Sparkの歴史
・Sparkの仕組み
・データ分析におけるSpark導入のメリット
・Sparkの導入方法
・参考文献

Sparkについて

OSS(Open Source Software)のインメモリ処理の分散並列基盤
※インメモリ(メインメモリ上で処理するため、ディスクストレージを介するよりも高速かつ安定しているらしい。メモリに乗り切らないようなケースでも、ディスクを利用するなどしてアプリケーションが問題なく動作するように作られている。)
ビッグデータ活用を支えるプラットフォームとして注目されている。
比較の際はHadoopで言う所の、MapReduceが引き合いに出されます。MapReduceがバッチ処理に特化した仕組みになっているのに対して、Sparkはバッチ処理だけでなく、対話的な処理や繰り返し処理・ストリーミング処理といった多様なデータ処理のパターンに対応できる仕組みになっていることから、Sparkに注目が集まっています。

分散処理エンジンとして、Apache Spark Coreがあり、以下の様な機能が用意されているようです。

Spark SQL・・・SQL処理用
Spark Streaming・・・ストリーミング処理
MLib・・・機械学習処理用
 基本API
 ・統計情報の取得
 ・相関関係の算出
 ・層化抽出
 ・仮説検定
 ・ランダムデータ生成
 機械学習アルゴリズム
 ・分類回帰(SVM、ロジット、線形回帰、ナイーブベイズ、決定木、ランダムフォレスト、isotonic回帰)
 ・協調フィルタリング(ALS(Alternating Least Squares))
 ・クラスタリング(K平均法、混合ガウシアン、PIC、LDA、ストリーミングK平均法)
 ・次元削減(SVD、PCA)
 ・特徴抽出&変換(TF-IDF、Word2Vec、etc…)
 ・頻出パターンマイニング(FP-Groth、アソシエーションルール、PrefixSpan)
GraphX・・・グラフの処理
※機械学習の機能があるのは嬉しいですね。スクラッチで作るのはペインフルなので。

データソースとしては、HDFS(Hadoop Distributed File System)だけでなく、Hive、HBase、PodgreSQL、MySQL、CSVファイルなども対応している。

Sparkの歴史

2009年にUC BarkeleyのMatei Zaharia氏がScalaを用いて開発

Google Trendで検索トレンドを見たところ、2014年頃から今に至るまで伸び続けています。
GoogleTrend

Sparkの仕組み

RDD(Resilient Distributed Dataset)
不変で並列実行可能な分割されたコレクションのことを指します。
これによりユーザーが並列分散処理を意識しなくとも、それができる環境が提供されます。
こちらのRDDを
map,filter,union,flatMapなどのTransformations
で生み出し、RDDのデータを
reduce,countなどのActions
を実行することで動かします。
RDDを任意の時点でメモリやディスクに保存することもできるため、処理の再実行のみならず、同じデータを複数回利用するなど、複雑なデータ処理に関しても効率的に実行できます。
SparkはScalaという言語で書かれていますが、Python・Java・Rなどに対してAPIが提供されているので、幅広いユーザーが扱うことが可能です。しかしながら、痒いところに手を届かせるためには、Scalaの勉強は避けられないのかもしれません。実際のところ、大体のデータマイニング系エンジニアの知り合いはScalaの勉強などをしていた印象があります。

データ分析におけるSpark導入のメリット

・並列分散の手法を意識せずに大規模データの処理が可能になる。
・Rなどの苦手とする反復処理を高速化させることができる。
・データ取得、前処理、分析処理、レポーティングを一つのプラットフォームで実行できる。
・コマンドラインによるインタラクティブな試行錯誤がしやすいので、データ分析のタスクに向いている。(PythonやRのように)

※レコメンドエンジンとしてSparkを用いるケースも観察されているようです。
 Yahoo!台湾

Sparkの導入方法

ここでの導入はあくまでもプログラミングを学ぶためであって、ローカル環境での実行環境の構築ではSparkのメリットを感じることはできません。
※SparkはJava仮想マシン上で動作するため、事前にOracle JDKをインストールしておきましょう。
Java SE Development Kit 8 Downloads

■Windowsの場合
Apache Sparkで始めるお手軽リアルタイムウインドウ集計
こちらで手順が詳しく書かれています。

■Macの場合
Macの方が導入はシンプルです。
Download Spark
からダウンロードします。(SparkRはSparkのver1.4以降からサポートされているそうです。)
HOW TO INSTALL APACHE SPARK ON MAC OS X YOSEMITE
Apache Spark 1.6.0 setup on Mac OS X Yosemite
そこで、

build/mvn -Pyarn -Phadoop-2.4 -Dhadoop.version=2.4.0 -DskipTests clean package

を実行すればコンパイルできるようです。(10分ほどかかる。私は30分もかかったが。。)

そのあと、パスを通します。
PATH=/usr/local/share/spark/bin:$PATH; export PATH

それから、pythonでsparkを動かしてみます。
pyspark

を実行して、試しにReadmeのテキスト行数をカウントさせるコマンドを入力します。これは、textFileメソッドを用いて、Readmeファイルを読み込み、linesという名前のRDDを生成していることを意味します。
lines = sc.textFile(“/usr/local/share/spark/README.md”)

そして、このコマンドで値が返されたら、うまく動いていることが確認できるでしょう。
print lines.count()

加えて、ワードカウントを行うコマンドを入力します。(スペース区切りで分解し、mapメソッドで(単語, 1)という組み合わせに変換したあと、最終的にreduceByKeyで同一の単語を足し上げるという処理をしています。)

counts = lines.flatMap(lambda line: line.split(” “)) \
.map(lambda word: (word, 1)) \
.reduceByKey(lambda a, b: a + b)

そして、これで結果を取得します。
counts.collect()

スクリーンショット 2016-05-02 12.21.50

参考文献

Spark を活用する:ビッグデータアプリケーション用の高速インメモリコンピューティング

Apache Spark入門 動かして学ぶ最新並列分散処理フレームワーク (NEXT ONE)

詳解 Apache Spark

[Apache Spark]RDDについて簡単にまとめてみた

Rユーザのためのspark入門

Spark 1.4.0 の SparkR を動かす

Spark Examples