データ分析がしたい

企業でデータ分析などやっています。主にRやPythonによるデータマイニング・機械学習関連の話題やその他備忘録について書いてます。

scikit-learnを用いたサンプルデータ生成

機械学習の勉強や新しいアルゴリズムのテストをする場合、irisなどのシステム組み込みのサンプルデータを利用するか、UCIリポジトリなどのネット上の公開データから良さげなものを探すというのが一般的だと思います。

しかしながら、irisなどの組み込みデータは一般にデータ数が少なく、分類問題として物足りなかったり、ネット上の公開データを利用するにしても適当なデータ数や特徴量数、問題設定や難度のデータを探すのが難しいですし、前処理が必要なデータも多く手軽に使えるサンプルデータとなると中々見つけられないといったことがあるかと思います。

そういった場合、適当なデータ数や難しさのデータを自分で生成して利用すると、後の計算コスト評価や機械学習アルゴリズムの理解において色々と便利です。
サンプルデータの作り方としては、何らかの統計モデルに基づいて作る方法もありますが、データの質にこだわらないのであればscikit-learnのサンプルジェネレータを使って作るのが手っ取り早いです。

本記事では特にscikit-learnのmake_classificationを利用した分類問題用データの作成方法について紹介しますが、scikit-learnには分類問題以外にも回帰やマルチラベル分類用のデータや疎なデータを生成する関数もあるので、興味があれば調べてみると良いと思います。(http://scikit-learn.org/stable/index.html

make_classificationの使い方

make_classification関数はsklearn.datasetsに含まれており、サンプル数や特徴量数などいくつかのパラメータを指定するだけで、簡単に分類問題用のサンプルデータ作成することができます。
関数定義は下記の通りです。

sklearn.datasets.make_classification(
    n_samples=100,n_features=20, n_informative=2, n_redundant=2,
    n_repeated=0, n_classes=2, n_clusters_per_class=2, weights=None,
    flip_y=0.01, class_sep=1.0, hypercube=True, shift=0.0, scale=1.0,
    shuffle=True, random_state=None)

パラメータは色々ありますが、例えば、データ数が1000で特徴量が20の2値分類データを作る場合、n_sample=1000、n_feature=20、n_class=2とすると(1000, 20)と(1000,)のnummpy.arrayオブジェクトを返します。そして、n_informative、n_redundant、n_clusters_per_class、flip_yなどでデータ構造や分類の難しさをコントロールする感じです。

下記に各パラメータの大まかな説明を記載しています。

パラメータ名 説明
n_samples 生成するサンプルの数。
n_features 生成する特徴量の数。
n_informative 目的変数のラベルと相関が強い特徴量(Informative fearture)の数。
n_redundant Informative featureの線形結合から作られる特徴量(Redundant fearture)の数。
n_repeated Infomative、Redundant featureのコピーからなる特長量の数(Repeated feature)。
n_classes 分類するクラス数。2なら2値分類問題、3以上なら多値分類問題のデータが作られる。
n_clusters_per_class 1クラスあたりのクラスタ数。後述する生成アルゴリズムを参照。
weights クラスの比率。不均衡データを作りたい場合に指定する。例えば、2値分類問題の場合、Noneとすると0と1が50%ずつだが、[0.9, 0.1] と与えると0が90%、1が10%になる。
flip_y クラスのフリップ率。例えば0.01とすると各クラスの1%の符号がランダムに変更される。
class_sep 生成アルゴリズムに関係するパラメータ。後述する生成アルゴリズムにおける超立法体の頂点の距離。
hypercube 生成アルゴリズムに関係するパラメータ。Trueにすると後述する生成アルゴリズムにおいて、超立方体の頂点にクラスタを配置する。
shift 全ての特徴量にshiftを加算する。Noneが指定された場合、[-class_sep, class_sep]の一様乱数を加算する。
scale 全ての特徴量にscaleを乗算する。Noneが指定された場合、 [1, 100]の一様乱数を乗算する。(scaleはshift後に行われる)
shuffle Trueにすると行と列をシャッフルする。
random_state 乱数を制御するパラメータ。Noneにすると毎回違うデータが生成されが、整数をシードとして渡すと毎回同じデータが生成される。乱数オブジェクトを渡すことも可能。

(参照:sklearn.datasets.make_classification — scikit-learn 0.16.1 documentation


具体的な使い方は次のコードのような感じです。

### classification sample
import numpy as np
from sklearn.datasets import make_classification
from sklearn.cross_validation import train_test_split

from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import roc_auc_score

# サンプルデータの生成
# 1000 samples、5(infomative) + 2(redundant) + 13(independent) =  20 feature のデータを生成
dat = make_classification(n_samples=1000, n_features=20, n_informative=5, n_redundant=2,
                            , n_classes=2, n_clusters_per_class=10)

X = dat[0]
y = dat[1]
print "X shape", X.shape
print "y shape", y.shape

# 学習用とテスト用データの分割
# 80%を学習、20%をテストに利用する
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)

# 学習モデルの構築とパフォーマンス評価
# ロジスティック回帰、ランダムフォレスト、KNNの3つのモデルを作成しそれぞれのAUCを計算
clf = LogisticRegression()
clf.fit(X_train, y_train)
print "LogisticRegression AUC =", roc_auc_score(y_test, clf.predict_proba(X_test)[:,1])

clf = RandomForestClassifier(n_estimators=500, random_state=123)
clf.fit(X_train, y_train)
print "RandomForestClassifier AUC =", roc_auc_score(y_test, clf.predict_proba(X_test)[:,1])

clf = KNeighborsClassifier(n_neighbors=10)
clf.fit(X_train, y_train)
print "KNeighborsClassifier AUC =", roc_auc_score(y_test, clf.predict_proba(X_test)[:,1])

結果

X shape (1000, 20)
y shape (1000,)
LogisticRegression AUC = 0.711838942308
RandomForestClassifier AUC = 0.802534054487
KNeighborsClassifier AUC = 0.793820112179

上記コードではmake_classificationで1000サンプル20特徴量のデータを生成し、train_test_splitで学習とテスト用データに分割し、ロジスティック回帰とランダムフォレスト、KNNの3種類の分類機に入れてAUCを出力しています。(random_stateは毎回同じデータを生成するために指定)

精度の調整に関してはざっくり調べた限り、n_cluster_per_class、flip_yを増やすと分類が難しくなり、class_sepを増やすと簡単になります。また、hyper_cubeをFalseにしたり、informative featureを多めに作って一部の特徴量を削除するといった方法でも分類が難しくなります。分類が簡単すぎてモデル間で差が見えない場合、この辺を調整すると良いと思います。

make_classificationのデータ生成アルゴリズムについて

以下ではmake_classificationで利用されているデータ生成アルゴリズムについて説明しますが、別に知らなくても利用可能なので、興味ある人だけお読みください。

make_classificationで利用されているデータ生成アルゴリズムは下記文献においてMADELONデータの生成に利用されているアルゴリズムのようです。
http://clopinet.com/isabelle/Projects/NIPS2003/Slides/NIPS2003-Datasets.pdf


具体的には下記のような手順でデータを生成しているようです。

  • N(0,1)に従うInformative featureを生成し、クラスタとクラスのラベルに対応づける。その際、各クラスはそれぞれn_cluster_per_class個のラベルに対応するクラスタを持つように構成する。
  • Informative featureからなる行列に-1から1の一様分布からなる要素を持つ行列Aをかけることで、Informative feature間に相関を入れる。
  • n_informative次元で各頂点が±class_sepに位置する超立法体を考え、その頂点に各クラスタを配置する。
  • Informative featureからなる行列に-1から1の一様分布からなる要素を持つ行列Bをかけることで、Redundant featureを生成する。
  • Informative featureとRedundant featureからランダムにRepeated featureを選択する。
  • N(0,1)分布に従うIndependent featureを生成する。
  • N(0,0.1)分布に従うランダムノイズを各特徴量に加算する。(make_classificationでは行ってなさそう)
  • flip_yの比率に応じてランダムにラベルを入れ替える。


簡単に言うと超立方体の頂点にガウス分布に従うデータを作って、それに関係する特徴量や色々なノイズを加えてデータを構成しているイメージです。例えば、n_informativeが2、n_cluster_per_classが2の場合は正方形の頂点に各ラベルに対応したガウス分布データを発生させて色々変換するといった感じです。

データを用いる上での注意点として、構成方法から容易に想像できるように頂点近傍に同種のラベルが集まりやすいアルゴリズムであるため、KNNなどの近傍型アルゴリズムにやや有利なデータが生成されやすいです。(上記のコードでKNNのパフォーマンスが良いのも恐らくこのためと想定されます。)
なので、生成したデータを用いてモデルのパフォーマンスを見るときは、上記のようなデータ構造を持つことを意識して評価した方が良いと思います。

RのtransitionPlot関数を用いた遷移図の作成

こちらの記事で紹介されているRのtransitionPlotを用いた図がとても綺麗でしたので、試しにこれを使って遷移図を作図してみました。

パッケージインストール

transitionPlotはGmiscパッケージに含まれていますがcranサーバから取得できないので、以下のようにリポジトリを指定してソースからインストールする必要があります。

reps = c("http://ftp.sunet.se/pub/lang/CRAN","http://cran.gforge.se")
install.packages("Gmisc", repos=reps, dependencies=TRUE, type="source")
library(Gmisc)
続きを読む

はてなブックマーク記事のレコメンドシステムを作成 PythonによるはてなAPIの活用とRによるモデルベースレコメンド

私は情報収集にはてなブックマークを多用しており、暇な時は結構な割合ではてなブックマークで記事を探してます。しかし、はてなブックマークは最新の記事を探すのは便利ですが、過去の記事を探すにはいまいち使えません。個人的には多少過去の記事でも自分が興味を持っている分野に関しては、レコメンドして欲しいと感じてます。

ありがたいことにはてなはAPIを公開しており、はてなブックマークの情報を比較的簡単に取得できます。そこでこのAPIを利用して自分に合った記事を見つけるようなレコメンド機能をRとPythonで作成してみたいと思います。

続きを読む

[R][データ分析]階層ベイズモデルのサンプルコード bayesmパッケージを利用

Rの階層ベイズモデルのサンプルコードが全然見当たらなかったので、自分で書くことにします。詳細を説明しだすとかなり面倒な領域なので、取り合えず使えるというレベルを目指します。

利用するパッケージは「bayesm」です。
階層ベイズに限らずベイズ推定用MCMCの実行はWinBUGSが一般的だと思いますが、Rのみで利用可能かつ事前分布に関する知識なしで利用可能なのが魅力的なので。

続きを読む

[R]SQLクエリでRのデータを加工・集計できるパッケージ「sqldf」

私が分析を行う際、データ加工や集計作業は基本的にSQLで行い、分析やモデル作成はRで行うことが多いです。
しかし、DBが使えないような場合やちょっとした集計などRでデータを加工・集計したい場合があります。

RでSQLで行うようなデータ加工・集計を行うには、基本的にsubsetやorder、merge、aggregateといった関数を利用します。
SQLとRの関係については、以下のページがわかりやすいです。
 http://d.hatena.ne.jp/a_bicky/20110529/1306667230
ただ、こういった関数はよく使い方を忘れてしまい、Webなりヘルプなりを使うたびに調べるなんてことが起こります。
正直、かなり面倒なわけです。

そこでSQLクエリを使って直接Rのデータを加工できたら良いなーと思うわけですが、
「sqldf」パッケージを用いれば直接SQLを利用してデータを扱うことが可能です。

基本的な使い方は、以下の通りです。

sqldf("SQLクエリ")

SQLクエリ」部分にクエリを書くわけですが、データテーブルとしてR内で生成したデータフレームを用いることができます。

以下、sqldfを利用したデータ加工のサンプルです。

install.packages("sqldf")
library(sqldf)

head(CO2)
summary(CO2)

# 10件抽出
sqldf("SELECT * FROM CO2 limit 10")

# WHEREとORDER BY
sqldf("SELECT Plant, conc, uptake FROM CO2 WHERE Plant IN ('Qn1','Mn1') ORDER BY uptake")

# GROUP BY
sqldf("SELECT Type, COUNT(*) AS CNT FROM CO2 GROUP BY Type")
sqldf("SELECT Plant, Type, COUNT(*) AS CNT, SUM(uptake) AS sum_uptake FROM CO2 GROUP BY Plant, Type")

# DISTINCT
sqldf("SELECT DISTINCT Plant FROM CO2")

set.seed(123)
dat1 <- CO2[sample(1:nrow(CO2),10),]
dat2 <- CO2[sample(1:nrow(CO2),10),]
dat1
dat2

# INNER JOIN
sqldf("SELECT * FROM dat1 AS A INNER JOIN dat2 AS B ON A.Plant = B.Plant")
# LEFT JOIN
sqldf("SELECT * FROM dat1 AS A LEFT JOIN dat2 AS B ON A.Plant = B.Plant")


データをどのように処理しているかというと内部でSQLiteを動かしているらしいので、DB依存の関数には注意しましょう。
もしくはオプションを用いて利用するDBをSQLiteから外部DBに変更することで対応可能です。
その辺の話やsqldfについての詳細はhelpを見るか、以下のページの解説がわかりやすいと思います。
 http://ameblo.jp/wdkz/entry-11115437591.html

注意点としては、以下のようなものがあります。
カラム名にカンマが含まれている場合エラーが起こる。
カラム名に日本語が含まれているとエラーが起こる場合がある。
・RPostgreSQLなどのDB関連パッケージがロードされている場合、内部DBとしてSQLiteよりそちらが優先される。
 (外部DBの利用を意図していない場合、ライブラリをunloadすれば良いです。コマンドはdetach("package:RPostgreSQL", unload=TRUE))


SQLでの集計に慣れている人は、このパッケージで作業が大分楽になるのではないでしょうか。

[Python][データ分析]Pythonによるデータ分析環境

統計分析用ツールとして便利なRですが統計処理に特化しており、プログラム言語としては若干使い難い点もあります。
例えば、Web解析やXMLパースなどはRよりPythonのライブラリの方が機能的に充実していますし、モデル等をPythonで開発するツールに組み込む場合も開発言語と同じ言語で実装した方が効率が良いでしょう。
また、速度的にもRのforはかなり遅くシンンプルなシミュレーションでも結構差が出たりします。

ということでPythonでもデータ分析したい場面が存在しますが、Pythonのリストはデータ分析用のデータ構造ではないため、データハンドリングや集計が非常にやり難いです。
例えば、行・列に対する演算やカテゴリでのグループ化集計といった作業は、非常に面倒な作業になります。
PythonからRを利用するRPy2ライブラリもありますが、PythonとR間でのデータのやり取りが必要になり結構面倒です。

そこで使えるのが、pandasライブラリです。
 http://pandas.pydata.org/

これはPython上でRのようにデータフレームを用いて処理が可能で、GROUP BY集計やapplyによる行・列を対象にした処理、ピボットテーブルまで使うことができるという優れ物です。
HPの記載にもありますが、Rを参考に作ったライブラリらしくデータハンドリング方法がRと非常に似ています。そのためRを使っている人ならすぐ使えるようになると思います。
pandasの使い方について今度まとめる予定です。

これに数値計算ライブラリであるscipyとグラフ表示用のmatplotlibを入れておけば基本的にやりたいことはできるのではないでしょうか。

scipy
 http://www.scipy.org/
matplotlib
 http://matplotlib.org/

あと、Pythonによる機械学習ライブラリとしてはscikit-learnが良いらしいです。
 http://scikit-learn.org/stable/
 http://sucrose.hatenablog.com/entry/2013/05/25/133021
このライブラリはまだ使ったことがないので、今後調査していきたいですね。

取り合えずこの辺が基本でしょうか。
正直分析はRがメインなのでほとんどPythonを用いないのですが、今後の状況次第ではPythonでの分析も選択肢として悪くないかもしれないですね。

[R]SQLのCASE WHEN 表現をRで書く

検定やモデル作成、ヒストグラムを書くといった場合に連続データをカテゴリデータに変換したい場合があります。
このときSQLでデータを変換する場合、例えば以下のようなクエリを書きます。

SELECT
     CASE WHEN revenue < 1000 THEN '1000未満'
          WHEN revenue < 5000 THEN '1000以上5000未満'
          WHEN revenue < 10000 THEN '5000以上10000未満'
          ELSE '10000以上'
     END AS revenue_category
     ,COUNT(*) AS CNT
FROM test_data
GROUP BY
     CASE WHEN revenue < 1000 THEN '1000未満'
          WHEN revenue < 5000 THEN '1000以上5000未満'
          WHEN revenue < 10000 THEN '5000以上10000未満'
          ELSE '10000以上'
     END

このような集計をRで実行しようとした場合、CASE WHENに対応するものとしてifelse関数を用います。
例えば上記と同様の処理をRで実行しようとした場合、以下のようになります。

revenue_category <- ifelse(revenue < 1000, "1000未満",
					ifelse(revenue < 5000, "1000以上5000未満",
					ifelse(revenue < 10000, "5000以上10000未満","10000以上")))
table(revenue_category)

単純にデータの傾向を見たい場合は代表値である平均や分散を見れば良いですが、
外れ値がある場合や分布が正規分布ではない場合、カテゴリ化して集計した方が安全です。
特にマーケティングにおいては、一部の人の売上が全体の80%を占めるという「ジップの法則」があるため、
平均や分散で評価が難しいことが多いです。
例えば、2つのキャンペーンを比べる際、外れ値が2,3個あるだけで全体傾向はほぼ同じであるのに、
平均値で評価すると効果が大きく異なっている、と判断しかねません。
(上記のような場合は中央値を見れば改善されますが、正規性を仮定できないデータではやはり分布を見た方が安全です。)

カテゴリ集計によって得られる結果は、基本的に「XXという行動をする人は(YYする人に比べて)ZZ円以上の割合が多い」
といったものになります。
必要があれば別軸での分解を行いながらこれらの結果を積み上げて、仮説を立案・検証していくわけです。

Rでカテゴリ分類ができれば、プログラムを組んでクロス集計やカテゴリ別検定といった作業を自動化でき、
作業効率改善に繋がるわけです。
全カラムのクロス集計とかカイ2乗検定とかカラムが少なければエクセルのピボットテーブルなどで可能ですが、
流石に1000カラムを超えるようなデータでは現実的に無理ですからね。