Rust で重回帰
簡単なアルゴリズムを実装しつつ Rust を学びたい。環境はこちら。
- Rust(Nightly) rustc 1.6.0-nightly (d5fde83ae 2015-11-12)
- rust-csv 0.14.3 : CSV の読み込みに利用
- nalgebra 0.3.1 : 行列/ベクトルの処理に利用
Python や R と比較すると、型のために Vec
や slice
が透過的に扱えず戸惑う。また、
Generics を使った場合に所有権まわりのエラーが出て解決できなかったため、浮動小数点 f64
のみの定義とした。もう少し理解できたら書き直したい。
extern crate csv; extern crate nalgebra; use std::vec::Vec; use nalgebra::{DVec, DMat, Inv, Mean, ColSlice, Iterable}; fn main() { // R の "trees" データを使用 let data = "Girth,Height,Volume 8.3,70,10.3 8.6,65,10.3 8.8,63,10.2 10.5,72,16.4 10.7,81,18.8 10.8,83,19.7 11.0,66,15.6 11.0,75,18.2 11.1,80,22.6 11.2,75,19.9 11.3,79,24.2 11.4,76,21.0 11.4,76,21.4 11.7,69,21.3 12.0,75,19.1 12.9,74,22.2 12.9,85,33.8 13.3,86,27.4 13.7,71,25.7 13.8,64,24.9 14.0,78,34.5 14.2,80,31.7 14.5,74,36.3 16.0,72,38.3 16.3,77,42.6 17.3,81,55.4 17.5,82,55.7 17.9,80,58.3 18.0,80,51.5 18.0,80,51.0 20.6,87,77.0"; let mut ncols = 0; // http://burntsushi.net/rustdoc/csv/ let mut reader = csv::Reader::from_string(data).has_headers(true); let mut x: Vec<f64> = vec![]; for row in reader.decode() { let vals: Vec<f64> = row.unwrap(); // 1 列目を目的変数, 残りは説明変数とする for i in 0 .. vals.len() { x.push(vals[i]); } // 説明変数の数 (各行は必ず同数のデータを含むと仮定) ncols = vals.len(); } // レコード数 let nrows: usize = x.len() / ncols; // ベクトルから nrows x ncols の DMat (動的サイズの行列) を作成 // http://nalgebra.org/doc/nalgebra/struct.DMat.html let dx = DMat::from_row_vec(nrows, ncols, &x); // 重回帰 let coefs = lm_fit(&dx); println!("結果 {:?}", &coefs); } fn lm_fit(data: &DMat<f64>) -> Vec<f64> { // 重回帰 // 列ごとの平均値を計算 let means = data.mean(); let nrows = data.nrows(); let nfeatures = data.ncols() - 1; // 偏差平方和積和行列 // DMat::from_fn で、行番号 i, 列番号 j を引数とする関数から各要素の値を生成できる let smx = DMat::from_fn(nfeatures, nfeatures, |i, j| sum_square(&data.col_slice(i + 1, 0, nrows), &data.col_slice(j + 1, 0, nrows), means[i+1], means[j+1])); // 偏差積和行列 // DMat と Vec では演算ができないため、こちらも一列の DMat として生成 let smy = DMat::from_fn(nfeatures, 1, |i, j| sum_square(&data.col_slice(i + 1, 0, nrows), &data.col_slice(0, 0, nrows), means[i], means[0])); // 偏回帰係数を計算し、Vec に変換 let mut res = (smx.inv().unwrap() * smy).to_vec(); // 切片を計算し、0 番目の要素として挿入 let intercept = (1..means.len()).fold(means[0], |m, i| m - res[i-1]*means[i]); res.insert(0, intercept); return res } fn sum_square(vec1: &DVec<f64>, vec2: &DVec<f64>, m1: f64, m2: f64) -> f64 { // 平方和 let mut val: f64 = 0.; for (v1, v2) in vec1.iter().zip(vec2.iter()) { val += (v1 - m1) * (v2 - m2); } return val; }
実行結果。左から切片と各説明変数の偏回帰係数。
# 結果 [10.81637050444831, -0.04548349100904309, 0.19517974893553702]
R の結果。
lm(Girth ~ Height + Volume, data = trees) # # Call: # lm(formula = Girth ~ Height + Volume, data = trees) # # Coefficients: # (Intercept) Height Volume # 10.81637 -0.04548 0.19518
2016/12/26追記 上では逆行列を使ったが、利用するパッケージに線形方程式を解く関数 / メソッドが用意されていればそちらを使った方が数値的に安定 / 高速。
岩波データサイエンス Vol.1
ご恵贈いただきました。 ありがとうございます! あわせてタスクもいただきました (下部)。
書籍のコンテンツ
各章ごとの内容は id:sfchaos さんが詳しく紹介されています。
まだ すべて読めていないのですが、以下 3 点がよいポイントだと思います。
- 理論 と サンプルプログラム 両方の記載がある
BUGS
,Stan
,PyMC3
と主要なパッケージが網羅されている- サンプルは単純な回帰だけでなく 時系列 / 空間ベイズを含む
補足 書籍には コラム "Pythonとは" という データ分析視点での Python 紹介があるのですが、中身は結構な pandas
推しでした。著者の方、いったい何者なんだ...。
Stan 入門
依頼により、著者の松浦さんが作成した RStan
サンプルの PyStan
版を作成させていただきました。
以下リポジトリの "Example..." フォルダ中に含まれる Jupyter Notebook ( .ipynb
) を開くと、 GitHub 上でプログラムと出力を確認することができます。
RStan
と PyStan
の API 差異については公式ドキュメントに記載がありますがアップデートがされていません。下表にサンプル中で利用したものを簡単にまとめます。
RStan |
PyStan |
概要 |
---|---|---|
stan(...) |
pystan.stan(...) |
モデルのコンパイルとサンプリングの実行 |
stan_model(...) |
pystan.StanModel(...) |
モデルのコンパイル |
sampling(...) |
StanModel.sampling(...) |
サンプリングの実行 |
extract(...) |
StanFit4model.extract() |
サンプリング結果を取得 |
traceplot(...) |
StanFit4model.plot(...) |
サンプリングの経過をプロット |
書籍や Rstan
のサンプルと合わせてご参照ください。
- 作者: 岩波データサイエンス刊行委員会
- 出版社/メーカー: 岩波書店
- 発売日: 2015/10/08
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (4件) を見る
pandas 0.17.0 の主要な変更点
先日 10/9 に pandas
0.17.0 がリリースされた。直近のバージョンアップの中では かなり機能追加が多いリリースとなった。
重要な変更は リリースノート にハイライトとして列挙しているのだが、これらはある程度 pandas
を使いこなしている方向けの記載となっている。
そのため、ここでは よりライトなユーザ向けに重要と思われる変更を書く。特に、ユーザ側のプログラムに影響がある以下の3点について記載する。
- ソート API の統合 (
sort_values
/sort_index
) - 重複削除 API の改善 (
drop_duplicates
/duplicated
) .plot
アクセサの追加
準備
import numpy as np import pandas as pd np.__version__ # '1.10.1' pd.__version__ # u'0.17.0'
1. ソート API の統合 ( sort_values
/ sort_index
)
これまで一貫性がなかったソート API が以下のように統合された。以前のメソッド / 引数は deprecate されており、将来のバージョン ( 通例では 0.19.0 以降 ) で削除される。
値によるソート ( sort_values
として統合):
0.16.2以前 | 0.17.0以降 |
---|---|
Series.order() |
Series.sort_values() |
Series.sort() |
Series.sort_values(inplace=True) |
DataFrame.sort(columns=...) |
DataFrame.sort_values(by=...) |
ラベル ( Index
) によるソート ( sort_index
として統合):
0.16.2以前 | 0.17.0以降 |
---|---|
Series.sort_index() |
Series.sort_index() |
Series.sortlevel(level=...) |
Series.sort_index(level=...) |
DataFrame.sort_index() |
DataFrame.sort_index() |
DataFrame.sortlevel(level=...) |
DataFrame.sort_index(level=...) |
DataFrame.sort() |
DataFrame.sort_index() |
それぞれ、引数として以下のオプションが利用できる。
inplace
: 破壊的にソートするかどうかを指定する。既定はFalse
(非破壊)axis
: ソートの方向(軸)を指定するascending
: ソート順序 (昇順/降順) を指定する。既定はTrue
(昇順)
s = pd.Series(list('afcedb'), index=[1, 3, 2, 4, 0, 5]) s # 1 a # 3 f # 2 c # 4 e # 0 d # 5 b # dtype: object # 値によるソート s.sort_values() # 1 a # 5 b # 2 c # 0 d # 4 e # 3 f # dtype: object # ラベル (index) によるソート s.sort_index() # 0 d # 1 a # 2 c # 3 f # 4 e # 5 b # dtype: object
補足 もともとの仕様は NumPy
(というか Python) の .sort
メソッドの挙動 ( 破壊的なソート) にあわせたものだった。今回の変更により NumPy
とは差異が発生する。
a = np.array([3, 2, 1]) a # array([3, 2, 1]) a.sort() a # array([1, 2, 3])
2. 重複削除 API の改善 ( drop_duplicates
/ duplicated
)
重複した値の削除を行う drop_duplicates
で、重複値の全削除ができるようになった。
0.16.2以前 | 0.17.0以降 |
---|---|
.drop_duplicates() |
.drop_duplicates() |
.drop_duplicates(take_last=True) |
.drop_duplicates(keep='last') |
- | .drop_duplicates(keep=False) |
削除方法の指定は これまでの take_last
から keep
オプション に変更された。
s = pd.Series(list('abebdcd')) s # 0 a # 1 b # 2 e # 3 b # 4 d # 5 c # 6 d # dtype: object # 重複がある場合、最初に出現した値を残す s.drop_duplicates() # 0 a # 1 b # 2 e # 4 d # 5 c # dtype: object # 重複がある場合、最後に出現した値を残す s.drop_duplicates(keep='last') # 0 a # 2 e # 3 b # 5 c # 6 d # dtype: object # 重複した値を残さない s.drop_duplicates(keep=False) # 0 a # 2 e # 5 c # dtype: object
また、値が重複しているかどうかを調べる .duplicated
でも keep
オプションが利用できる。
# 重複がある場合、最初に出現した値以外を True s.duplicated() # 0 False # 1 False # 2 False # 3 True # 4 False # 5 False # 6 True # dtype: bool # 重複がある場合、最後に出現した値以外を True s.duplicated(keep='last') # 0 False # 1 True # 2 False # 3 False # 4 True # 5 False # 6 False # dtype: bool # 重複 全てを True s.duplicated(keep=False) # 0 False # 1 True # 2 False # 3 True # 4 True # 5 False # 6 True # dtype: bool
Index
の値から重複を全削除する場合は Index.duplicated
を以下のように利用すればよい。
s = pd.Series(np.arange(6), index=list('abcbda')) s # a 0 # b 1 # c 2 # b 3 # d 4 # a 5 # dtype: int64 s[~s.index.duplicated(keep=False)] # c 2 # d 4 # dtype: int64
3. .plot
アクセサの追加
DataFrame
や Series
をプロットする際のグラフの種類が、 .plot.bar()
や .plot.hist()
のように .plot
をアクセサとして指定できるようになった。
これまでのメソッド呼び出し ( .plot(kind='bar')
) は引き続き利用できる。
df = pd.DataFrame(np.random.randn(5, 10)) # df.plot(kind='box') と同じ df.plot.box()
その他
また、上記ほどではないが重要な内容を列挙する。気になるものがあれば詳細はリンク先で。
- GIL 解放 (
Dask
利用時のパフォーマンス向上) (詳細) pandas.io.data
を deprecate し、別パッケージpandas-datareader
として分離 (詳細)pd.to_datetime
のエラー処理の変更、日時パース時の挙動統一 (詳細)Index
同士の比較演算時の一部挙動の変更 (詳細)- ターミナル上での日本語データ表示時の位置補正オプションの追加 (詳細)
まとめ
pandas
0.17.0 での特に重要な変更点 3 点を記載した。
- ソート API の統合 (
sort_values
/sort_index
) - 重複削除 API の改善 (
drop_duplicates
/duplicated
) .plot
アクセサの追加
これら以外の変更点については リリースノート を一読お願いします。
10/18 編集 コメントご指摘により誤記を修正
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (11件) を見る
PyConJP 2015: pandas/Daskについてお話させていただきました
10日、11日と PyCon JP に参加させていただきました。ご参加いただいた皆様、スタッフの皆様ありがとうございました。資料はこちらになります。
pandas internals
パフォーマンス向上のための pandas 内部実装の説明といくつかの TIPS について。そのうち翻訳するかもしれません。
Dask: 軽量並列分散フレームワーク (LT)
元ネタ
以下のエントリをベースに、それぞれ新しい内容を追加しています。
Python XGBoost + pandas 連携の改善
一部 こちらの続き。その後 いくつかプルリクを送り、XGBoost
と pandas
を連携させて使えるようになってきたため、その内容を書きたい。
できるようになったことは 以下 3 点。
DMatrix
でのラベルと型の指定pd.DataFrame
からのDMatrix
の作成xgb.cv
の結果をpd.DataFrame
として取得
補足 XGBoost では PyPI の更新をスクリプトで不定期にやっているようで、同一バージョンに見えても枝番が振られていたりして見分けにくい。記載は本日時点のこのコミットの情報。
%matplotlib inline import numpy as np import xgboost as xgb from sklearn import datasets import matplotlib.pyplot as plt plt.style.use('ggplot') xgb.__version__ # '0.4'
1. DMatrix
でのラベルと型の指定
これは pandas
関係ない。DMatrix
に以下 2 つのプロパティが追加され、任意のラベル/型が指定できるようにした (指定しない場合はこれまでと同じ挙動)。
feature_names
: 変数のラベルを指定。ラベルは英数字のみ。feature_types
: 変数の型を指定。指定できるのは、q
: 量的変数 (既定)i
: ダミー変数int
: 整数型float
: 浮動小数点型。
まずは feature_names
のみを指定。
iris = datasets.load_iris() features = ['SepalLength', 'SepalWidth', 'PetalLength', 'PetalWidth'] dm = xgb.DMatrix(iris.data, label=iris.target, feature_names=features) dm.feature_names # ['SepalLength', 'SepalWidth', 'PetalLength', 'PetalWidth']
DMatrix
で指定した feature_names
は xgb.Booster
へも引き継がれ、.get_dump
時 やプロットに利用できる。
params={'objective': 'multi:softprob', 'eval_metric': 'mlogloss', 'eta': 0.3, 'num_class': 3} np.random.seed(1) bst = xgb.train(params, dm, num_boost_round=18) xgb.plot_importance(bst)
次に feature_types
。これは内部的な学習には全く関係がなく、.get_dump
や xgb.plot_tree
をする際の ラベルのフォーマットを決めるだけ。例えば以下のようなデータがあったとする。
data = np.array([[1, 0.1, 0], [2, 0.2, 1], [2, 0.2, 0]]) data # array([[ 1. , 0.1, 0. ], # [ 2. , 0.2, 1. ], # [ 2. , 0.2, 0. ]])
このとき、1 列目を int
, 2 列目を float
, 3 列目をダミー変数 i
として扱いたければ以下のように指定すると、 plot_tree
したときの見た目が若干変わる。全変数 既定 (q
: 量的変数) でもとくに不自由ないので、強いこだわりのある方以外は指定する必要はない。
dm = xgb.DMatrix(data, feature_names=['A', 'B', 'C'], feature_types=['int', 'float', 'i']) dm.feature_types # ['int', 'float', 'i']
補足 これまで、変数のラベル/型の指定は "featmap.txt" のような設定ファイルから別途指定する必要があった。この指定が Python で直接できるようになった。
2 pd.DataFrame
からの DMatrix
の作成
pandas.DataFrame
から DMatrix
が作成できるようにした。このとき、feature_names
と feature_types
は DataFrame
の定義から適当に設定される。
補足 DMatrix
の label
には pd.Series
が渡せる (前からできた)。
import pandas as pd # 表示行数を指定 pd.set_option('display.max_rows', 8) # 一行に表示する文字数を指定 pd.set_option('display.width', 120) # pd.DataFrame を作成 train = pd.DataFrame(iris.data, columns=['SepalLength', 'SepalWidth', 'PetalLength', 'PetalWidth']) train # SepalLength SepalWidth PetalLength PetalWidth # 0 5.1 3.5 1.4 0.2 # 1 4.9 3.0 1.4 0.2 # 2 4.7 3.2 1.3 0.2 # 3 4.6 3.1 1.5 0.2 # .. ... ... ... ... # 146 6.3 2.5 5.0 1.9 # 147 6.5 3.0 5.2 2.0 # 148 6.2 3.4 5.4 2.3 # 149 5.9 3.0 5.1 1.8 # # [150 rows x 4 columns] # DMatrix を作成 dm = xgb.DMatrix(train, label=pd.Series(iris.target)) dm.feature_names # ['SepalLength', 'SepalWidth', 'PetalLength', 'PetalWidth'] dm.feature_types # ['q', 'q', 'q', 'q']
DMatrix
に変換できるのは、列が int64
, float64
, bool
の DataFrame
のみ。他の型を自動でダミー変数に変換したりはしない。
df = pd.DataFrame([[1, 0.1, False], [2, 0.2, True], [2, 0.2, False]], columns=['A', 'B', 'C']) df # A B C # 0 1 0.1 False # 1 2 0.2 True # 2 2 0.2 False df.dtypes # A int64 # B float64 # C bool # dtype: object dm = xgb.DMatrix(df) dm.feature_names # ['A', 'B', 'C'] dm.feature_types # ['int', 'q', 'i'] # object 型を含むためエラー df = pd.DataFrame([[1, 0.1, 'x'], [2, 0.2, 'y'], [2, 0.2, 'z']], columns=['A', 'B', 'C']) df # A B C # 0 1 0.1 x # 1 2 0.2 y # 2 2 0.2 z df.dtypes # A int64 # B float64 # C object # dtype: object xgb.DMatrix(df) # ValueError: DataFrame.dtypes must be int, float or bool
3. xgb.cv
の結果を pd.DataFrame
として取得
クロスバリデーションを行う xgb.cv
はこれまで実行結果を文字列で返していたため、その結果をパースもしくは目視で確認する必要があった。この返り値を pd.DataFrame
もしくは np.ndarray
とし、プログラムで処理がしやすいようにした。
train = pd.DataFrame(iris.data, columns=['SepalLength', 'SepalWidth', 'PetalLength', 'PetalWidth']) dm = xgb.DMatrix(train, label=pd.Series(iris.target)) cv = xgb.cv(params, dm, num_boost_round=50, nfold=10) cv # test-mlogloss-mean test-mlogloss-std train-mlogloss-mean train-mlogloss-std # 0 0.753459 0.027033 0.737631 0.003818 # 1 0.552303 0.048738 0.526929 0.005102 # 2 0.423481 0.066469 0.390115 0.005873 # 3 0.339942 0.082163 0.295637 0.006148 # .. ... ... ... ... # 46 0.208001 0.259161 0.018070 0.001759 # 47 0.208355 0.261166 0.017898 0.001724 # 48 0.208468 0.261520 0.017755 0.001703 # 49 0.208566 0.260967 0.017617 0.001686 # # [50 rows x 4 columns] cv.sort(columns='test-mlogloss-mean') # test-mlogloss-mean test-mlogloss-std train-mlogloss-mean train-mlogloss-std # 11 0.175506 0.177823 0.058496 0.004715 # 13 0.175514 0.194641 0.045222 0.004117 # 14 0.176150 0.203265 0.040493 0.004050 # 12 0.176277 0.186902 0.051074 0.004423 # .. ... ... ... ... # 3 0.339942 0.082163 0.295637 0.006148 # 2 0.423481 0.066469 0.390115 0.005873 # 1 0.552303 0.048738 0.526929 0.005102 # 0 0.753459 0.027033 0.737631 0.003818 # # [50 rows x 4 columns]
pandas
がインストールされていない、もしくは as_pandas=False
を指定した場合、返り値は np.ndarray
となる。
cv = xgb.cv(params, dm, num_boost_round=50, nfold=10, as_pandas=False) cv # array([[ 0.7534586 , 0.02703308, 0.7376308 , 0.00381775], # ... # [ 0.2085664 , 0.26096721, 0.0176169 , 0.00168606]]) np.argmin(cv, axis=0) # array([11, 0, 49, 49])
まとめ
XGBoost と pandas
をより便利に使うため、以下 3 点の修正を行った。
DMatrix
でのラベルと型の指定pd.DataFrame
からのDMatrix
の作成xgb.cv
の結果をpd.DataFrame
として取得
これ使って Kaggle で大勝利したい...と思ったのですが、他の日本勢に使っていただいたほうが勝率高いはずなのでシェアします。
Python Dask で 並列 DataFrame 処理
はじめに
先日のエントリで少し記載した Dask
について、その使い方を書く。Dask
を使うと、NumPy
や pandas
の API を利用して並列計算/分散処理を行うことができる。また、Dask
は Out-Of-Core (データ量が多くメモリに乗らない場合) の処理も考慮した実装になっている。
上にも書いたが、Dask
は NumPy
や pandas
を置き換えるものではない。数値計算のためのバックエンドとして NumPy
や pandas
を利用するため、むしろこれらのパッケージが必須である。
Dask
は NumPy
や pandas
の API を完全にはサポートしていないため、並列 / Out-Of-Core 処理が必要な場面では Dask
を、他では NumPy
/ pandas
を使うのがよいと思う。pandas
とDask
のデータはそれぞれ簡単に相互変換できる。
補足 とはいえ都度の変換は手間なので、pandas
の処理実行時に Dask
を利用するオプションをつける という検討はされている。
インストール
pip
もしくは conda
で。
pip install dask
準備
まずは必要なパッケージを import
する。
import numpy as np import pandas as pd import dask.dataframe as dd np.__version__ # '1.9.3' pd.__version__ # '0.16.2' # バージョン表示のためにインポート。 import dask dask.__version__ # '0.7.1'
pandas
から Dask
への変換
サンプルデータは すでにメモリ上にある pd.DataFrame
とする。
df = pd.DataFrame({'X': np.arange(10), 'Y': np.arange(10, 20), 'Z': np.arange(20, 30)}, index=list('abcdefghij')) df # X Y Z # a 0 10 20 # b 1 11 21 # c 2 12 22 # d 3 13 23 # e 4 14 24 # f 5 15 25 # g 6 16 26 # h 7 17 27 # i 8 18 28 # j 9 19 29
pandas
のデータ構造から Dask
に変換するには dd.from_pandas
。2つめの引数で データをいくつのパーティションに分割するかを指定している。結果は dask.dataframe.DataFrame
( dd.DataFrame
) となる。
divisions
はデータがどこで分割されたかを示す。表示から、1つ目のパーティションには .index
が "a" 〜 "e" までのデータが、2つ目のには "f" 〜 "j" までのデータが含まれていることがわかる。
重要 dd.DataFrame
の処理全般について、行の順序は保証されない。各パーティションには divisions
で示される .index
を持つ行が含まれるが、パーティション内が常にソートされているとは限らない。
ddf = dd.from_pandas(df, 2) ddf # dd.DataFrame<from_pandas-b08addf72f0693a9fa1bb6c21d91f3d4, divisions=('a', 'f', 'j')> # DataFrame の列名 ddf.columns # ('X', 'Y', 'Z') # DataFrame の index ddf.index # dd.Index<from_pandas-b08addf72f0693a9fa1bb6c21d91f3d4-index, divisions=('a', 'f', 'j')> # DataFrame の divisions (パーティションの分割箇所) ddf.divisions # ('a', 'f', 'j') # DataFrame のパーティション数 ddf.npartitions # 2
dd.DataFrame
からデータを取得する (計算処理を実行する) には .compute()
。結果、元の pd.DataFrame
が得られる。
ddf.compute() # X Y Z # a 0 10 20 # b 1 11 21 # c 2 12 22 # d 3 13 23 # e 4 14 24 # f 5 15 25 # g 6 16 26 # h 7 17 27 # i 8 18 28 # j 9 19 29
内部処理
ここから、dd.DataFrame
でどういった処理ができるのか、内部動作とあわせて記載する。といっても難しいことは全くやっていない。
まず、データ全体について 1 加算する処理を考える。これは 各パーティションごとに 1 加算して連結するのと同じ。
ddf + 1 # dd.DataFrame<elemwise-5b9ae0407254158903343113fac41421, divisions=('a', 'f', 'j')> (ddf + 1).compute() # 略
次に、列ごとの合計をとる処理。これは、各パーティションごとに列の合計をとって連結し、もう一度 合計をとる処理と同じ。
列の合計をとるため、結果は Series
になる。ddf.sum()
の時点では .compute()
が呼ばれていないため実際の計算は行われていないが、Dask
が結果の型 や divisions
を推測し正しく表示している。
ddf.sum() # dd.Series<dataframe-sum--second-7ba12c9d58c17f61406b84b6c30d7a26, divisions=(None, None)> ddf.sum().compute() # X 45 # Y 145 # Z 245 # dtype: int64
Dask
ではこのような形で、計算処理をパーティションごとに並列 / Out-Of-Core 実行できる形に読み替えている。これらの処理は内部的には Computational Graph ( Dask Graph ) として表現され、.compute()
によって実行される。
各処理の Dask Graph は、.visualize()
メソッドを利用して確認できる。Graph 上で縦につながっていない処理同士は並列で実行できる。
ddf.sum().visualize()
各列の平均をとる場合、内部的には各列の .sum()
と 各列の .count()
をそれぞれ計算して除算。
ddf.mean().compute() # X 4.5 # Y 14.5 # Z 24.5 # dtype: float64 ddf.mean().visualize()
DataFrame
同士の演算や、演算をチェインすることもできる。互いのパーティションが異なる場合はそれらが一致するよう調整が行われる。
((ddf - (ddf * 2)) == - ddf).visualize()
また、累積関数 ( cumxxx
) や ウィンドウ関数 ( rolling_xxx
) なども利用できる。
ddf.cumsum().compute() # X Y Z # a 0 10 20 # b 1 21 41 # c 3 33 63 # d 6 46 86 # e 10 60 110 # f 15 75 135 # g 21 91 161 # h 28 108 188 # i 36 126 216 # j 45 145 245 ddf.cumsum().visualize()
concat
, join
などの 連結 / 結合もできる。通常の演算と同じく、dd.DataFrame
同士のパーティションは適当に調整される。
df2 = pd.DataFrame({'A': np.arange(5), 'B': np.arange(10, 15)}, index=list('adefg')) df2 # A B # a 0 10 # d 1 11 # e 2 12 # f 3 13 # g 4 14 ddf2 = dd.from_pandas(df2, 2) ddf2 # dd.DataFrame<from_pandas-667963fc37e22688843f02da80df5963, divisions=('a', 'f', 'g')> ddf.join(ddf2).compute() # X Y Z A B # a 0 10 20 0 10 # b 1 11 21 NaN NaN # c 2 12 22 NaN NaN # d 3 13 23 1 11 # e 4 14 24 2 12 # f 5 15 25 3 13 # g 6 16 26 4 14 # h 7 17 27 NaN NaN # i 8 18 28 NaN NaN # j 9 19 29 NaN NaN ddf.join(ddf2).visualize()
サポートされている処理の一覧は以下のAPIドキュメントを。一部利用できない引数が明記されていないが、次バージョンにて改訂。
9/26 追記 処理結果については、行の順序以外は pandas
の処理と一致するはず。例外は quantile
のような percentile をとる処理。これらは Out-Of-Core 処理のための近似アルゴリズムを使っており、正確な値とずれることがある。
実データでの利用例
こちらが良エントリ (英語)。
- Analyzing Reddit Comments with Dask and Castra
- Out-of-Core Dataframes in Python: Dask and OpenStreetMap
まとめ
Dask
を利用して DataFrame
を並列処理する方法を記載した。手順は、
dd.from_pandas
を利用してpd.DataFrame
をdd.DataFrame
へ変換。- 実行したいメソッド / 演算を
dd.DataFrame
に対して適用。 .compute()
で計算を実行し、結果を取得する。計算処理はDask
にて自動的に並列化される。
最後、pandas
0.16.2 時点では並列処理による速度向上は大きくはない。これは Python の GIL (Global Interpreter Lock ) により並列実行できる処理が限定されているため。今月中にリリース予定の pandas
0.17.0 では いくつかの処理で Cython から明示的に GIL 解放するよう実装を変更しており、並列化による速度向上は大きくなる。
Python 次世代の多次元配列パッケージ群
このところ、たびたび NumPy
後継が...とか 並列処理が...という話を聞くので、この秋 注目の多次元配列パッケージをまとめたい。
バックエンド系
NumPy
のように数値計算処理を自前で実装しているパッケージ。
DyND
Blaze
プロジェクトのひとつ。C++ 実装 + Python バインディング。GitHub にいくつか Example があがっているが、複合型やカテゴリカル型、GroupBy 操作がサポートされていて熱い。ラベルデータも NumPy
より簡単に実装できそうだ。
並列分散系
自身では直接 数値計算処理を行わず、バックエンド ( 主に NumPy
)を利用して並列/分散処理を行うパッケージ。1 物理PC/複数コアでの並列計算を主用途とし、NumPy
, pandas
では少し苦しいが PySpark
などを使うほどじゃない...という場合に利用するもの。
Dask
Blaze
プロジェクトのひとつ。Pure-Python 実装。主用途は単一物理PCでの複数コア計算だが、複数PCに処理を分散させることもできる。
NumPy
の API で並列計算を行う dask.array
のほか、toolz
相当の操作を行う dask.bag
、pandas
相当の処理を行う dask.dataframe
など一連のデータ構造が揃っている。
開発者はベースパッケージである NumPy
、toolz
、pandas
のコミッタとの兼任が多い。自分もコミット権限をいただいており、dask.dataframe
への API 追加を行なっている。
DistArray
NumPy
の API で並列計算を行う DistArray
のみをサポート、ほかのデータ構造はない。Pure-Python 実装。見た感じ、並列処理の基本的な考え方は Dask
と同じようだ。Enthought が好きな方はこちらを使えばよいかと。
bolt
NumPy
の API でローカル計算/分散処理を行う。Pure-Python実装。ローカル計算は NumPy
をそのまま使い (並列処理しない)、分散処理は Spark
で行う。海外の方が期待の新星っぽい扱いをしていたので気になったのだが、現時点で実装されている API は多くはない。
ほかにもいくつかプロジェクトがあるが、自分としては Dask
を流行らせたいので、その記事を書きます。
9/24 追記 書きました。
NumPy でつくる俺々データ構造
はじめに
Python での数値計算の基盤をなす NumPy 、直感的なスライスやブロードキャスト、関数のベクトル適用など大変便利だ。
import numpy as np np.__version__ # '1.9.2' np.array([1, 2, 3]) # array([1, 2, 3]) np.array([1, 2, 3])[:2] # array([1, 2]) np.array([1, 2, 3]) + 1 # array([2, 3, 4])
が、用途によっては NumPy 標準ではその機能を実現できない場合がある。例えば、
- 配列とメタデータをひとつのクラスで扱いたい
- 配列への入力や型を制約/検証したい
- 自作クラスを NumPy の Universal Functions (ufunc) に対応させたい
- 新しい型 (
dtype
) を作りたい
こういったとき、NumPy の動作を継承した自作のクラスが作れるとうれしい。方法として、大きく二つがある。
- NumPy
ndarray
を継承する - NumPy Array Interface を実装する
どちらを使うべきかは目的による。サブクラス化しないほうが柔軟な処理ができるが、互換をとるための実装は多くなる。pandas
の主要データ構造では 柔軟性を重視して ndarray
の継承をやめ、現在は 2 の方法を取っている。
ndarray
のサブクラスを作成する方法は以下のドキュメントに記載されている。きちんとやりたい場合はこちらを読むのがよい。
このエントリでは 2. Array Interface を使って最低限の動作をさせる方法を書く。
クラスの定義
サンプルとして、名前をもつ配列クラス NamedArray
を考える。__init__
での初期化時に、名前をあらわす name
と配列の値 values
を渡すことにする。
class NamedArray(object): def __init__(self, name, values): self.name = name self.values = np.array(values) def __repr__(self): return 'NamedArray: {0}: {1}'.format(self.name, self.values) def __len__(self): return len(self.values) n = NamedArray('x', [1, 2, 3]) n # NamedArray: x: [1 2 3] len(n) # 3
このクラスに ndarray
と互換の挙動をさせたい。まず、NamedArray
インスタンスから ndarray
が直接作成できるとうれしい。が、NamedArray
を np.array
に渡すと shape
のない ndarray
(0-dimmentional array) が作成されてしまう。
np.array(n) # array(NamedArray: x: [1 2 3], dtype=object) np.array(n).shape # ()
理由は np.array
が NamedArray
をスカラーとして扱うため。期待の動作をさせるためには、np.array
に NamedArray
クラスが配列であることを教える必要がある。
具体的には NamedArray
に Array Interface ( __array__
) を定義すればよい。以下のように __array__
メソッドを追加したクラス NPNamedArray
を作成すると、期待通り values
の値から ndarray
が作成されていることがわかる。
class NPNamedArray(NamedArray): def __array__(self): return self.values n = NPNamedArray('x', [1, 2, 3]) np.array(n) # array([1, 2, 3]) np.array(n).shape # (3,)
ufunc 処理のオーバーライド
また、Array Interface をもつインスタンスには ufunc が適用できるようになる。引数は ufunc 内で ndarray
に変換されるため、返り値は ndarray
になる。
np.sin(n)
# array([ 0.84147098, 0.90929743, 0.14112001])
このときの内部の動きをみてみる。ufunc 適用後の処理は __array_wrap__
メソッドによって制御することができる。
__array_wrap__
メソッドには out_arr
と context
2 つの引数を定義する必要がある。これらに何が渡されているのか、print
して確認する。
class NPNamedArray(NamedArray): def __array__(self): return self.values def __array_wrap__(self, out_arr, context=None): print('out_arr: ', out_arr) print('context: ', context) return out_arr n = NPNamedArray('x', [1, 2, 3]) np.sin(n) # ('out_arr: ', array([ 0.84147098, 0.90929743, 0.14112001])) # ('context: ', (<ufunc 'sin'>, (NamedArray: x: [1 2 3],), 0)) # array([ 0.84147098, 0.90929743, 0.14112001])
out_arr
には ufunc が適用された結果の ndarray
が、 context
には ufunc を呼び出した際のシグネチャ ( ufunc 自体, ufunc の引数, ufunc のドメイン ) が tuple
で渡ってきている。__array_wrap__
の返り値を NPNamedArray
に書き換えてやれば、期待する動作ができそうだ。
補足 np.isreal
や np.ones_like
など、Universal Functions に記載されていても __array_wrap__
を通らないものもある。
補足 ufunc のドメインは NumPy の場合 (おそらく) 全部 0。
class NPNamedArray(NamedArray): def __array__(self): return self.values def __array_wrap__(self, out_arr, context=None): return NPNamedArray(self.name, out_arr) n = NPNamedArray('x', [1, 2, 3]) np.sin(n) # NamedArray: x: [ 0.84147098 0.90929743 0.14112001]
ここまでの時点で、NPNamedArray
には __add__
メソッドを定義していないため、通常の加算はエラーになる。
n + 1 # TypeError: unsupported operand type(s) for +: 'NPNamedArray' and 'int'
が、ufunc である np.add
や ndarray
を含む演算は利用できるようになる。これは先日のエントリに記載した通り ndarray
側に Array Interface をもつクラスとの演算の定義があるため。
np.add(n, 1) # NamedArray: x: [2 3 4] n + np.array([5, 5, 5]) # NamedArray: x: [6 7 8]
また、__array_wrap__
の context
引数の値によって、特定の ufunc の処理をオーバーライドすることができる。例えば np.sin
を適用したときに任意の値を返すこともできる。
class NPNamedArray(NamedArray): def __array__(self): return self.values def __array_wrap__(self, out_arr, context=None): if context[0] is np.sin: return u'返したい値' return NPNamedArray(self.name, out_arr) n = NPNamedArray('x', [1, 2, 3]) np.sin(n) # u'ひまわり'
とはいえ、上のように個々の関数について条件分岐を作成するのは手間だ。こんなとき、 ufunc 自体のプロパティを参照することで より簡単に書ける場合がある。
例えば、特定の型への処理をサポートする ufunc のみを利用する場合は .types
プロパティを参照すればよい。このプロパティは ufunc がどのような引数/返り値の型を取るかを示す Array-protocol Type String を返す。
参考
np.sin.types # ['e->e', 'f->f', 'd->d', 'g->g', 'F->F', 'D->D', 'G->G', 'O->O'] np.cos.types # ['e->e', 'f->f', 'd->d', 'g->g', 'F->F', 'D->D', 'G->G', 'O->O'] np.add.types # ['??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l', 'LL->L', # 'qq->q', 'QQ->Q', 'ee->e', 'ff->f', 'dd->d', 'gg->g', 'FF->F', 'DD->D', 'GG->G', # 'Mm->M', 'mm->m', 'mM->M', 'OO->O']
まとめ
NumPy 互換の動作を実現する方法のうち、Array Interface を利用した方法を記載した。自分で一から実装するよりは簡単だと思う。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (12件) を見る
Python でパイプ演算子を使いたい
ネタ記事です。/ This is a joke post which makes no practical sense.
はじめに
Python pandas
では主要な操作を以下のようにメソッドチェインの形で書くことができる。
# Python (pandas) df.assign(x=df['y'] + df['z']).groupby('x').sum()
pandas
v0.16.2 で DataFrame
と Series
に .pipe
というメソッドが追加され、このチェインを外部の関数/メソッドに対して連結できるようになった。利用例は以下のリンクを。
補足 matplotlib
でも v1.5.0 で ラベルデータ対応 が追加され、各関数が .pipe
から利用できるようになる予定。
このメソッドチェインによる "処理の連結" を、R ( {magrittr}
) のようにパイプ演算子 %>%
を使って統一的に書きたい。
# R (magrittr + dplyr) df %>% mutate(x = y + z) %>% group_by(x) %>% summarize_each(funs(sum))
R と異なり、 Python では自作の演算子を定義することはできない。Python の演算子の中では 右方向ビットシフト演算子 >>
がパイプ演算子っぽく見えなくもない。これを使って R のような動作をさせたい。
演算子のオーバーライド
まずはサンプルデータを作成する。
import numpy as np import pandas as pd pd.__version__ # '0.16.2' df = pd.DataFrame({'X': [1, 2, 3, 4, 5], 'Y': ['a', 'b', 'a', 'b', 'a']}) df # X Y # 0 1 a # 1 2 b # 2 3 a # 3 4 b # 4 5 a
普通に .pipe
を使うとこのような感じになる。
func = lambda x: x df.pipe(func) # X Y # 0 1 a # 1 2 b # 2 3 a # 3 4 b # 4 5 a
この処理を >>
演算子を使って書きたい。そのままでは、演算方法が未定義のため TypeError
となる。
df >> func
# TypeError: unsupported operand type(s) for >>: 'DataFrame' and 'function'
Python では各演算子に対応する処理はマジックメソッドとして実装する。そのため、>>
に対応する __rshift__
に .pipe
と同じ処理をさせればよい。また、R のように 右辺には関数を渡したいので、pandas
のメソッドを適当にラップする関数 pipify
を定義する。
pd.DataFrame.__rshift__ = pd.DataFrame.pipe def pipify(f): def wrapped(*args, **kwargs): return lambda self: f(self, *args, **kwargs) return wrapped head = pipify(pd.DataFrame.head)
できた。
df >> head(n=2) # X Y # 0 1 a # 1 2 b
groupby
も同じように書ける。集約関数は DataFrame
や DataFrameGroupBy
双方から使うため以下のように定義する。
groupby = pipify(pd.DataFrame.groupby) def sum(*args, **kwargs): return lambda self: self.sum(*args, **kwargs)
また、groupby
した結果からも >>
で演算を連結したい。そのためにDataFrameGroupBy.__rshift__
を定義する。
def _pipe(self, func): return func(self) pd.core.groupby.DataFrameGroupBy.__rshift__ = _pipe
できた。
df >> groupby('Y') >> sum() # X # Y # a 9 # b 6 df >> sum() # X 15 # Y ababa # dtype: object
右辺側クラスによるオーバーライド
DataFrame
が左辺にある場合、パイプ演算は上記のように定義できた。加えて、DataFrame
が右辺にある場合も同じようなことがしたい。
以下のような記法で、左辺の辞書から DataFrame
を作成したい。これも (当然) そのままではできない。
data = {'X': [1, 2, 3, 4, 5], 'Y': ['a', 'b', 'a', 'b', 'a']} data >> pd.DataFrame() # TypeError: unsupported operand type(s) for >>: 'dict' and 'DataFrame'
Python では、左辺側に 右辺を処理できるマジックメソッドがない場合は、右辺側のマジックメソッド ( DataFarme.__rrshift__
) が利用される。これを適当に定義すればよい。
def _get(a, b): return pd.DataFrame(b) pd.DataFrame.__rrshift__ = _get
できた。
data = {'X': [1, 2, 3, 4, 5], 'Y': ['a', 'b', 'a', 'b', 'a']} data >> pd.DataFrame() # X Y # 0 1 a # 1 2 b # 2 3 a # 3 4 b # 4 5 a
NumPy に対するオーバーライド
DataFrame
は np.ndarray
からも作成することができる。が、>>
ではうまく動かない。
data = np.array([[1, 2], [3, 4]]) pd.DataFrame(data) # 0 1 # 0 1 2 # 1 3 4 data >> pd.DataFrame() # TypeError: ufunc 'right_shift' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
ndarray
のマジックメソッドは スカラーや ndarray
ライクな引数を汎用的に処理できるように定義されている。具体的には、引数が ndarray
のサブクラス、もしくは Array Interface ( __array__
) をもつとき、ndarray
によって処理が行われる。
DataFrame
は ndarray
のサブクラスではないが Array Interface をもつ。右辺が不正なため通常の例外が送出されてしまい、右辺側のマジックメソッドに処理が移らない。
参考 Subclassing ndarray — NumPy v1.11 Manual
issubclass(pd.DataFrame, np.ndarray) # False df.__array__() # array([[1, 'a'], # [2, 'b'], # [3, 'a'], # [4, 'b'], # [5, 'a']], dtype=object)
ここで DataFrame
側のマジックメソッドを利用したい場合、上記リンクに記載されている __array_priority__
プロパティを設定すればよい。DataFrame
側の優先度を上げると期待した処理が行われることがわかる。
pd.DataFrame.__array_priority__ = 1000 data >> pd.DataFrame() # 0 1 # 0 1 2 # 1 3 4
まとめ
R 使ったほうがいい。Use R!
補足 こんなパッケージもあります。
2016/10/18 追記
続き(?)です。
Python XGBoost の変数重要度プロット / 可視化の実装
Gradient Boosting Decision Tree の C++ 実装 & 各言語のバインディングである XGBoost
、かなり強いらしいという話は伺っていたのだが自分で使ったことはなかった。こちらの記事で Python 版の使い方が記載されていたので試してみた。
その際、Python でのプロット / 可視化の実装がなかったためプルリクを出した。無事 マージ & リリースされたのでその使い方を書きたい。まずはデータを準備し学習を行う。
import numpy as np import xgboost as xgb from sklearn import datasets import matplotlib.pyplot as plt plt.style.use('ggplot') xgb.__version__ # '0.4' iris = datasets.load_iris() dm = xgb.DMatrix(iris.data, label=iris.target) np.random.seed(1) params={'objective': 'multi:softprob', 'eval_metric': 'mlogloss', 'eta': 0.3, 'num_class': 3} bst = xgb.train(params, dm, num_boost_round=18)
1. 変数重要度のプロット
Python 側には R のように importance matrix を返す関数がない。GitHub 上でも F score を見ろ、という回答がされていたので F score をそのままプロットするようにした。
xgb.plot_importance(bst)
棒グラフの色、タイトル/軸のラベルは以下のように変更できる。
xgb.plot_importance(bst, color='red', title='title', xlabel='x', ylabel='y')
color
にリストを渡せば棒ごとに色が変わる。色の順序は matplotlib
の barh
と同じく下からになる。また、ラベルを消したい場合は None
を渡す。
xgb.plot_importance(bst, color=['r', 'r', 'b', 'b'], title=None, xlabel=None, ylabel=None)
XGBoost
は内部的に変数名を保持していない。変数名でプロットしたい場合は 一度 F score を含む辞書を取得して、キーを差し替えてからプロットする。
bst.get_fscore() # {'f0': 17, 'f1': 16, 'f2': 95, 'f3': 59} iris.feature_names # ['sepal length (cm)', # 'sepal width (cm)', # 'petal length (cm)', # 'petal width (cm)'] mapper = {'f{0}'.format(i): v for i, v in enumerate(iris.feature_names)} mapped = {mapper[k]: v for k, v in bst.get_fscore().items()} mapped # {'petal length (cm)': 95, # 'petal width (cm)': 59, # 'sepal length (cm)': 17, # 'sepal width (cm)': 16} xgb.plot_importance(mapped)
2. 決定木のプロット
以下二つの関数を追加した。graphviz
が必要なためインストールしておくこと。
to_graphviz
: 任意の決定木をgraphviz
インスタンスに変換する。IPython
上であればそのまま描画できる。plot_tree
:to_graphviz
で取得したgraphviz
インスタンスをmatplotlib
のAxes
上に描画する。
IPython
から実行する。num_trees
で指定した番号に対応する木が描画される。
xgb.to_graphviz(bst, num_trees=1)
エッジの色分けが不要なら明示的に黒を指定する。
xgb.to_graphviz(bst, num_trees=2, yes_color='#000000', no_color='#000000')
IPython
を使っていない場合や、サブプロットにしたい場合には plot_tree
を利用する。
_, axes = plt.subplots(1, 2) xgb.plot_tree(bst, num_trees=2, ax=axes[0]) xgb.plot_tree(bst, num_trees=3, ax=axes[1])
何かおかしいことをやっていたら 本体の方で issue お願いします。
10/3追記 その後の修正を以下にしました。変数名の指定などが簡単になっています。