Theano で Deep Learning <2> : 多層パーセプトロン
Python Theano を使って Deep Learning の理論とアルゴリズムを学ぶ会、第二回。
目次
DeepLearning 0.1 より、
| 第一回 MNIST データをロジスティック回帰で判別する | 英 |
| 第二回 多層パーセプトロン (今回) | 英 |
| 第三回 畳み込みニューラルネットワーク | 英 |
| 第四回 Denoising オートエンコーダ | 英 |
| 第五回 多層 Denoising オートエンコーダ | 英 |
| 第六回の準備1 networkx でマルコフ確率場 / 確率伝搬法を実装する | - |
| 第六回の準備2 ホップフィールドネットワーク | - |
| 第六回 制約付きボルツマンマシン | 英 |
| Deep Belief Networks | 英 |
| Hybrid Monte-Carlo Sampling | 英 |
| Recurrent Neural Networks with Word Embeddings | 英 |
| LSTM Networks for Sentiment Analysis | 英 |
| Modeling and generating sequences of polyphonic music with the RNN-RBM | 英 |
多層パーセプトロンとは
ロジスティック回帰では判別関数が直線のため、入力データが線形分離できない場合は判別率が落ちていた。
多層パーセプトロンでは、入力を非線形変換/次元変換することによってデータを分離しやすくし、そのような場合の判別率を上げようとするもの。
モデル
多層パーセプトロンを DeepLearning 0.1 Multilayer Perceptron の式をベースに図示するとこんな感じになる。
: 入力層 - 隠れ層の間で適用される係数行列。次元は 入力データの説明変数の数
x 隠れ層のユニット数
: 入力層 - 隠れ層の間で適用される重みベクトル。次元は 隠れ層のユニット数
: 隠れ層の活性化関数。シグモイド関数もしくは
(ハイパボリックタンジェント) をよく使う。
: 隠れ層 - 出力層の間で適用される係数行列 = ロジスティック回帰の係数行列。次元は 隠れ層のユニット数 x 出力データのクラス数
: 隠れ層 - 出力層の間で適用される重みベクトル = ロジスティック回帰の重みベクトル。次元は 出力データのクラス数
: 出力層の活性化関数。2クラスの場合は シグモイド、多クラスの場合はソフトマックス関数 (ロジスティック回帰と同じ)。
つまり 3層パーセプトロンでは、
また、多層パーセプトロンの各層のユニット数 = その層に渡ってくるデータの次元 と考えればよい。
補足 PRML ではバイアス項をダミーのユニットを作って扱っているため、ユニット数 / 係数行列の次元など定義が異なる。
Pattern Recognition and Machine Learning (Information Science and Statistics)
- 作者: Christopher Bishop
- 出版社/メーカー: Springer
- 発売日: 2010/02/15
- メディア: ハードカバー
- 購入: 5人 クリック: 67回
- この商品を含むブログ (29件) を見る
ロジスティック回帰から多層パーセプトロンへ
まず、隠れ層をあらわす HiddenLayer クラスを定義する。
係数行列の初期値の決め方
前回の LogisticRegression クラスと同じく、HiddenLayer クラスは 係数行列 W ( 上式 ) , バイアス
b ( 上式 ) を
Theano の共有変数として宣言する。
係数行列 W の初期値は以下の範囲の一様乱数からとると (特にネットワークの層が深い場合に) 学習の進みがよいらしい 。
活性化関数が
tanh:
活性化関数が
sigmoid:
※ ,
は、それぞれ係数行列の前層 / 後続層のユニット数
そのため、HiddenLayer クラスでは 指定範囲の一様乱数を生成する numpy.random.uniform を使って係数行列 W の初期値を設定している。宣言中で使われている以下の変数は HiddenLayer.__init__ の引数として渡ってきているもの。
rng: 乱数のシードn_in: 入力層のユニット数n_out: 隠れ層のユニット数
W_values = numpy.asarray(
rng.uniform(
low=-numpy.sqrt(6. / (n_in + n_out)),
high=numpy.sqrt(6. / (n_in + n_out)),
size=(n_in, n_out)
),
dtype=theano.config.floatX
)
if activation == theano.tensor.nnet.sigmoid:
W_values *= 4
W = theano.shared(value=W_values, name='W', borrow=True)
例えば、前層のユニット数が 3, 後続層のユニット数が 4 で 活性化関数が tanh の場合、係数行列の初期値は以下のようになる。
import numpy n_in = 3 n_out = 4 numpy.random.uniform( low=-numpy.sqrt(6. / (n_in + n_out)), high=numpy.sqrt(6. / (n_in + n_out)), size=(n_in, n_out) ) # [[ 0.01572058 -0.46004291 -0.23671422 -0.54857948] # [ 0.081257 0.87090393 -0.72051252 -0.38176917] # [-0.09591302 0.04587085 -0.42426934 -0.35021306]]
また、バイアス b の初期値は 0ベクトルでよいようだ。
b_values = numpy.zeros((n_out,), dtype=theano.config.floatX)
b = theano.shared(value=b_values, name='b', borrow=True)
続けて、隠れ層からの出力を計算する式を HiddenLayer.output プロパティとして定義している。宣言中で使われている以下の変数は HiddenLayer.__init__ の引数として渡ってきているもの。
input: 入力データをあらわすシンボル (TensorVariable型)。activation: 活性化関数。デフォルトはtheano.tensor.tanh(np.tanhのテンソル版 )。
また、HiddenLayer はこの時点では実際のデータを受け取っておらず、Theano のシンボルから式を作っているだけ。ここも考え方は 前回同様。
lin_output = T.dot(input, self.W) + self.b self.output = ( lin_output if activation is None else activation(lin_output) )
MLP クラスの定義
続けて 多層パーセプトロン自体をあらわすクラス MLP を定義している。このクラスは 以下二つのインスタンスをプロパティとして持つ。
正則化 (Regularization)
MLP クラスでは 学習の際に 負の対数尤度だけでなく、以下 二つのノルムを正則化項として最小化するように設定している。一般には、これらを罰則 ( Penalty ) として最小化することで過学習をさけ、学習器の汎化性能 (未知のデータへのあてはまりの度合い) を上げることができる。
- L1 ノルム : 係数行列
- L2 ノルム : 係数行列
# L1 ノルム self.L1 = ( abs(self.hiddenLayer.W).sum() + abs(self.logRegressionLayer.W).sum() ) # L2 ノルム self.L2_sqr = ( (self.hiddenLayer.W ** 2).sum() + (self.logRegressionLayer.W ** 2).sum() )
負の対数尤度と判別誤差
他、モデル自体の負の対数尤度 ( negative_log_likelihood ) / 判別誤差 ( errors ) はモデルの出力に対して計算すればよいため、出力層である LogisticRegression での計算結果をそのまま使っている。
self.negative_log_likelihood = (
self.logRegressionLayer.negative_log_likelihood
)
self.errors = self.logRegressionLayer.errors
損失関数の定義
これら 負の対数尤度と正則化項を足し合わせて損失関数を計算する式 cost を作っている。L1_reg, L2_reg はそれぞれ L1ノルム, L2ノルムの正則化項の重み。既定値は L1_reg=0.00, L2_reg=0.0001 のため、サンプルをそのまま実行した場合は L2 ノルムのみが正則化項として考慮されることになる。
cost = (
classifier.negative_log_likelihood(y)
+ L1_reg * classifier.L1
+ L2_reg * classifier.L2_sqr
)
勾配の計算
ここが今回のポイント。多層パーセプトロンの学習は 誤差逆伝播法 ( Back propagation) というアルゴリズムで行われ、一般には以下のような処理 (順逆モデリング) が必要になる。
- 学習器の出力値と、真の値との差異 = 出力誤差を計算
- 1 で求めた出力誤差を小さくする
- 1 で求めた出力誤差に対して "隠れ層 - 出力層間の処理の逆演算"を行い、"出力誤差が隠れ層からでたと仮定した場合の出力値 = 隠れ層誤差" を求める
- 3 で求めた隠れ層誤差を小さくする
、、、めんどくさい、、、。
が、Theano では式表現に対して偏微分を行い偏導関数を求めることができる。そのため、上で定義した損失関数 cost について、最適化したいパラメータ HiddenLayer.W, HiddenLayer.b, LogisticRegression.W, LogisticRegression.b それぞれで偏微分した偏導関数を求めておけば、学習の際はそれらを使ってパラメータごとに勾配計算 / 更新をかけていくことができる。
つまり、今後 パーセプトロンの層 / パラメータが増えた場合でも、モデル全体の損失関数を偏微分することによって各パラメータを個別に学習できることになる。これは強力だな、、、。
classifier.params:HiddenLayer.W,HiddenLayer.b,LogisticRegression.W,LogisticRegression.bからなるリストlearning_rate: 学習率。デフォルトは 0.01
gparams = [T.grad(cost, param) for param in classifier.params] updates = [ (param, param - learning_rate * gparam) for param, gparam in zip(classifier.params, gparams) ]
続けて定義されている theano.function の読み方は 前回と同じ。以下 一連の処理を行う関数を作っている。
indexを引数として受け取り、givensの指定によりindexの位置のデータを切り出して シンボルx,yに入れ、x,yをoutputsの式 (負の対数尤度 + 正則化項) に与えて計算結果を得て、updatesで指定された更新式 ( 損失関数の偏導関数からの勾配計算 ) によって共有変数 ( 各係数行列、バイアス ) を更新する
全部まとめて
プログラムの量としては Theano のおかげで、え？こんだけでいいの？という感じだ。最後に置いてあるスクリプトを前回と同じディレクトリに入れて実行すればよい。
パラメータ調整のコツ
また、元文書では補足的に以下のパラメータ調整についてカッコ内のようなことが記載されている。自分でパラメータ調整したい方は読むといいですね。
- 隠れ層の活性化関数 (
tanhがいい感じ ) - 係数行列の初期値 (上の内容と同じ)
- 学習率 (
で試せ、時間で減衰させてみろ )
- 隠れ層のユニット数 ( データによる。モデルが複雑な場合は増やせ )
隠れ層からの出力の可視化
最後。元文書にはない内容だが、今回のサンプルについて 入力が隠れ層においてどういった感じで分離されるのか、以下の方法で可視化してみた。
- 学習データ (
train_set_x,train_set_y) を対象 - MNIST のラベル 0 〜 9 について、隣り合う 9組 ( 0と1, 1と2 ... ) を比較
- 隠れ層からの出力データ 500 次元に主成分分析をかけて 2次元に写像
- ランダムサンプリングした 300 レコードを描画
入力 - 隠れ層 間の変換によって各クラスが 主成分 ( = 分散大 ) の方向に離れていけば、分離超平面がみえるはずだ。
学習開始前 (隠れ層出力)
50 エポック学習後 (隠れ層出力)
曇りのない目で見てみると、なんか少し分離しやすそうになったかも？といえなくもない感じですね！
プログラム
プロット部分のみ gist に置いた ( 他は元文書のスクリプトと同じなので )。
Theano に関してのポイントは 1 点だけ。隠れ層からの出力を計算する HiddenLayer.output は TensorVariable 型のため、実際に計算するには一度 theano.function を通すこと。
# 隠れ層の出力 z_data を計算 apply_hidden = theano.function(inputs=[x_symbol], outputs=classifier.hiddenLayer.output) z_data = apply_hidden(x_data.get_value()) labels = y_data.eval()
Visualize Multilayer Perceptron Example in deeplearning.net · GitHub
11/30追記: 隠れ層での変換についてもう少しわかりやすい例をつくった。
12/07追記 続きはこちら。
- 作者: 岡谷貴之
- 出版社/メーカー: 講談社
- 発売日: 2015/04/08
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (4件) を見る
Python pandas の算術演算 / 集約関数 / 統計関数まとめ
概要
恒例の pandas 記事。今回は 基本的な算術演算についてまとめた。このあたりの挙動は numpy と一緒で直感的だと思うが、知っていないとハマるポイントがいくつかあるので。
準備
サンプルは DataFrame のみ。だが内容は Series でも同じ ( 行/列 2次元のデータに関するに記載は除く)。
import numpy as np import pandas as pd df = pd.DataFrame({'C1': [1, 1, 1], 'C2': [1, 1, 1], 'C3': [1, 1, 1]}) df # C1 C2 C3 # 0 1 1 1 # 1 1 1 1 # 2 1 1 1
四則演算
まずは基本的な例。
# スカラーの加算はデータ全体に対して適用 ( ブロードキャスト ) される df + 1 # C1 C2 C3 # 0 2 2 2 # 1 2 2 2 # 2 2 2 2 # np.array の加算は 各行に対する加算 df + np.array([1, 2, 3]) # C1 C2 C3 # 0 2 3 4 # 1 2 3 4 # 2 2 3 4 # 要素の長さが違う場合は ValueError df + np.array([1, 2]) # ValueError: Wrong number of items passed 2, place
とはいえ、演算をデータ全体へブロードキャストさせたいことはマレだと思う。特定の列 あるいは 特定セルに対して演算を行う場合は、以下の記事でまとめたデータ選択記法を使って計算結果を代入する。
補足 データへの演算自体は非破壊的な処理だが、代入によって元のデータが変更されることに注意。
# 1列目に [1, 2, 3] を足す df['C1'] = df['C1'] + np.array([1, 2, 3]) df # C1 C2 C3 # 0 2 1 1 # 1 3 1 1 # 2 4 1 1 # 3行目 / 3列目の値に 5 を足す df.iloc[2, 2] += 5 df # C1 C2 C3 # 0 2 1 1 # 1 3 1 1 # 2 4 1 6 # 複数列に対する演算も OK # (裏側では DataFrame 選択 -> 計算 -> 代入 をしているだけなので) df[['C1', 'C2']] -= 5 df # C1 C2 C3 # 0 -3 -4 1 # 1 -2 -4 1 # 2 -1 -4 6
算術演算メソッド
演算の挙動をもうすこし制御したい場合は、DataFrame.add。DataFrame や Series は、add のほかにも Python 標準モジュールの operators と対応する算術演算メソッド / 論理演算メソッドを持つ。利用できるメソッドの一覧はこちら。オプションはどれも重要なので順番に。
列に対するブロードキャスト (axis=0)
上記のような1列の加算ではなく、列全体に対してブロードキャストをしたい場合は axis=0 を指定。
axis=0もしくはaxis='index'で列に対する演算axis=1もしくはaxis='columns'で行に対する演算 (デフォルト)
# データは元に戻す df = pd.DataFrame({'C1': [1, 1, 1], 'C2': [1, 1, 1], 'C3': [1, 1, 1]}) df # C1 C2 C3 # 0 1 1 1 # 1 1 1 1 # 2 1 1 1 df.add(np.array([1, 2, 3]), axis=0) # C1 C2 C3 # 0 2 2 2 # 1 3 3 3 # 2 4 4 4 df.add(np.array([1, 2, 3]), axis='index') # 略
欠測値 ( NaN ) 要素へのパディング (fill_value)
データに 欠測値 ( NaN ) が含まれるとき、その要素への演算の結果も NaN になる。これは numpy の挙動と同じ。
# numpy の挙動 np.nan + 1 # nan # NaN を含む DataFrame を定義する df_nan = pd.DataFrame({'C1': [1, np.nan, np.nan], 'C2': [np.nan, 1, np.nan], 'C3': [np.nan, np.nan, 1]}) df_nan # C1 C2 C3 # 0 1 NaN NaN # 1 NaN 1 NaN # 2 NaN NaN 1 # NaN を含むセルの演算結果は NaN df.add(df_nan) # C1 C2 C3 # 0 2 NaN NaN # 1 NaN 2 NaN # 2 NaN NaN 2
これに気づかず四則演算していると、なんか合計があわないな？ということになってしまう。この挙動を変えたい場合は fill_value。演算の実行前に fill_value で指定した値がパディングされる。ただし、演算対象の要素が 両方とも NaN の場合は NaN のまま。
df.add(df_nan, fill_value=0) # C1 C2 C3 # 0 2 1 1 # 1 1 2 1 # 2 1 1 2 # データの順序が変わっても有効 ( fill_value は演算対象 両方のデータに適用される) df_nan.add(df, fill_value=0) # C1 C2 C3 # 0 2 1 1 # 1 1 2 1 # 2 1 1 2 # 要素が 両方とも NaN の場合はパディングされない df_nan.add(df_nan, fill_value=0) # C1 C2 C3 # 0 2 NaN NaN # 1 NaN 2 NaN # 2 NaN NaN 2 # DataFrame について、対象データが array, Series の場合の処理は未実装 (NotImplementedError) なので注意 df.add(np.array([1, np.nan, 1]), fill_value=0) # NotImplementedError: fill_value 0 not supported
最後のオプション、 level は index が複数の階層 (レベル) 持つ場合 ( MultiIndex を持つ場合 ) の制御を行う。このオプションの挙動は近日公開予定 (?) の MultiIndex 編で。
行列積
* で DataFrame 同士の積をとった場合は、各要素同士の積になる。
df2 = pd.DataFrame({'C1': [2, 3],
'C2': [4, 5]})
df2
# C1 C2
# 0 2 4
# 1 3 5
df3 = pd.DataFrame({'C1': [1, -2],
'C2': [-3, 4]})
df3
# C1 C2
# 0 1 -3
# 1 -2 4
df2 * df3
# C1 C2
# 0 2 -12
# 1 -6 20
行列の積をとりたい場合は DataFrame.dot。ただし、行列の積をとるためには元データの columns と 引数の index のラベルが一致している必要がある (説明 次節)。そのため、上のように columns / index が一致しないデータから直接 行列の積をとることはできない。そんなときは DataFrame.values プロパティで引数側の内部データを numpy.array として取り出せばよい。
# NG! df2.dot(df3) # ValueError: matrices are not aligned # values プロパティで内部データを numpy.array として取り出せる df3.values # array([[ 1, -3], # [-2, 4]], dtype=int64) # OK! df2.dot(df3.values) # 0 1 # 0 -6 10 # 1 -7 11
引数側を転置すれば、元データの columns と 引数の index のラベルが一致する。このようなデータはそのまま DataFrame.dot で積をとれる。
df3.T # 0 1 # C1 1 -2 # C2 -3 4 df2.dot(df3.T) # 0 1 # 0 -10 12 # 1 -12 14
補足 そもそも Python には行列積の演算子がない。が、現在 PEP-465 で行列積として @ 演算子が提案されている。
ラベルによる演算の制御
ここまでは index, columns が一致するデータのみを対象にしてきた。実際には演算を行うデータの index, columns が異なることもある。その場合の挙動には少し注意が必要。
Series, DataFrame 同士の演算は、 index, columns のラベルが一致する要素同士で行われる。例えば以下のように index と columns がずれているとき、対応しない要素は NaN となる。
df4 = pd.DataFrame({'C1': [1, 1, 1],
'C2': [1, 1, 1],
'C3': [1, 1, 1]})
df4
# C1 C2 C3
# 0 1 1 1
# 1 1 1 1
# 2 1 1 1
df5 = pd.DataFrame({'C2': [1, 1, 1],
'C3': [1, 1, 1],
'C4': [1, 1, 1]},
index=[1, 2, 3])
df5
# C2 C3 C4
# 1 1 1 1
# 2 1 1 1
# 3 1 1 1
df4 + df5
# C1 C2 C3 C4
# 0 NaN NaN NaN NaN
# 1 NaN 2 2 NaN
# 2 NaN 2 2 NaN
# 3 NaN NaN NaN NaN
どう直せばよいかは状況による。index, columns のずれ自体は OK で、 NaN でのパディングが気に入らない場合は fill_value。
df4.add(df5, fill_value=0) # C1 C2 C3 C4 # 0 1 1 1 NaN # 1 1 2 2 1 # 2 1 2 2 1 # 3 NaN 1 1 1
index, columnsに関係なく、対応する位置の要素同士で演算したい場合は DataFrame.values を使って取得した numpy.array のみを渡せばよい。
df4 + df5.values # C1 C2 C3 # 0 2 2 2 # 1 2 2 2 # 2 2 2 2
もしくは、以下のようにして データの index もしくは columns を揃えてから演算する。
df4のindexをdf5にそろえる場合は、df4.indexに代入。df5のindexをdf4にそろえる場合は、df5.indexに代入。もしくは、DataFrame.reset_index(drop=True)でindexを 0 から振りなおし
# df4 の index を df5 にそろえる df4.index = [1, 2, 3] df4 # C1 C2 C3 # 1 1 1 1 # 2 1 1 1 # 3 1 1 1 # df5 の index を df4 にそろえる df5.reset_index(drop=True) # C2 C3 C4 # 0 1 1 1 # 1 1 1 1 # 2 1 1 1 # reset_index デフォルトでは、元の index をカラムとして残す df5.reset_index() # index C2 C3 C4 # 0 1 1 1 1 # 1 2 1 1 1 # 2 3 1 1 1
集約関数
Series, DataFrame は 合計を取得する sum, 平均を計算する mean などの一連の集約関数に対応するメソッドを持つ。一覧は こちら。
df6 = pd.DataFrame({'C1': ['A', 'B', 'C'],
'C2': [1, 2, 3],
'C3': [4, 5, 6]})
df6
# C1 C2 C3
# 0 A 1 4
# 1 B 2 5
# 2 C 3 6
集約系の関数は既定では列方向に対して適用される。基本的には 数値型のみが集約対象になる。が、一部の関数では数値以外の型に対してもがんばってくれるものもある。
# mean は数値型のみ集約 df6.mean() # C2 2 # C3 5 # dtype: float64 # sum は 文字列も集約 df6.sum() # C1 ABC # C2 6 # C3 15 # dtype: object # 数値以外が不要な場合は numeric_only=True df6.sum(numeric_only=True) # C2 6 # C3 15 # dtype: int64 # 行方向へ適用したい場合は axis = 1 # このとき、数値型以外は除外して関数適用される df6.sum(axis=1) # 0 5 # 1 7 # 2 9 # dtype: int64 # 明示的に含めようとするとエラー df6.sum(axis=1, numeric_only=False) # TypeError: cannot concatenate 'str' and 'long' objects
複数列のデータをまとめて関数適用したい、という場合もたまーにある。pandas には直接対応するメソッドはないが、以下のようにすればできる。
# values.flatten() で 2列分の値を ひとつの numpy.array として取得 df6[['C2', 'C3']].values.flatten() # array([1, 4, 2, 5, 3, 6], dtype=int64) # 集約関数適用 np.mean(df6[['C2', 'C3']].values.flatten()) # 3.5
補足 numpy.array.flatten() は1次元化した array のコピーを返す。
統計関数
また、分散を計算する var, 標準偏差を計算する std などの統計関数に対応するメソッドもある。一覧は上と同じく こちら。関数の適用方向など、挙動は集約関数と一緒。
1点 覚えておくべき箇所は、pandas では分散 / 標準偏差について不偏推定量の計算がデフォルトになっている。これは numpy の挙動 ( 標本統計量を返す ) とは異なる。この挙動は pandas, numpy ともに ddof オプションで制御できる。
pandas: 不偏推定量の計算 (ddof=True) がデフォルト。numpy: 標本統計量の計算 (ddof=False) がデフォルト。
補足 ddof = Delta Degrees of Freedom
df7 = pd.DataFrame({'C1': [1, 2, 3, 4],
'C2': [4, 5, 6, 7],
'C3': [2, 3, 3, 2]})
# C1 C2 C3
# 0 1 4 2
# 1 2 5 3
# 2 3 6 3
# 3 4 7 2
# 不偏分散
df7.var()
# C1 1.666667
# C2 1.666667
# C3 0.333333
# dtype: float64
# 標本分散
df7.var(ddof=False)
# C1 1.25
# C2 1.25
# C3 0.25
# dtype: float64
# 標本分散 (numpy)
np.var(df7)
# C1 1.25
# C2 1.25
# C3 0.25
# dtype: float64
# 不偏標準偏差
df7.std()
# C1 1.290994
# C2 1.290994
# C3 0.577350
# dtype: float64
# 標本標準偏差
df7.std(ddof=False)
# C1 1.118034
# C2 1.118034
# C3 0.500000
# dtype: float64
# 標本標準偏差 (numpy)
np.std(df7)
# C1 1.118034
# C2 1.118034
# C3 0.500000
# dtype: float64
基本統計量の表示 ( describe )
最後。基本統計量をまとめて計算したい場合は DataFrame.describe()。
df7.describe() # C1 C2 C3 # count 4.000000 4.000000 4.00000 # mean 2.500000 5.500000 2.50000 # std 1.290994 1.290994 0.57735 # min 1.000000 4.000000 2.00000 # 25% 1.750000 4.750000 2.00000 # 50% 2.500000 5.500000 2.50000 # 75% 3.250000 6.250000 3.00000 # max 4.000000 7.000000 3.00000
まとめ
pandas での算術演算、集約関数、統計関数の挙動をまとめた。ポイントは、
DataFrameへの演算は適用される方向に注意。演算方向を指定する場合、列方向ならaxis=0, 行方向はaxis=1。- 演算したい要素に
NaNが含まれる場合、必要に応じてfill_value。 - 演算したい要素同士のラベルは一致させること。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (9件) を見る
Theano で Deep Learning <1> : MNIST データをロジスティック回帰で判別する
概要
ここ数年 Deep Learning 勢の隆盛いちじるしい。自分が学生の頃は ニューラルネットワークはオワコン扱いだったのに、、、どうしてこうなった？自分もちょっと触ってみようかな、と記事やらスライドやら読んでみても、活性化関数が〜 とか、 制約付き何とかマシンが〜（聞き取れず。何か中ボスっぽい名前）とか、何言っているのかよくわからん。
巷には 中身がわかっていなくてもある程度 使えるパッケージもいくつかあるようだが、せっかくなので少しは勉強したい。
Python 使って できんかな？と思って探してみると、すでに Theano というPython パッケージの開発チームが作った DeepLearning Documentation 0.1 という大部の聖典 (バイブル) があった。
当然だがこの文書では Theano の機能をいろいろ使っているため、ぱっと見では 何をやってんだかよくわからない。ということで、Theano の仕様を調べたり numpy で翻訳したりしながら読み解いてみることにした。
聖典の目次
かなり楽観的な見通しだが、 1章/週 のペースで進めても 3ヶ月近くかかるな、、。
| 第一回 MNIST データをロジスティック回帰で判別する (今回) | 英 |
| 第二回 多層パーセプトロン | 英 |
| 第三回 畳み込みニューラルネットワーク | 英 |
| 第四回 Denoising オートエンコーダ | 英 |
| 第五回 多層 Denoising オートエンコーダ | 英 |
| 第六回の準備1 networkx でマルコフ確率場 / 確率伝搬法を実装する | - |
| 第六回の準備2 ホップフィールドネットワーク | - |
| 第六回 制約付きボルツマンマシン | 英 |
| Deep Belief Networks | 英 |
| Hybrid Monte-Carlo Sampling | 英 |
| Recurrent Neural Networks with Word Embeddings | 英 |
| LSTM Networks for Sentiment Analysis | 英 |
| Modeling and generating sequences of polyphonic music with the RNN-RBM | 英 |
全体の流れ / スクリプトは元文書を参照。以下の各セクション名には対応するセクションへのリンクを貼っている。
Theano とは
こんなことができます。
- シンボル(変数名) / 値 (スカラー, ベクトル, 行列) / 式を 数式表現そのもの (構文木) として扱える。
- 式の定義 / 途中計算結果の再利用が楽
- 機械学習で頻出する勾配計算を解析的に実行できる
- 高速な計算実行
インストール
pip install theano
準備
import numpy import theano import theano tensor as T
モデル
ロジスティック回帰が 入力をどのようにクラス分類するか？という方法が書いてある。ロジスティック回帰 そのものについてはこちらでまとめた。
プログラム中で まず 目に付くのは以下のような theano.shared 表現。Theano では以下 二つのメモリ空間を明示的にわけて扱う。
- Python 空間 ( user space ) : Python スクリプト上で定義された普通の変数が存在する。
Theano空間 :theano.sharedによって定義された変数 (以降、共有変数) が存在する。物理的に Python 空間とは別に確保され、Theanoの式表現(後述)から参照できる。また、GPU へのデータ転送を高速にする。
theano.shared で宣言された変数は 実データとして numpy のデータをもつ。
# この表現で作られるデータは W = theano.shared( value=numpy.zeros( (n_in, n_out), dtype=theano.config.floatX ), name='W', borrow=True ) # numpy だとこちらと同じ W = numpy.zeros((n_in, n_out), dtype=numpy.float64)
theano.shared による共有変数の定義
上でも指定されている引数 borrow は Python 空間上で定義されたデータの実体を 共有変数でも共有するかどうかを決める。
# Python 空間で変数定義 x = numpy.array([1, 1, 1]) # theano.shared で 共有変数化 xs = theano.shared(x) # 返り値は int32 / vector の TensorSharedVariable 型 xs # <TensorType(int32, vector)> type(xs) # theano.tensor.sharedvar.TensorSharedVariable # 内部の値を得るためには get_value xs.get_value() # array([1, 1, 1]) # デフォルトでは 対象オブジェクトのコピーが 共有変数の実体になるため、 # Python 空間の変数を変更しても 共有変数の実体には影響しない x[1] = 2 x # array([1, 2, 1]) xs.get_value() # array([1, 1, 1]) # 同じ実体を共有したい場合は borrow=True xs = theano.shared(x, borrow=True) xs.get_value() # array([1, 2, 1]) # Python 空間での変更が 共有変数へも反映される x[1] = 3 x # array([1, 3, 1]) xs.get_value() # array([1, 3, 1])
theano.tensor 上の関数の呼び出し
Theano では シンボル(変数名) / スカラー / ベクトル / 行列などはテンソル型として抽象化される。テンソル型に対する処理は theano.tensor 上で定義された関数を使って行う。
例えば、多項ロジスティック回帰では 入力があるクラスに分類される確率はソフトマックス関数を用いて書けた。Theano にはソフトマックス関数を計算する theano.tensor.nnet.softmax があって、
d = theano.shared(value=numpy.matrix([[0.1, 0.2, 0.3], [0.3, 0.2, 0.1], [0.1, 0.2, 0.3]], dtype=theano.config.floatX), name='d', borrow=True) sm = T.nnet.softmax(d) # 返り値は TensorVariable 型 type(sm) # <class 'theano.tensor.var.TensorVariable'> # TensorVariable から結果を取り出す (関数を評価する) には eval() sm.eval() # [[ 0.30060961 0.33222499 0.3671654 ] # [ 0.3671654 0.33222499 0.30060961] # [ 0.30060961 0.33222499 0.3671654 ]]
ロジスティック回帰の予測値は上記の確率を最大にするラベルになる。値が最大となるラベルをとる処理は theano.tensor.argmax。これは numpy.argmax のテンソル版。以下の通り、theano.tensor の関数は通常 (共有変数化されていない) の numpy.array も受け取ることができる。その場合も返り値は TensorVariable 型になる。
# 3 番目の要素 ( index でみて 2 ) が最も大きい am = T.argmax(numpy.array([1, 2, 3, 2, 1])) type(am) # theano.tensor.var.TensorVariable am.eval() # array(2L, dtype=int64)
損失関数の定義
これはそのまま。すべて theano.tensor の関数であることには注意。
LogisticRegression のクラス作成
LogisticRegression.errors で誤差 (ここでの誤差は真のクラスと 判別結果の差)を取っている。論理演算 != についても テンソル版である theano.tensor.neq が定義されているのでそちらを使う。
# 真のクラス y = numpy.array([1, 0, 1, 1, 1]) # 予測値が2つ誤判別しているとする y_pred = numpy.array([1, 1, 1, 0, 1]) # theano.tensor T.mean(T.neq(y_pred, y)).eval() # 0.4 # numpy numpy.mean(y_pred != y) # 0.4
モデルの学習
ロジスティック回帰の学習について。勾配 (Gradient) の計算は Theano だと theano.grad を使って簡単にできるようだ。ここからの一連の処理が今回のキモ。
Theano での勾配の計算 (微分)
Theano での勾配の計算については公式ドキュメントに 詳細 がある。Theano では、定義された式にあらわれるシンボル / 定数を構文木として保持している。そのため、勾配計算の際は 構文木解析によって導関数を求めて実行することができる。とりあえず を定義して微分してみると、
# シンボル x を宣言 (この時点で x に値はない) x = T.dscalar('x') # y = x ** 2 という式を宣言 # シンボルに対して演算子から式を作っていけるところがカッコイイ y = x ** 2 # y を x について微分 (まだ値はない) gy = T.grad(y, x) # gy の定義 (導関数) を確認 theano.pp(gy) # '((fill((x ** TensorConstant{2}), TensorConstant{1.0}) * TensorConstant{2}) * (x ** (TensorConstant{2} - TensorConstant{1})))'
TensorConstant は {} 内の値をもつスカラー。また、fill は 第一引数を 第二引数で埋める、という意味になる。順番に数値を埋めてみると、、、おお、、、ちゃんと になっている。
# '((fill((x ** TensorConstant{2}), TensorConstant{1.0}) * TensorConstant{2}) * (x ** (TensorConstant{2} - TensorConstant{1})))' # = '((fill((x ** 2), 1.0) * 2) * (x ** (2 - 1)))' # = '(2 * x)'
続けて、式から関数を作成し、実際の値を渡して計算 (式の評価) を行う。え、式にそのまま引数渡せないの？と思うが、どうも function の実行時に式表現をコンパイルしているのでダメらしい。
# 式から関数を作る。function の第一引数は関数に与える引数, 第二引数に関数化する式 f = function([x], gy) # x = 4 のときの x ** 2 の傾きを求める f(4) # array(8.0) # x = 8 のときの x ** 2 の傾きを求める f(8) # array(16.0)
さて、元の文書にもどってロジスティック回帰を定義 / 学習を始めるところをみてみる。ここでのテンソル型の式から学習器を作っていく流れは格好いい。
# ミニバッチ確率勾配降下法で指定する index (スカラー)を入れるシンボル index = T.lscalar() # 入力データの行列を示すシンボル x = T.matrix('x') # ラベルのベクトルを示すシンボル y = T.ivector('y') # ロジスティック回帰の処理も theano.tensor の関数で定義されているため、 # シンボル x が実データなしで渡せる classifier = LogisticRegression(input=x, n_in=28 * 28, n_out=10) # 引数 y もシンボル。返り値 cost では負の対数尤度を計算する式が返ってくる cost = classifier.negative_log_likelihood(y) # w, b の勾配を計算する式 g_W = T.grad(cost=cost, wrt=classifier.W) g_b = T.grad(cost=cost, wrt=classifier.b)
で、最後にこれらを theano.functionで関数化する。引数の意味は、
inputs: 関数への入力 (引数)。outputs: 関数化される式。updates: 共有変数を更新する式。givens: 引数 (inputs) -> 共有変数へのマッピングを行う辞書。
つまり、
# w, b の更新を定義する式 updates = [(classifier.W, classifier.W - learning_rate * g_W), (classifier.b, classifier.b - learning_rate * g_b)] train_model = theano.function( inputs=[index], outputs=cost, updates=updates, givens={ x: train_set_x[index * batch_size: (index + 1) * batch_size], y: train_set_y[index * batch_size: (index + 1) * batch_size] } )
上の定義は
indexを引数として受け取り、givensの指定によりindexの位置のデータを切り出して シンボルx,yに入れ、x,yをoutputsの式 (負の対数尤度) に与えて計算結果を得て、updatesで指定された更新式によって共有変数を更新する
という一連の処理を行う関数を作っている。うーん、これは訓練されていないと ちゃっと読めないな、、、。
そして、この時点においても train_model では ロジスティック回帰の学習を 関数として定義しただけで、実計算はまだ行われていない。実際の学習 ( 続く while ループ中 ) では この定義に対して index だけを渡してやれば、 必要な値が計算 / 関連する値が自動的に更新されていく。
全部まとめて
以降は特になし。MNISTデータをローカルに保存した後、まとめ にあるスクリプトをコピペして実行すればよい。
11/30追記: つづきはこちら。
- 作者: 岡谷貴之
- 出版社/メーカー: 講談社
- 発売日: 2015/04/08
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (4件) を見る
ロジスティック回帰 (勾配降下法 / 確率的勾配降下法) を可視化する
いつの間にかシリーズ化して、今回はロジスティック回帰をやる。自分は行列計算ができないクラスタ所属なので、入力が3次元以上 / 出力が多クラスになるとちょっときつい。教科書を読んでいるときはなんかわかった感じになるんだが、式とか字面を追ってるだけだからな、、、やっぱり自分で手を動かさないとダメだ。
また、ちょっとした事情により今回は Python でやりたい。Python のわかりやすい実装ないんかな？と探していたら 以下の ipyton notebook を見つけた。
http://nbviewer.ipython.org/gist/mitmul/9283713
こちらのリンク先に2クラス/多クラスのロジスティック回帰 (確率的勾配降下法) のサンプルがある。ありがたいことです。理論的な説明も書いてあるので ロジスティック回帰って何？という方は上を読んでください (放り投げ)。
この記事ではロジスティック回帰で使われる勾配法について書きたい。
補足 この記事ではデータを pandas で扱う / matplotlib でアニメーションさせる都合上、元のプログラムとは処理/変数名などだいぶ変わっているため見比べる場合は注意。
補足 実際の分析で使いたい場合は scikit-learn の SGDClassifier 。
サンプルデータのロードと準備
iris で、"Species" を目的変数にする。
補足 rpy2 を設定している方は rpy2から、そうでない方は こちら から .csv でダウンロードして読み込み (もしくは read_csv のファイルパスとして直接 URL 指定しても読める)。
import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib.animation as animation # iris の読み込みはどちらかで # rpy2 経由で R から iris をロード # import pandas.rpy.common as com # iris = com.load_data('iris') # csv から読み込み # http://aima.cs.berkeley.edu/data/iris.csv names = ['Sepal.Length', 'Sepal.Width', 'Petal.Length', 'Petal.Width', 'Species'] iris = pd.read_csv('iris.csv', header=None, names=names) np.random.seed(1) # 描画領域のサイズ FIGSIZE = (5, 3.5)
ロジスティック回帰 (2クラス)
まず iris データの一部を切り出して、説明変数 2次元 / 目的変数 2クラスのデータにする。
# 2 クラスにするため、setosa, versicolor のデータのみ抽出 data = iris[:100] # 説明変数は 2つ = 2 次元 columns = ['Petal.Width', 'Petal.Length'] x = data[columns] # データ (説明変数) y = data['Species'] # ラベル (目的変数)
ロジスティック回帰での目的変数は 0, 1 のみからなるベクトルである必要があるので、以下のようにして変換。論理演算で bool の pd.Series を作り、int 型に変換すればよい。
# ラベルを0, 1の列に変換 y = (y == 'setosa').astype(int) y # 1 1 # 2 1 # ... # 99 0 # 100 0 # Name: Species, Length: 100, dtype: int64
プロットしてみる。線形分離できるのでサンプルとしては OK だろう。
def plot_x_by_y(x, y, colors, ax=None): if ax is None: # 描画領域を作成 fig = plt.figure(figsize=FIGSIZE) # 描画領域に Axes を追加、マージン調整 ax = fig.add_subplot(1, 1, 1) fig.subplots_adjust(bottom=0.15) x1 = x.columns[0] x2 = x.columns[1] for (species, group), c in zip(x.groupby(y), colors): ax = group.plot(kind='scatter', x=x1, y=x2, color=c, ax=ax, figsize=FIGSIZE) return ax plot_x_by_y(x, y, colors=['red', 'blue']) plt.show()
最適化法
ロジスティック回帰における予測値 は、
で計算される。それぞれの意味は、
: 真のラベル
の予測値ベクトル。次元は (入力データ数, 1)
: シグモイド関数。計算結果を (0, 1) 区間に写像する
: 入力データ。次元は (入力データ数, 説明変数の数)
: 係数ベクトル。次元は(クラス数 = 2)
: バイアス (スカラー)
ロジスティック回帰では、予測値 と真のラベル
の差から計算される損失関数を最小化する回帰直線を求めることになる。この回帰直線を求める方法として勾配法が知られている。勾配法とは、損失関数から 勾配 (Gradient) を求めることによって損失関数をより小さくする
,
を逐次更新していく手法 (詳細は上記リンク)。
ここでは以下 3 種類の勾配法を試す。上のリンクでやっているのは 3つ目の ミニバッチ確率的勾配降下法 / MSGD になる。
- 勾配降下法 (Gradient Descent)
- 確率的勾配降下法 (Stochastic Gradient Descent / SGD)
- ミニバッチ確率的勾配降下法 (Minibatch SGD / MSGD)
※ MSGD は SGD と特に区別されない場合も多い。
まず、共通で使う関数を定義する。
def p_y_given_x(x, w, b): # x, w, b から y の予測値 (yhat) を計算 def sigmoid(a): return 1.0 / (1.0 + np.exp(-a)) return sigmoid(np.dot(x, w) + b) def grad(x, y, w, b): # 現予測値から勾配を計算 error = y - p_y_given_x(x, w, b) w_grad = -np.mean(x.T * error, axis=1) b_grad = -np.mean(error) return w_grad, b_grad
勾配降下法 (Gradient Descent)
勾配法の中ではもっともシンプルな考え方で、全データ (x, y) を使って勾配を求める。実装は以下のようになる。
引数の意味は、
x,y,w,b: 上記数式に対応eta: 学習率 (勾配をどの程度 w, b に反映させるか)num: イテレーション回数
補足 num の既定値はそれぞれの手法がある程度収束する値にしている。ちゃんとやる場合は収束判定すべき。
def gd(x, y, w, b, eta=0.1, num=100): for i in range(1, num): # 入力をまとめて処理 w_grad, b_grad = grad(x, y, w, b) w -= eta * w_grad b -= eta * b_grad e = np.mean(np.abs(y - p_y_given_x(x, w, b))) yield i, w, b, e
実行してみる。w, b の初期値は全て 0 のベクトル / スカラーにした。
# w, b の初期値を作成 w, b = np.zeros(2), 0 gen = gd(x, y, w, b) # gen はジェネレータ gen # <generator object gd at 0x11108e5f0> # 以降、gen.next() を呼ぶたびに 一回 勾配降下法を実行して更新した結果を返す。 # タプルの中身は (イテレーション回数, w, b, 誤差) gen.next() # (1, array([-0.027 , -0.06995]), 0.0, 0.47227246182037463) gen.next() # (2, array([-0.04810306, -0.12007078]), 0.0054926687253766763, 0.45337584157628485) gen.next() # (3, array([-0.06522081, -0.15679859]), 0.014689105435764377, 0.44035208019270661) gen.next() # (4, array([-0.07963026, -0.18443466]), 0.026392079742058178, 0.43101983106700503) # ...
実行のたびに w, b が更新されるとともに、誤差の値が小さくなっていることがわかる。ここでの誤差 は 誤判別率 ではなく、真のクラスと クラス所属確率の差。
確率的勾配降下法 (Stochastic Gradient Descent / SGD)
勾配降下法では勾配の計算を全データを行列とみて行うため、データ量が増えると計算が厳しい。確率的勾配降下法では、入力をランダムに (確率的に) 並べ替えて、データ一行ずつから勾配を求めて係数/バイアスを更新する。一行ずつの処理になるため計算量の問題は解消する。実装は以下のようになる。
def sgd(x, y, w, b, eta=0.1, num=4): for i in range(1, num): for index in range(x.shape[0]): # 一行ずつ処理 _x = x.iloc[[index], ] _y = y.iloc[[index], ] w_grad, b_grad = grad(_x, _y, w, b) w -= eta * w_grad b -= eta * b_grad e = np.mean(np.abs(y - p_y_given_x(x, w, b))) yield i, w, b, e
確率的に処理するため、データをランダムに並べ替える。ランダムに並べ替えられた結果の表示から、x, y の index が一致している = 対応関係が維持されていることを確認。
# index と同じ長さの配列を作成し、ランダムにシャッフル indexer = np.arange(x.shape[0]) np.random.shuffle(indexer) # x, y を シャッフルされた index の順序に並び替え x = x.iloc[indexer, ] y = y.iloc[indexer, ] x.head() # Petal.Width Petal.Length # 81 1.1 3.8 # 85 1.5 4.5 # 34 0.2 1.4 # 82 1.0 3.7 # 94 1.0 3.3 y.head() # 81 0 # 85 0 # 34 1 # 82 0 # 94 0 # Name: Species, dtype: int64
実行してみよう。
# w, b の初期値を作成 w, b = np.zeros(2), 0 gen = sgd(x, y, w, b) # 以降、gen.next() を呼ぶたびに 一行分のデータで更新した結果を返す。 # タプルの中身は (イテレーション回数, w, b, 誤差) gen.next() # (1, array([-0.055, -0.19 ]), -0.050000000000000003, 0.43367318357380596) gen.next() # (1, array([-0.09571092, -0.31213277]), -0.07714061571720722, 0.4071586323043157) gen.next() # (1, array([-0.08310602, -0.22389845]), -0.014116100082250033, 0.42197958122210588)
gen.next によって 一行処理した結果が返ってくる。一行あたりの処理をみると損失は増加していることもある (が徐々に減っていく)。
補足 イテレーション回数は全データに対する回数のため、データが一回りするまで増えない。
ミニバッチ確率的勾配降下法 (Minibatch SGD / MSGD)
確率的勾配降下法について、一行ずつではなく 計算できる程度の固まり (Minibatch) 単位で処理していく方法。実装は以下のようになる。
※ 事前にデータをランダムに並べ替えるのは確率的勾配降下法と同様。
引数はこれまでから一つ増えて、
batch_size: Minibatch として切り出すデータ数 (行数)
def msgd(x, y, w, b, eta=0.1, num=25, batch_size=10): for i in range(1, num): for index in range(0, x.shape[0], batch_size): # index 番目のバッチを切り出して処理 _x = x[index:index + batch_size] _y = y[index:index + batch_size] w_grad, b_grad = grad(_x, _y, w, b) w -= eta * w_grad b -= eta * b_grad e = np.mean(np.abs(y - p_y_given_x(x, w, b))) yield i, w, b, e
実行。
# w, b の初期値を作成 w, b = np.zeros(2), 0 gen = msgd(x, y, w, b) # 読み方/補足は 確率的勾配降下法と同じ。 gen.next() # (1, array([ 0.011 , 0.0725]), 0.050000000000000003, 0.52633544830099133) gen.next() # (1, array([ 0.02251519, 0.13792904]), 0.096115503575706834, 0.54788728897754557)
可視化
以上 3 つの処理を見比べてわかるとおり、差異は勾配を計算するのに利用するデータの切り出し処理のみ。それぞれどんな動きをしているのか可視化してみる。
def plot_logreg(x, y, fitter, title): # 描画領域を作成 fig = plt.figure(figsize=FIGSIZE) # 描画領域に Axes を追加、マージン調整 ax = fig.add_subplot(1, 1, 1) fig.subplots_adjust(bottom=0.15) # 描画オブジェクト保存用 objs = [] # 回帰直線描画用の x 座標 bx = np.arange(x.iloc[:, 0].min(), x.iloc[:, 0].max(), 0.1) # w, b の初期値を作成 w, b = np.zeros(2), 0 # 勾配法の関数からジェネレータを生成 gen = fitter(x, y, w, b) # ジェネレータを実行し、勾配法 1ステップごとの結果を得る for i, w, b, e in gen: # 回帰直線の y 座標を計算 by = -b/w[1] - w[0]/w[1]*bx # 回帰直線を生成 l = ax.plot(bx, by, color='gray', linestyle='dashed') # 描画するテキストを生成 wt = """Iteration = {0} times w = [{1[0]:.2f}, {1[1]:.2f}] b = {2:.2f} error = {3:.3}""".format(i, w, b, e) # axes 上の相対座標 (0.1, 0.9) に text の上部を合わせて描画 t = ax.text(0.1, 0.9, wt, va='top', transform=ax.transAxes) # 描画した line, text をアニメーション用の配列に入れる objs.append(tuple(l) + (t, )) # データ, 表題を描画 ax = plot_x_by_y(x, y, colors=['red', 'blue'], ax=ax) ax.set_title(title) # アニメーション開始 ani = animation.ArtistAnimation(fig, objs, interval=1, repeat=False) plt.show() plot_logreg(x, y, gd, title='Gradient descent') plot_logreg(x, y, sgd, title='Stochastic gradient descent') plot_logreg(x, y, msgd, title='Minibatch SGD')
勾配降下法 (Gradient Descent)
勾配降下法では 全データから勾配をもとめるので、イテレーション時に損失が増えることはない。ただしデータ数が多くなると使えない。
確率的勾配降下法 (Stochastic Gradient Descent / SGD)
確率的勾配降下法では 一行ずつで勾配をもとめるので、処理される行によっては損失が増えることもある。が、全体を通じて徐々に損失は減っていく。
ミニバッチ確率的勾配降下法 (Minibatch SGD / MSGD)
ミニバッチ確率的勾配降下法では ある程度の固まりから勾配をもとめるので、確率的勾配降下法と比べると 更新処理ごとの損失のばらつきは小さい。
多項ロジスティック回帰 (多クラス)
続けて、目的変数のクラス数が 2より大きい場合 = 多項ロジスティック回帰を行う。データは同じく iris を使う。説明変数は 2次元のままとする。
data = iris
x = data[columns]
y = data['Species']
線形分離できないクラスもあるが、まあやってみよう。
多項ロジスティック回帰の場合は の次元が変わる。伴って、予測値を求める式は
となり関連する変数の次元も以下のように変わる。
: ソフトマックス関数。
の計算結果が各クラスに所属する確率を計算する。詳細は上k(略)
上記のとおり y を多項ロジスティック回帰の形式に変換する。
# ラベルを カテゴリに対応した 0, 1 からなる3列に変換 y = pd.DataFrame({'setosa': (y == 'setosa').astype(int), 'versicolor': (y == 'versicolor').astype(int), 'virginica': (y == 'virginica').astype(int)}) y.head() # setosa versicolor virginica # 1 1 0 0 # 2 1 0 0 # 3 1 0 0 # 4 1 0 0 # 5 1 0 0
続けて、多項ロジスティック回帰で使う関数を定義する。
def p_y_given_x(x, w, b): def softmax(x): return (np.exp(x).T / np.sum(np.exp(x), axis=1)).T return softmax(np.dot(x, w.T) + b) def grad(x, y, w, b): error = p_y_given_x(x, w, b) - y w_grad = np.zeros_like(w) b_grad = np.zeros_like(b) for j in range(w.shape[0]): w_grad[j] = (error.iloc[:, j] * x.T).mean(axis=1) b_grad[j] = error.iloc[:, j].mean() return w_grad, b_grad, np.mean(np.abs(error), axis=0)
ミニバッチ確率的勾配降下法 (Minibatch SGD / MSGD)
勾配法については考え方は一緒なので、MSGD のみ。ほかの手法の場合は ループ箇所とデータ切り出し処理を "2クラスのロジスティック回帰" の場合と同様に書き換えればよい。
def msgd2(x, y, w, b, eta=0.1, num=800, batch_size=10): for i in range(num): for index in range(0, x.shape[0], batch_size): # index 番目のバッチを切り出して処理 _x = x[index:index + batch_size] _y = y[index:index + batch_size] w_grad, b_grad, error = grad(_x, _y, w, b) w -= eta * w_grad b -= eta * b_grad e = np.sum(np.mean(np.abs(y - p_y_given_x(x, w, b)))) yield i, w, b, e
可視化
可視化部分も似たようなコードなので、全体は gist で。各クラスを分割する位置に少しずつ回帰直線が寄っていく。
補足 行列 w の出力がずれていて済まぬ、、済まぬ、、。 tex 記法がうまく動かなかったんだ、、。
まとめ
ロジスティック回帰、とくに勾配法についてまとめた。目的変数のクラス数によってデータの次元 / 予測値算出のロジックは若干異なる。が、おおまかな流れは一緒。
プログラム
これまでのシリーズ
- KalmanFilter の動きを可視化する 一次元版 - StatsFragments
- KalmanFilter の動きを可視化する 二次元版 - StatsFragments
- 動的時間伸縮法 / DTW (Dynamic Time Warping) を可視化する - StatsFragments
11/28追記/修正: 不要な変数定義とその説明を削除。一部グラフのタイトルにスペルミスがあったため修正。
pandas でメモリに乗らない 大容量ファイルを上手に扱う
概要
分析のためにデータ集めしていると、たまに マジか！? と思うサイズの CSV に出くわすことがある。なぜこんなに育つまで放っておいたのか、、、？ このエントリでは普通には開けないサイズの CSV を pandas を使ってうまいこと処理する方法をまとめたい。
サンプルデータ
たまには実データ使おう、ということで WorldBankから GDPデータを落とす。以下のページ右上の "DOWNLOAD DATA" ボタンで CSV を選択し、ローカルに zip を保存する。解凍した "ny.gdp.mktp.cd_Indicator_en_csv_v2.csv" ファイルをサンプルとして使う。
http://data.worldbank.org/indicator/NY.GDP.MKTP.CD?page=1
補足 pandas の Remote Data Access で WorldBank のデータは直接 落っことせるが、今回は ローカルに保存した csv を読み取りたいという設定で。
chunksize を使って ファイルを分割して読み込む
まず、pandas で普通に CSV を読む場合は以下のように pd.read_csv を使う。サンプルデータではヘッダが 3行目から始まるため、冒頭の 2行を skiprows オプションで読み飛ばす。
import pandas as pd fname = 'ny.gdp.mktp.cd_Indicator_en_csv_v2.csv' df = pd.read_csv(fname, skiprows=[0, 1]) df.shape # (258, 59) df.head() # Country Name Country Code Indicator Name Indicator Code 1961 # 0 Aruba ABW GDP (current US$) NY.GDP.MKTP.CD NaN # 1 Andorra AND GDP (current US$) NY.GDP.MKTP.CD NaN # 2 Afghanistan AFG GDP (current US$) NY.GDP.MKTP.CD 5.488889e+08 # 3 Angola AGO GDP (current US$) NY.GDP.MKTP.CD NaN # 4 Albania ALB GDP (current US$) NY.GDP.MKTP.CD NaN # # [5 rows x 59 columns]
では読み取り対象ファイルが メモリに乗らないほど大きい場合はどうするか？ read_csv に chunksize オプションを指定することでファイルの中身を 指定した行数で分割して読み込むことができる。chunksize には 1回で読み取りたい行数を指定する。例えば 50 行ずつ読み取るなら、chunksize=50。
reader = pd.read_csv(fname, skiprows=[0, 1], chunksize=50)
chunksize を指定したとき、返り値は DataFrame ではなく TextFileReader インスタンスとなる。
type(reader) # pandas.io.parsers.TextFileReader
TextFileReader を for でループさせると、ファイルの中身を指定した行数ごとに DataFrame として読み取る。TextFileReader は 現時点での読み取り位置 (ポインタ) を覚えており、ファイルの中身をすべて読み取るとループが終了する。
for r in reader: print(type(r), r.shape) # (<class 'pandas.core.frame.DataFrame'>, (50, 59)) # (<class 'pandas.core.frame.DataFrame'>, (50, 59)) # (<class 'pandas.core.frame.DataFrame'>, (50, 59)) # (<class 'pandas.core.frame.DataFrame'>, (50, 59)) # (<class 'pandas.core.frame.DataFrame'>, (50, 59)) # (<class 'pandas.core.frame.DataFrame'>, (8, 59))
もしくは、TextFileReader.get_chunk で、現在の読み取り位置 (ポインタ) から N 行を読み取る。
# 上のループで全ての要素は読み取り済みになっているので、再度 reader を初期化 reader = pd.read_csv(fname, skiprows=[0, 1], chunksize=50) # 先頭から 5行を読み込み reader.get_chunk(5) # Country Name Country Code Indicator Name Indicator Code 1961 # 0 Aruba ABW GDP (current US$) NY.GDP.MKTP.CD NaN # 1 Andorra AND GDP (current US$) NY.GDP.MKTP.CD NaN # 2 Afghanistan AFG GDP (current US$) NY.GDP.MKTP.CD 5.488889e+08 # 3 Angola AGO GDP (current US$) NY.GDP.MKTP.CD NaN # 4 Albania ALB GDP (current US$) NY.GDP.MKTP.CD NaN # # [5 rows x 59 columns] # 次の4行 (6行目 - 9行目)を読み込み reader.get_chunk(4) # Country Name Country Code Indicator Name Indicator Code 1961 # 0 Andean Region ANR GDP (current US$) NY.GDP.MKTP.CD NaN # 1 Arab World ARB GDP (current US$) NY.GDP.MKTP.CD NaN # 2 United Arab Emirates ARE GDP (current US$) NY.GDP.MKTP.CD NaN # 3 Argentina ARG GDP (current US$) NY.GDP.MKTP.CD NaN # # [4 rows x 59 columns]
補足 chunksize オプションは、 read_csv と read_table で利用できる 。
補足 pandas の read_csv はかなり速いので、パフォーマンス面でも pandas を使うのはよいと思う。
A new high performance, memory-efficient file parser engine for pandas | Wes McKinney
python - Fastest way to parse large CSV files in Pandas - Stack Overflow
でも、分割された DataFrame を扱うのはめんどうなんだけど？
分析に必要なレコードは 全データのごく一部、なんてこともよくある。chunk させて読み込んだ 各グループに対して前処理をかけると、サイズが劇的に減ってメモリに乗る。そんなときは、前処理をかけた後、各 DataFrame を結合してひとつにしたい。
DataFrame の結合には pd.concat。TextFileReader から読み込んだ 各 DataFrame はそれぞれが 0 から始まる index を持つ。pd.concat は既定では index を維持して DataFrame を結合するため、そのままでは index が重複してしまう。 重複を避けるためには、結合後に index を振りなおす = ignore_index=True を指定する必要がある。
よく使う表現は、
def preprocess(x): # 不要な行, 列のフィルタなど、データサイズを削減する処理 return x reader = pd.read_csv(fname, skiprows=[0, 1], chunksize=50) df = pd.concat((preprocess(r) for r in reader), ignore_index=True) df.shape # (258, 59) df.head() # Country Name Country Code Indicator Name Indicator Code 1961 # 0 Aruba ABW GDP (current US$) NY.GDP.MKTP.CD NaN # 1 Andorra AND GDP (current US$) NY.GDP.MKTP.CD NaN # 2 Afghanistan AFG GDP (current US$) NY.GDP.MKTP.CD 5.488889e+08 # 3 Angola AGO GDP (current US$) NY.GDP.MKTP.CD NaN # 4 Albania ALB GDP (current US$) NY.GDP.MKTP.CD NaN # # [5 rows x 59 columns]
メモリ上の DataFrame のサイズを確認したい
DataFrame を適切な大きさに分割するために、DataFrame のメモリ上の占有サイズが知りたい。これは v0.15.0 以降であれば可能。
DataFrame 全体のメモリ上のサイズを表示するには DataFrame.info()。 表示された情報の最終行、 "memory size" をみる。ただし、この表示には object 型のカラムが占めているメモリは含まれないので注意。"+" がついているのは実際の占有領域が表示よりも大きいことを示す。
df.info() # <class 'pandas.core.frame.DataFrame'> # Int64Index: 258 entries, 0 to 257 # Data columns (total 59 columns): # Country Name 258 non-null object # Country Code 258 non-null object # Indicator Name 258 non-null object # Indicator Code 258 non-null object # 1961 127 non-null float64 # ... # 2012 222 non-null float64 # 2013 216 non-null float64 # 2014 0 non-null float64 # Unnamed: 58 0 non-null float64 # dtypes: float64(55), object(4) # memory usage: 116.9+ KB
ちゃんと表示する場合は DataFrame.memory_usage(index=True)。index=True は index の占有メモリも表示に含めるオプション。表示はカラム別になるため、全体のサイズを見たい場合は sum をとる。
df.memory_usage(index=True) # Index 2064 # Country Name 1032 # Country Code 1032 # Indicator Name 1032 # Indicator Code 1032 # 1961 2064 # ... # 2013 2064 # 2014 2064 # Unnamed: 58 2064 # Length: 60, dtype: int64 df.memory_usage(index=True).sum() # 119712
DataFrame を複製せずに処理する
DataFrame のメソッドを呼び出した場合、一般には DataFrame をコピーしてから処理を行って返す (非破壊的な処理を行う)。大容量の DataFrame ではこのコピーによって 実メモリのサイズを超え MemoryError になる / もしくはコピーに非常に時間がかかることがある。
例えば、欠損値の補完を fillna で行う処理をみてみる。fillna の返り値では NaN が 0 でパディングされているが、元データは変更されていない。fillna では、コピーした DataFrame に対して処理を行っていることが "1961" カラムの値を比較するとわかる。
df.fillna(0).head() # Country Name Country Code Indicator Name Indicator Code 1961 # 0 Aruba ABW GDP (current US$) NY.GDP.MKTP.CD 0.000000e+00 # 1 Andorra AND GDP (current US$) NY.GDP.MKTP.CD 0.000000e+00 # 2 Afghanistan AFG GDP (current US$) NY.GDP.MKTP.CD 5.488889e+08 # 3 Angola AGO GDP (current US$) NY.GDP.MKTP.CD 0.000000e+00 # 4 Albania ALB GDP (current US$) NY.GDP.MKTP.CD 0.000000e+00 # dfはそのまま df.head() # Country Name Country Code Indicator Name Indicator Code 1961 # 0 Aruba ABW GDP (current US$) NY.GDP.MKTP.CD NaN # 1 Andorra AND GDP (current US$) NY.GDP.MKTP.CD NaN # 2 Afghanistan AFG GDP (current US$) NY.GDP.MKTP.CD 5.488889e+08 # 3 Angola AGO GDP (current US$) NY.GDP.MKTP.CD NaN # 4 Albania ALB GDP (current US$) NY.GDP.MKTP.CD NaN # # [5 rows x 59 columns]
このコピー処理を避けたい場合は inplace=True を指定する。こいつが指定されると、各メソッドは DataFrame のコピーを行わず、呼び出し元の DataFrame を直接処理する (破壊的な処理を行う)。また、inplace=True が指定されたメソッドは返り値をもたない ( None を返す ) ため、代入を行ってはならない。
# 返り値はないため、代入してはダメ df.fillna(0, inplace=True) # df 自体が変更された df.head() # Country Name Country Code Indicator Name Indicator Code 1961 # 0 Aruba ABW GDP (current US$) NY.GDP.MKTP.CD 0.000000e+00 # 1 Andorra AND GDP (current US$) NY.GDP.MKTP.CD 0.000000e+00 # 2 Afghanistan AFG GDP (current US$) NY.GDP.MKTP.CD 5.488889e+08 # 3 Angola AGO GDP (current US$) NY.GDP.MKTP.CD 0.000000e+00 # 4 Albania ALB GDP (current US$) NY.GDP.MKTP.CD 0.000000e+00 # # [5 rows x 59 columns]
やりたいメソッドが inplace オプションを持つかどうかは API ガイド で確認。
eval による処理の最適化
numexpr のインストール
pd.eval を使うには numexpr パッケージが必要なので、入っていなければインストールする。
pip install numexpr
eval のメリット
pd.eval では文字列表現された式を受け取って評価できる。pd.eval を使うと以下 2つのメリットがある。
- 大容量の
DataFrameを効率的に処理できる - 数値演算を一括して処理できる
補足 numexpr のソースは読んだことがないので詳細不明だが、 pandas では連続するオペレータは逐次処理されていく。例えば "A + B + C" であれば まず "A + B" の結果から DataFrame を作って、さらに "+ C" して DataFrame を作る。DataFrame 生成にはそれなりにコストがかかるので この辺の内部処理が最適化されるのかな、と思っている。
補足 逆に、小さい処理では eval から numexpr を呼び出すコストの方が メリットよりも大きくなるため、pd.eval を使う場合は通常処理とのパフォーマンス比較を行ってからにすべき。今回のサンプルくらいのデータサイズであれば eval よりも 通常の演算のほうが速い。
eval でサポートされる式表現は、公式ドキュメントにあるとおり、
- 数値演算オペレータ。ただし、シフト演算
<<と>>を除く。 - 比較演算オペレータ。連結も可能 (
2 < df < df2とか) - 論理演算オペレータ (
notも含む)。 listとtupleのリテラル表現 ([1, 2]や(1, 2))- プロパティアクセス (
df.a) __getitem__でのアクセス (df[0])
eval の利用
例えばサンプルデータの "2011" 〜 "2013" カラムの値を足し合わせたいときはこう書ける。
pd.eval("df['2011'] + df['2012'] + df['2013']") # 0 2.584464e+09 # 1 0.000000e+00 # 2 5.910162e+10 # 3 3.411516e+11 # 4 3.813926e+10 # ... # 256 6.218183e+10 # 257 3.623061e+10 # Length: 258, dtype: float64
また、pd.eval ではなく DataFrame.eval として呼んだ場合、式表現中はデータのカラム名を変数名としてもつ名前空間で評価される ( DataFrame.query と同じ )。DataFrame.query については以下の記事参照。
まとめ
pandas で大容量ファイルを扱うときは、
chunksizeで分割して読み込み- メソッド呼び出しは
inplace=True - オペレータを使う処理は
pd.eval
これさえ覚えておけば、大容量ファイルが送りつけられた場合でも安心。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (9件) を見る
Python pandas データ選択処理をちょっと詳しく <後編>
概要
こちらの続き。これで pandas でのデータ選択についてはひとまず終わり。
Python pandas データ選択処理をちょっと詳しく <前編> - StatsFragments
Python pandas データ選択処理をちょっと詳しく <中編> - StatsFragments
サンプルデータの準備
データは 前編と同じものを使う。ただし変数名は変えた。
import pandas as pd s1 = pd.Series([1, 2, 3], index = ['I1', 'I2', 'I3']) df1 = pd.DataFrame({'C1': [11, 21, 31], 'C2': [12, 22, 32], 'C3': [13, 23, 33]}, index = ['I1', 'I2', 'I3']) s1 # I1 1 # I2 2 # I3 3 # dtype: int64 df1 # C1 C2 C3 # I1 11 12 13 # I2 21 22 23 # I3 31 32 33
where でのデータ選択
これまでみたとおり、Series から __getitem__ すると、条件に該当する要素のみが返ってくる。
s1[s1 > 2] # I3 3 # dtype: int64
が、ときには元データと同じ長さのデータがほしい場合がある。Series.where を使うと 条件に該当しない ラベルは NaN でパディングし、元データと同じ長さの結果を返してくれる。
where の引数はデータと同じ長さの bool
型の numpy.array もしくは Series である必要がある。まあ Series への論理演算では長さは変わらないので、以下のような使い方なら特別 意識する必要はないかも。
s1.where(s1 > 2) # I1 NaN # I2 NaN # I3 3 # dtype: float64
NaN 以外でパディングしたいんだよ、ってときは 第二引数に パディングに使う値を渡す。NaN でなく 0 でパディングしたければ、
s1.where(s1 > 2, 0) # I1 0 # I2 0 # I3 3 # dtype: int64
また、第二引数にはデータと同じ長さの numpy.array や Series も渡せる。このとき、パディングはそれぞれ対応する位置にある値で行われ、第一引数の条件に該当しないデータを 第二引数で置換するような動きになる。つまり if - else のような表現だと考えていい。
# 第一引数の条件に該当しない s1 の 1, 2番目の要素が array の 1, 2 番目の要素で置換される s1.where(s1 > 2, np.array([4, 5, 6])) # I1 4 # I2 5 # I3 3 # dtype: int64 # 置換用の Series を作る s2 = pd.Series([4, 5, 6], index = ['I1', 'I2', 'I3']) s2 # I1 4 # I2 5 # I3 6 # dtype: int64 # 第一引数の条件に該当しない s1 の 1, 2番目の要素が Series s2 の 1, 2 番目の要素で置換される s1.where(s1 > 2, s2) # I1 4 # I2 5 # I3 3 # dtype: int64
DataFrame でも同様。
df1.where(df1 > 22) # C1 C2 C3 # I1 NaN NaN NaN # I2 NaN NaN 23 # I3 31 32 33 # 0 でパディング df1.where(df1 > 22, 0) # C1 C2 C3 # I1 0 0 0 # I2 0 0 23 # I3 31 32 33 # 置換用の DataFrame を作る df2 = pd.DataFrame({'C1': [44, 54, 64], 'C2': [45, 55, 65], 'C3': [46, 56, 66]}, index = ['I1', 'I2', 'I3']) df2 # C1 C2 C3 # I1 44 45 46 # I2 54 55 56 # I3 64 65 66 # df1 のうち、22以下の値を df2 の値で置換 df1.where(df1 > 22, df2) # C1 C2 C3 # I1 44 45 46 # I2 54 55 23 # I3 31 32 33
where がことさら便利なのは以下のようなケース。
DataFrameに新しいカラムを作りたい。 (あるいは既存の列の値を置き換えたい)- 新しいカラムは、"C2" 列の値が 30を超える場合は "C1" 列の値を使う。
- それ以外は "C3" 列の値を使う。
これが一行でかける。
df1['C4'] = df1['C1'].where(df1['C2'] > 30, df1['C3']) df1 # C1 C2 C3 C4 # I1 11 12 13 13 # I2 21 22 23 23 # I3 31 32 33 31
mask でのデータ選択
where の逆の操作として mask がある。こちらは第一引数の条件に該当するセルを NaN でマスクする。ただし第二引数の指定はできないので 使いどころは限られる。
2015/05/06追記 v0.16.1以降で、DataFrame.mask は where と同じ引数を処理できるようになった。
df1.mask(df1 > 22) # C1 C2 C3 C4 # I1 11 12 13 13 # I2 21 22 NaN NaN # I3 NaN NaN NaN NaN # v0.16.1以降 df1.mask(df1 > 22, 0) # C1 C2 C3 # I1 11 12 13 # I2 21 22 0 # I3 0 0 0
query でのデータ選択
最後にquery。query で何か新しい処理ができるようになるわけではないが、__getitem__ と同じ操作がよりシンプル な表現で書ける。
numexpr のインストール
query を使うには numexpr パッケージが必要なので、入っていなければインストールする。
pip install numexpr
サンプルデータの準備
また、さきほどの例で列追加したのでサンプルデータを作り直す。
df1 = pd.DataFrame({'C1': [11, 21, 31],
'C2': [12, 22, 32],
'C3': [13, 23, 33]},
index = ['I1', 'I2', 'I3'])
df1
# C1 C2 C3
# I1 11 12 13
# I2 21 22 23
# I3 31 32 33
query の利用
__getitem__ を利用したデータ選択では、論理演算の組み合わせで bool の Series を作ってやる必要がある。そのため、[] 内で元データへの参照 ( 下の例では df1 )が繰り返しでてくるし、複数条件の組み合わせの際は 演算順序の都合上 () がでてきたりと式が複雑になりがち。これはわかりにくい。
df1[df1['C1'] > 20] # C1 C2 C3 # I2 21 22 23 # I3 31 32 33 df1[df1['C2'] < 30] # C1 C2 C3 # I1 11 12 13 # I2 21 22 23 df1[(df1['C1'] > 20) & (df1['C2'] < 30)] # C1 C2 C3 # I2 21 22 23
同じ処理は query を使うとすっきり書ける。query の引数にはデータ選択に使う条件式を文字列で渡す。この式が評価される名前空間 = query 名前空間の中では、query を呼び出したデータの列名が あたかも変数のように参照できる。
df1.query('C1 > 20 & C2 < 30') # C1 C2 C3 # I2 21 22 23
ただし、query 名前空間の中で使える表現は限られるので注意。例えば以下のような メソッド呼び出しはできない。
# NG! df1.query('C1.isin([11, 21])') # NotImplementedError: 'Call' nodes are not implemented
同じ処理を行う場合は in 演算子を使う。
# in を使えば OK df1.query('C1 in [11, 21]') # C1 C2 C3 # I1 11 12 13 # I2 21 22 23
また、numexpr で利用できる関数 の呼び出しは query 名前空間上では使えないっぽい。
df1.query('C1 > sqrt(400)') # NotImplementedError: 'Call' nodes are not implemented
そのため、query に渡す式表現では論理表現 以外を含めないほうがよい。条件によって式表現を変えたい、なんて場合は 式表現を都度 文字列として連結 + 生成するか、ローカル変数に計算結果を入れて式表現に渡す。
ローカル変数を式表現中に含める際は、変数名を @ ではじめる。
x = 20 # NG! df1.query('C1 > x') # UndefinedVariableError: name 'x' is not defined # OK! df1.query('C1 > @x') # C1 C2 C3 # I2 21 22 23 # I3 31 32 33
query 名前空間上で index の値を参照する場合は、式表現中で index と指定する。
df1.query('index in ["I1", "I2"]') # C1 C2 C3 # I1 11 12 13 # I2 21 22 23
index が列名と重複した場合は 列名が優先。
df_idx = pd.DataFrame({'index': [1, 2, 3]}, index=[3, 2, 1])
df_idx
# index
# 3 1
# 2 2
# 1 3
df_idx.query('index >= 2')
# index
# 2 2
# 1 3
まとめ
whereを使うと 元データと同じ形のデータが取得できる。また、if - else的表現が一行で書ける。- 複数の条件式を組み合わせる場合は
queryを使うとシンプルに書ける。
全三回で pandas でのデータ選択処理をまとめた。公式ドキュメント、またAPIガイドには ここで記載しなかったメソッド / 例もあるので、より尖った使い方をしたい人は読んでおくといいと思う。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (9件) を見る
Python pandas データ選択処理をちょっと詳しく <中編>
こちらの続き。
上の記事では bool でのデータ選択について 最後にしれっと書いて終わらせたのだが、一番よく使うところなので中編として補足。
まず __getitem__ や ix の記法では、次のような指定によって 行 / 列を選択することができた。
index,columnsのラベルを直接指定しての選択index,columnsの番号(順序)を指定しての選択index,columnsに対応するboolのリストを指定しての選択
ここでは上記の選択方法をベースとして、ユースケースごとに Index や Series のプロパティ / メソッドを使ってできるだけシンプルにデータ選択を行う方法をまとめる。
補足 一部の内容はこちらの記事ともかぶる。下の記事のほうが簡単な内容なので、必要な方はまずこちらを参照。
簡単なデータ操作を Python pandas で行う - StatsFragments
準備
今回のサンプルは DataFrame だけで。
import pandas as pd import numpy as np df = pd.DataFrame({'N1': [1, 2, 3, 4, 5, 6], 'N2': [10, 20, 30, 40, 50, 60], 'N3': [6, 5, 4, 3, 2, 1], 'F1': [1.1, 2.2, 3.3, 4.4, 5.5, 6.6], 'F2': [1.1, 2.2, 3.3, 4.4, 5.5, 6.6], 'S1': ['A', 'b', 'C', 'D', 'E', 'F'], 'S2': ['A', 'X', 'X', 'X', 'E', 'F'], 'D1': pd.date_range('2014-11-01', freq='D', periods=6)}, index=pd.date_range('2014-11-01', freq='M', periods=6), columns=['N1', 'N2', 'N3', 'F1', 'F2', 'S1', 'S2', 'D1']) df # N1 N2 N3 F1 F2 S1 S2 D1 # 2014-11-30 1 10 6 1.1 1.1 A A 2014-11-01 # 2014-12-31 2 20 5 2.2 2.2 b X 2014-11-02 # 2015-01-31 3 30 4 3.3 3.3 C X 2014-11-03 # 2015-02-28 4 40 3 4.4 4.4 D X 2014-11-04 # 2015-03-31 5 50 2 5.5 5.5 E E 2014-11-05 # 2015-04-30 6 60 1 6.6 6.6 F F 2014-11-06
補足 後の例のために index は月次の日付型にした。pd.date_range の使い方は以下の記事参照。
Python pandas で日時関連のデータ操作をカンタンに - StatsFragments
以下の例では __getitem__ を使ってシンプルに書けるところは __getitem__ で書く。もちろん ix, (ラベルや番号なら loc, iloc) を使っても書ける。
index, columns のラベルを特定の条件で選択
前編で index, columns のラベルそのものを直接指定すればデータ選択できるのはわかった。が、ラベルが特定の条件を満たすとき (ラベルが特定の文字で始まるとか) だけ選択するにはどうすれば？というのがここでの話。
補足 内部実装の話になるが、 index, columns は どちらも pd.Index 型のクラスが使われている ( DatetimeIndex は Index のサブクラス)。index, columns とも裏側にあるオブジェクトは同一のため、このセクションで記載する方法は 行 / 列が入れ替わっても使える。
df.index # <class 'pandas.tseries.index.DatetimeIndex'> # [2014-11-30, ..., 2015-04-30] # Length: 6, Freq: M, Timezone: None df.columns # Index([u'N1', u'N2', u'N3', u'F1', u'F2', u'S1', u'S2', u'D1'], dtype='object')
ラベルに関数適用して選択したい
pd.Index.map。たとえば 大文字の "N" から始まるラベル名のみ抽出したいなら
df.columns.map(lambda x: x.startswith('N')) # array([ True, True, True, False, False, False, False, False], dtype=bool) df.ix[:, df.columns.map(lambda x: x.startswith('N'))] # N1 N2 N3 # 2014-11-30 1 10 6 # 2014-12-31 2 20 5 # 2015-01-31 3 30 4 # 2015-02-28 4 40 3 # 2015-03-31 5 50 2 # 2015-04-30 6 60 1
ということで map を使えば index, columns のラベルに対してあらゆる関数を適用してデータ選択できる。
さらに、よく使うと思われるケースではより簡便な方法が用意されている。
2015/05/06追記 v0.16.1 で Index にも .str アクセサが追加され、文字列処理関数を直接呼び出せるようになった。そのため、上の例は df.ix[:, df.columns.str.startswith('N')] とも書ける。.str アクセサについては以降の記載を参照。
リストに含まれるラベルだけ選択したい
たとえば選択したいラベルのリストがあり、そこに含まれるものだけ選択したいなんてことがある。選択候補リストに余計なラベルが含まれていると、__getitem__ では KeyError になり、ix では NaN (値のない) 列ができてしまう。
df[['N1', 'N2', 'N4']] # KeyError: "['N4'] not in index" df.ix[:, ['N1', 'N2', 'N4']] # N1 N2 N4 # 2014-11-30 1 10 NaN # 2014-12-31 2 20 NaN # 2015-01-31 3 30 NaN # 2015-02-28 4 40 NaN # 2015-03-31 5 50 NaN # 2015-04-30 6 60 NaN
いや NaN とかいいから 存在する列だけが欲しいんだけど、、、というときに pd.Index.isin。
df.columns.isin(['N1', 'N2', 'N4']) # array([ True, True, False, False, False, False, False, False], dtype=bool) df.ix[:, df.columns.isin(['N1', 'N2', 'N4'])] # N1 N2 # 2014-11-30 1 10 # 2014-12-31 2 20 # 2015-01-31 3 30 # 2015-02-28 4 40 # 2015-03-31 5 50 # 2015-04-30 6 60
ラベルをソートして選択したい
pd.Index.order。columns をアルファベット順に並べ替えて、前から3つを取得したければ、
df.columns.order() # Index([u'D1', u'F1', u'F2', u'N1', u'N2', u'N3', u'S1', u'S2'], dtype='object') df.columns.order()[:3] # Index([u'D1', u'F1', u'F2'], dtype='object') df[df.columns.order()[:3]] # D1 F1 F2 # 2014-11-30 2014-11-01 1.1 1.1 # 2014-12-31 2014-11-02 2.2 2.2 # 2015-01-31 2014-11-03 3.3 3.3 # 2015-02-28 2014-11-04 4.4 4.4 # 2015-03-31 2014-11-05 5.5 5.5 # 2015-04-30 2014-11-06 6.6 6.6
特定の年, 月, etc... のデータだけ選択したい
DatetimeIndex へのプロパティアクセスを使う。使えるプロパティはこちら。
index が 2015年の日付になっている行のみ抽出するときは、pd.DatetimeIndex.year で 年のみを含む numpy.array を作って論理演算する。
df.index.year # array([2014, 2014, 2015, 2015, 2015, 2015], dtype=int32) df.index.year == 2015 # array([False, False, True, True, True, True], dtype=bool) df[df.index.year == 2015] # N1 N2 N3 F1 F2 S1 S2 D1 # 2015-01-31 3 30 4 3.3 3.3 C X 2014-11-03 # 2015-02-28 4 40 3 4.4 4.4 D X 2014-11-04 # 2015-03-31 5 50 2 5.5 5.5 E E 2014-11-05 # 2015-04-30 6 60 1 6.6 6.6 F F 2014-11-06
補足 前編でまとめた原則には書かなかったが、DatetimeIndex (と PeriodIndex という日時関連の別クラス ) を index に持つ Series, DataFrame では、 例外的に __getitem__ の引数として日時-like な文字列が使えたりもする。詳しくは こちら。そのため、同じ処理は以下のようにも書ける。
df['2015'] # N1 N2 N3 F1 F2 S1 S2 D1 # 2015-01-31 3 30 4 3.3 3.3 C X 2014-11-03 # 2015-02-28 4 40 3 4.4 4.4 D X 2014-11-04 # 2015-03-31 5 50 2 5.5 5.5 E E 2014-11-05 # 2015-04-30 6 60 1 6.6 6.6 F F 2014-11-06
ラベルが重複したデータを削除したい
2015/05/06修正 v0.16向けの内容に変更
これは v0.16 以降であれば簡単。例示のため index に "A" が3つ重複したデータを作る。
df_dup = pd.DataFrame({'N1': [1, 2, 3, 4],
'N2': [6, 5, 4, 3],
'S1': ['A', 'B', 'C', 'D']},
index=['A', 'A', 'A', 'B'])
df_dup
# N1 N2 S1
# A 1 6 A
# A 2 5 B
# A 3 4 C
# B 4 3 D
Index.duplicated() で Index の値が重複しているかどうかがわかるため、これを使って列選択すればよい。重複している値が True となっているため ~ で論理否定をとって行選択すると、
df_dup.index.duplicated() # array([False, True, True, False], dtype=bool) df_dup[~df_dup.index.duplicated()] # N1 N2 S1 # A 1 6 A # B 4 3 D
このとき、各重複グループの一番最初のデータは削除されない。この挙動は take_last オプションで変更できる。
take_last=False (Default): 重複しているグループの最初のデータを残すtake_last=True: 重複しているグループの最後のデータを残す。
df_dup[~df_dup.index.duplicated(take_last=True)] # N1 N2 S1 # A 3 4 C # B 4 3 D
いやいや重複データは全削除したいんですけど？という場合は Dummy 列で groupby して filter。
df_dup.groupby(["Dummy"]).filter(lambda x:x.shape[0] == 1) # N1 N2 S1 Dummy # B 4 3 D B df_dup.groupby(["Dummy"]).filter(lambda x:x.shape[0] == 1)[['N1', 'N2', 'S1']] # N1 N2 S1 # B 4 3 D
補足 直感的でないな、と思った方が多いと思うが 開発者でもタスクとしては認識している ので、、、。
列, 行の値から特定の条件で選択
上のセクションでは index, columns 自体のラベルからデータ選択する方法を書いた。ここからは データの中身 (行 / 列 / もしくは各セルの値 ) をもとにデータ選択する方法を記載する。
補足 DataFrame では 実データは列持ち (各列が特定の型のデータを保持している) なので、ここからの方法では行/列の方向を意識する必要がある。
特定の型の列のみ取り出す
DataFrame.dtypes プロパティで各カラムの型が取得できるので、それらに対して論理演算をかける。
df.dtypes # N1 int64 # N2 int64 # N3 int64 # F1 float64 # F2 float64 # S1 object # S2 object # D1 datetime64[ns] # dtype: object df.dtypes == np.float64 # N1 False # N2 False # N3 False # F1 True # F2 True # S1 False # S2 False # D1 False # dtype: bool df.ix[:, df.dtypes == np.float64] # F1 F2 # 2014-11-30 1.1 1.1 # 2014-12-31 2.2 2.2 # 2015-01-31 3.3 3.3 # 2015-02-28 4.4 4.4 # 2015-03-31 5.5 5.5 # 2015-04-30 6.6 6.6
値が特定の条件を満たす行/列を選択したい
たいていは 普通の演算でいける。"N1" カラムの値が偶数の行だけ抽出するには、
df['N1'] % 2 == 0 # 2014-11-30 False # 2014-12-31 True # 2015-01-31 False # 2015-02-28 True # 2015-03-31 False # 2015-04-30 True # Freq: M, Name: N1, dtype: bool df[df['N1'] % 2 == 0] # N1 N2 N3 F1 F2 S1 S2 D1 # 2014-12-31 2 20 5 2.2 2.2 b X 2014-11-02 # 2015-02-28 4 40 3 4.4 4.4 D X 2014-11-04 # 2015-04-30 6 60 1 6.6 6.6 F F 2014-11-06
各列の合計値が 50 を超えるカラムを抽出するには、
df.sum() # N1 21.0 # N2 210.0 # N3 21.0 # F1 23.1 # F2 23.1 # dtype: float64 indexer = df.sum() > 50 indexer # N1 False # N2 True # N3 False # F1 False # F2 False # dtype: bool df[indexer.index[indexer]] # N2 # 2014-11-30 10 # 2014-12-31 20 # 2015-01-31 30 # 2015-02-28 40 # 2015-03-31 50 # 2015-04-30 60
各行 の値に関数適用して選択したいときは apply。apply に渡す関数は 行 もしくは 列を Series として受け取って処理できるものでないとダメ。apply での関数の適用方向は axis オプションで決める。
axis=0: 各列への関数適用axis=1: 各行への関数適用
"N1" カラムと "N2" カラムの積が 100 を超える行だけをフィルタする場合、
df.apply(lambda x: x['N1'] * x['N2'], axis=1) # 2014-11-30 10 # 2014-12-31 40 # 2015-01-31 90 # 2015-02-28 160 # 2015-03-31 250 # 2015-04-30 360 # Freq: M, dtype: int64 df.apply(lambda x: x['N1'] * x['N2'], axis=1) > 100 # 2014-11-30 False # 2014-12-31 False # 2015-01-31 False # 2015-02-28 True # 2015-03-31 True # 2015-04-30 True # Freq: M, dtype: bool df[df.apply(lambda x: x['N1'] * x['N2'], axis=1) > 100] # N1 N2 N3 F1 F2 S1 S2 D1 # 2015-02-28 4 40 3 4.4 4.4 D X 2014-11-04 # 2015-03-31 5 50 2 5.5 5.5 E E 2014-11-05 # 2015-04-30 6 60 1 6.6 6.6 F F 2014-11-06
補足 上ではあえて apply を使ったが、各列同士は直接 要素の積をとれるため別に apply が必須ではない。
df[df['N1'] * df['N2'] > 100] # N1 N2 N3 F1 F2 S1 S2 D1 # 2015-02-28 4 40 3 4.4 4.4 D X 2014-11-04 # 2015-03-31 5 50 2 5.5 5.5 E E 2014-11-05 # 2015-04-30 6 60 1 6.6 6.6 F F 2014-11-06
リストに含まれる値だけ選択したい
index の例と同じ。Series も isin メソッドを持っているので、"S1" カラムの値が "A" もしくは "D" の列を選択するときは、
df['S1'].isin(['A', 'D']) # 2014-11-30 True # 2014-12-31 False # 2015-01-31 False # 2015-02-28 True # 2015-03-31 False # 2015-04-30 False # Freq: M, Name: S1, dtype: bool df[df['S1'].isin(['A', 'D'])] # N1 N2 N3 F1 F2 S1 S2 D1 # 2014-11-30 1 10 6 1.1 1.1 A A 2014-11-01 # 2015-02-28 4 40 3 4.4 4.4 D X 2014-11-04
値をソートして選択したい
DataFrame.sort。ソート順序の変更など、オプションの詳細はこちら。
"N2" カラムの値が大きいものを 上から順に 3行分 取得するには、ソートして 行番号でスライスすればよい。
df.sort('N2', ascending=False)[:3] # N1 N2 N3 F1 F2 S1 S2 D1 # 2015-04-30 6 60 1 6.6 6.6 F F 2014-11-06 # 2015-03-31 5 50 2 5.5 5.5 E E 2014-11-05 # 2015-02-28 4 40 3 4.4 4.4 D X 2014-11-04
特定の年, 月, etc... のデータだけ選択したい
日時型のカラムに対しては、dt アクセサを利用して日時型の各プロパティにアクセスできる。使えるプロパティはこちら。
"D1" カラムの日付が 2日, 3日, 5日の行だけ取得したければ、dt アクセサ + isin で、
df['D1'] # 2014-11-30 2014-11-01 # 2014-12-31 2014-11-02 # 2015-01-31 2014-11-03 # 2015-02-28 2014-11-04 # 2015-03-31 2014-11-05 # 2015-04-30 2014-11-06 # Freq: M, Name: D1, dtype: datetime64[ns] df['D1'].dt.day # 2014-11-30 1 # 2014-12-31 2 # 2015-01-31 3 # 2015-02-28 4 # 2015-03-31 5 # 2015-04-30 6 # Freq: M, dtype: int64 df['D1'].dt.day.isin([2, 3, 5]) # 2014-11-30 False # 2014-12-31 True # 2015-01-31 True # 2015-02-28 False # 2015-03-31 True # 2015-04-30 False # Freq: M, dtype: bool df[df['D1'].dt.day.isin([2, 3, 5])] # N1 N2 N3 F1 F2 S1 S2 D1 # 2014-12-31 2 20 5 2.2 2.2 b X 2014-11-02 # 2015-01-31 3 30 4 3.3 3.3 C X 2014-11-03 # 2015-03-31 5 50 2 5.5 5.5 E E 2014-11-05
補足 dt アクセサは集計でもかなり強力なのでオススメ。
Python pandas アクセサ / Grouperで少し高度なグルーピング/集計 - StatsFragments
Python 組み込みの文字列関数を使ってデータ選択したい
object型のカラムに対しては、str アクセサを利用して、Python 組み込みの 文字列関数を実行した結果を Series として取得できる。使えるメソッドはこちら。したがって、文字列関数を実行するだけならわざわざ apply を使わなくても済む。
str.lower を使って値が小文字の列を選択してみる。
# "S1" カラムの値を小文字化 df['S1'].str.lower() # 2014-11-30 a # 2014-12-31 b # 2015-01-31 c # 2015-02-28 d # 2015-03-31 e # 2015-04-30 f # Freq: M, Name: S1, dtype: object df['S1'].str.lower() == df['S1'] # 2014-11-30 False # 2014-12-31 True # 2015-01-31 False # 2015-02-28 False # 2015-03-31 False # 2015-04-30 False # Freq: M, Name: S1, dtype: bool df[df['S1'].str.lower() == df['S1']] # N1 N2 N3 F1 F2 S1 S2 D1 # 2014-12-31 2 20 5 2.2 2.2 b X 2014-11-02
2015/05/06追記 v0.16以降で .str から呼び出せる関数が追加されている。追加関数には str.islower もあるため、上の例は以下のように書ける。
df[df['S1'].str.islower()] # N1 N2 N3 F1 F2 S1 S2 D1 # 2014-12-31 2 20 5 2.2 2.2 b X 2014-11-02
値が重複したデータを削除したい
DataFrame.drop_duplicates。重複のうち、最初のひとつ もしくは 最後のひとつ どちらかを残す挙動は Index のときと同様。
df_dup2 = pd.DataFrame({'N1': [1, 1, 3, 1],
'N2': [1, 1, 4, 4],
'S1': ['A', 'A', 'B', 'C']})
df_dup2
# N1 N2 S1
# 0 1 1 A
# 1 1 1 A
# 2 3 4 B
# 3 1 4 C
# 完全に重複している行を、最初の一つを残して削除
df_dup2.drop_duplicates()
# N1 N2 S1
# 0 1 1 A
# 2 3 4 B
# 3 1 4 C
# N1 カラムの値が重複している行を、最初の一つを残して削除
df_dup2.drop_duplicates(subset=['N1'])
# N1 N2 S1
# 0 1 1 A
# 2 3 4 B
重複している値 全てを削除したければ、Index の場合と同じく groupby して filter。
df_dup2.groupby(["S1"]).filter(lambda x:x.shape[0] == 1) # N1 N2 S1 # 2 3 4 B # 3 1 4 C
欠測値 ( NaN ) のデータを削除したい
DataFrame.dropna。チェックするデータの方向 (行方向 or 列方向), 条件などのオプションの詳細はこちら。
サンプルとして、いくつかの欠測値 ( NaN ) を含むデータを作成。
df_na = pd.DataFrame({'N1': [1, 2, np.nan, 4],
'N2': [6, 5, 4, np.nan],
'S1': [np.nan, 'B', 'C', 'D']})
df_na
# N1 N2 S1
# 0 1 6 NaN
# 1 2 5 B
# 2 NaN 4 C
# 3 4 NaN D
# NaN がひとつでもある行を削除
df_na.dropna()
# N1 N2 S1
# 1 2 5 B
# N2 カラムに NaN がある行を削除
df_na.dropna(subset=['N2'])
# N1 N2 S1
# 0 1 6 NaN
# 1 2 5 B
# 2 NaN 4 C
# 0 or 1 行目に NaN がある列を削除
df_na.dropna(axis=1, subset=[0, 1])
# N1 N2
# 0 1 6
# 1 2 5
# 2 NaN 4
# 3 4 NaN
まとめ
pandas で行 / 列からデータ選択する方法をざっとまとめた。少し複雑なデータ選択をするときにやるべき流れは、
- 対象となる 行 / 列のデータと選択条件を確認する
- API ガイドにやりたい処理に直接 対応するものがないか探す。
- 直接 対応するものがなければ、全体をいくつかの単純な処理に分解して、最終的に
__getitem__やixに渡せる形にできるか考える。 - 自作関数を使ったほうが早そうだな、、、というときには
Index.maporDataFrame.apply。
次こそ where と query。
11/18 追記:
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (9件) を見る
Python pandas データ選択処理をちょっと詳しく <前編>
概要
書いていて長くなったため、まず前編として pandas で データを行 / 列から選択する方法を少し詳しく書く。特に、個人的にはけっこう重要だと思っている loc と iloc について 日本語で整理したものがなさそうなので。
サンプルデータの準備
import pandas as pd s = pd.Series([1, 2, 3], index = ['I1', 'I2', 'I3']) df = pd.DataFrame({'C1': [11, 21, 31], 'C2': [12, 22, 32], 'C3': [13, 23, 33]}, index = ['I1', 'I2', 'I3']) s # I1 1 # I2 2 # I3 3 # dtype: int64 df # C1 C2 C3 # I1 11 12 13 # I2 21 22 23 # I3 31 32 33
__getitem__ での選択
Series, もしくは DataFrame に対して、__getitem__ でアクセスした場合の挙動は以下のようになる。
Series:indexからの選択DataFrame:columnsからの選択 (例外あり、下表 * の箇所)
なんかあいまいな言い方なのは 下表のとおりいくつかの形式で index, columns を指定できるから。使える形式は Series と DataFrame で少し違うので注意。
Series |
DataFrame |
|---|---|
index の名前 (もしくは名前のリスト) |
columns の名前 (もしくは名前のリスト) |
index の順序 (番号) (もしくは番号のリスト) |
columns の順序 (番号)のリスト (番号のみはダメ) ただし、slice を渡した場合は index からの選択 * |
indexと同じ長さの bool のリスト |
indexと同じ長さの bool のリスト * |
| - | データと同じ次元の bool の DataFrame |
※ 上で "リスト" と書いている箇所は numpy.array, pd.Series, slice なんかのリスト - like でもよい。
11/15追記: 元データと同じ次元の DataFrame が渡せることの記載漏れ、またDataFrame へ slice を渡した場合の挙動が間違っていたので修正。
s[0] # 1 s['I1'] # 1 df['C1'] # I1 11 # I2 21 # I3 31 # Name: C1, dtype: int64 # 番号を数値として渡すとNG! df[1] # KeyError: 1 # 番号のリストならOK (columns からの選択) df[[1]] # C2 # I1 12 # I2 22 # I3 32 # 番号のスライスもOK (index からの選択) df[1:2] # C1 C2 C3 # I2 21 22 23 # NG! df['I1'] # KeyError: 'I1' s[[True, False, True]] # I1 1 # I3 3 # dtype: int64 df[[True, False, True]] # C1 C2 C3 # I1 11 12 13 # I3 31 32 33 # bool の DataFrame を作る df > 21 # C1 C2 C3 # I1 False False False # I2 False True True # I3 True True True # bool の DataFrame で選択 df[df > 21] # C1 C2 C3 # I1 NaN NaN NaN # I2 NaN 22 23 # I3 31 32 33
引数による 返り値の違い
引数を値だけ (ラベル, 数値)で渡すと次元が縮約され、Series では単一の値, DataFrame では Series が返ってくる。元のデータと同じ型の返り値がほしい場合は 引数をリスト型にして渡せばいい。
# 返り値は 値 s['I1'] # 1 # 返り値は Series s[['I1']] # I1 1 # dtype: int64 # 返り値は Series df['C1'] # I1 11 # I2 21 # I3 31 # Name: C1, dtype: int64 # 返り値は DataFrame df[['C1']] # C1 # I1 11 # I2 21 # I3 31
index, columns を元にした選択 ( ix, loc, iloc )
ix プロパティを使うと DataFrame で index, columns 両方を指定してデータ選択を行うことができる ( Series の挙動は __getitem__ と同じ)。
引数として使える形式は __getitem__ と同じだが、ix では DataFrame も以下すべての形式を使うことができる。
- 名前 (もしくは名前のリスト)
- 順序 (番号) (もしくは番号のリスト)
index, もしくはcolumnsと同じ長さのboolのリスト
ix はメソッドではなくプロパティなので、呼び出しは以下のようになる。
Series.ix[?]:?にはindexを特定できるものを指定DataFrame.ix[?, ?]:?にはそれぞれindex,columnsの順に特定できるものを指定
# 名前による指定 s.ix['I2'] # 2 df.ix['I2', 'C2'] # 22 # 順序による指定 s.ix[1] # 2 df.ix[1, 1] # 22 # 名前のリストによる指定 s.ix[['I1', 'I3']] # I1 1 # I3 3 # dtype: int64 df.ix[['I1', 'I3'], ['C1', 'C3']] # C1 C3 # I1 11 13 # I3 31 33 # bool のリストによる指定 s.ix[[True, False, True]] # I1 1 # I3 3 # dtype: int64 df.ix[[True, False, True], [True, False, True]] # C1 C3 # I1 11 13 # I3 31 33 # 第一引数, 第二引数で別々の形式を使うこともできる df.ix[1:, "C1"] # I2 21 # I3 31 # Name: C1, dtype: int64
DataFrame.ix の補足
DataFrameで第二引数を省略した場合は index への操作になる。
df.ix[1] # C1 21 # C2 22 # C3 23 # Name: I2, dtype: int64
DataFrame で columns に対して操作したい場合、以下のように第一引数を空にするとエラーになる ( R に慣れているとやりがち、、 )。第一引数には : を渡す必要がある (もしくは ixを使わず 直接 __getitem__ する )。
df.ix[, 'C3'] # SyntaxError: invalid syntax df.ix[:, 'C3'] I1 13 I2 23 I3 33 Name: C3, dtype: int64
引数による 返り値の違い
引数の型による返り値の違いは __getitem__ の動きと同じ。Series については挙動もまったく同じなので、ここでは DataFrame の場合だけ例示。
# 返り値は 値 df.ix[1, 1] # 22 # 返り値は Series df.ix[[1], 1] # I2 22 # Name: C2, dtype: int64 # 返り値は DataFrame df.ix[[1], [1]] # C2 # I2 22
ということで、だいたいの場合は ix を使えばよしなにやってくれる。
とはいえ
ラベル名 もしくは 番号 どちらかだけを指定してデータ選択したい場合もある。例えば、
index,columnsがint型であるindex,columnsに重複がある
df2 = pd.DataFrame({1: [11, 21, 31],
2: [12, 22, 32],
3: [13, 23, 33]},
index = [2, 2, 2])
df2
# 1 2 3
# 2 11 12 13
# 2 21 22 23
# 2 31 32 33
内部的には ix は ラベルを優先して処理を行うため、上記のデータについては以下のような動作をする。
indexに2を指定すると、3行目ではなくindexのラベルが 2 の行を選択columnsに[1, 2]を指定すると、2, 3列目ではなくcolumnsのラベルが 1, 2 の列を選択
df2.ix[2, [1, 2]] 1 2 2 11 12 2 21 22 2 31 32
つまり 上記のようなデータ、(特にデータによって index, columns の値が変わるとか、、) では ix を使うと意図しない挙動をする可能性がある。明示的に index, columns を番号で指定したい！というときには iloc を使う。
df2.iloc[2, [1, 2]] 2 32 3 33 Name: 2, dtype: int64 # 3列目は存在しないので NG! df2.iloc[2, 3] # IndexError: index out of bounds
同様に、明示的にラベルのみで選択したい場合は loc。
df2.loc[2, [1, 2]] 1 2 2 11 12 2 21 22 2 31 32 # ラベルが 1 の index は存在しないので NG! df.loc[1, 2] # KeyError: 'the label [1] is not in the [index]'
ix, iloc, loc については文法 / 挙動は基本的に一緒で、使える引数の形式のみが異なる。整理すると、
| 使える引数の形式 | ix |
iloc |
loc |
|---|---|---|---|
| 名前 (もしくは名前のリスト) | ○ | - | ○ |
| 順序 (番号) (もしくは番号のリスト) | ○ | ○ | - |
index, もしくは columns と同じ長さの bool のリスト |
○ | ○ | ○ |
一般に、pandas でスクリプトを書くような場合には index, columns をどちらの形式で指定して選択するか決まっているはず。そんなときは loc , iloc を使っておいたほうが安全だと思う。
プロパティによるアクセス
Series の index, DataFrame の columns はプロパティとしてもアクセスできる。が、オブジェクト自体のメソッド/プロパティと衝突した場合は メソッド/プロパティが優先されるので使わないほうがよい。例えば 上の ix を列名に持つデータがあったとすると、
s.I1 # 1 df3 = pd.DataFrame({'C1': [11, 21, 31], 'C2': [12, 22, 32], 'ix': [13, 23, 33]}, index = ['I1', 'I2', 'I3']) df3 # C1 C2 ix # I1 11 12 13 # I2 21 22 23 # I3 31 32 33 df3.C1 # I1 11 # I2 21 # I3 31 # Name: C1, dtype: int64 # NG! df3.ix <pandas.core.indexing._IXIndexer at 0x10bb31890>
bool のリストによってデータ選択ができるということは
pandas.Series や numpy.array への演算は原則 リスト内の各要素に対して適用され、結果は真偽値の numpy.array (もしくは Series ) になる。また、真偽値の numpy.array 同士で論理演算することもできる。
df.columns == 'C3' # array([False, False, True], dtype=bool) df.columns.isin(['C1', 'C2']) # array([ True, True, False], dtype=bool) (df.columns == 'C3') | (df.columns == 'C2') # array([False, True, True], dtype=bool)
そのため、上記のような条件式をそのまま行列選択時の引数として使うことができる。
df.ix[df.index.isin(['I1', 'I2']), df.columns == 'C3'] # C3 # I1 13 # I2 23
上の例ではあまりありがたみはないと思うが、外部関数で bool のリストを作ってしまえばどんなに複雑な行列選択でもできるのは便利。
まとめ
pandas で行 / 列からデータ選択するとき、
- 手元で対話的にちょっと試す場合は
ixが便利。 - ある程度の期間使うようなスクリプトを書く場合は 少し面倒でも
iloc,locが安全。
後編はもう少し複雑なデータ選択( where, query ) あたりを予定。
11/15追記: 取消線した内容とは違うが続きを書いた。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (9件) を見る
Python pandas で日時関連のデータ操作をカンタンに
概要
Python で日時/タイムスタンプ関連の操作をする場合は dateutil や arrow を使っている人が多いと思うが、 pandas でもそういった処理がわかりやすく書けるよ、という話。
pandas の本領は多次元データの蓄積/変形/集約処理にあるが、日時操作に関連した強力なメソッド / ユーティリティもいくつか持っている。今回は それらを使って日時操作を簡単に行う方法を書いてく。ということで DataFrame も Series もでてこない pandas 記事のはじまり。
※ ここでいう "日時/タイムスタンプ関連の操作" は文字列パース、日時加算/減算、タイムゾーン設定、条件に合致する日時のリスト生成などを想定。時系列補間/リサンプリングなんかはまた膨大になるので別途。
インストール
以下サンプルには 0.15での追加機能も含まれるため、0.15 以降が必要。
pip install pandas
準備
import pandas as pd
初期化/文字列パース
pd.to_datetime が便利。返り値は Python 標準のdatetime.datetime クラスを継承したpd.Timestamp インスタンスになる。
dt = pd.to_datetime('2014-11-09 10:00') dt # Timestamp('2014-11-09 10:00:00') type(dt) # pandas.tslib.Timestamp import datetime isinstance(dt, datetime.datetime) # True
pd.to_datetime は、まず pandas 独自の日時パース処理を行い、そこでパースできなければ dateutil.parser.parse を呼び出す。そのため結構無茶なフォーマットもパースできる。
pd.to_datetime('141109 1005') # Timestamp('2014-11-09 10:05:00')
リスト-likeなものを渡すと DatetimeIndex ( pandas を 普段 使ってない方は 日時のリストみたいなものだと思ってください) が返ってくる。
pd.to_datetime(['2014-11-09 10:00', '2014-11-10 10:00']) # <class 'pandas.tseries.index.DatetimeIndex'> # [2014-11-09 10:00:00, 2014-11-10 10:00:00] # Length: 2, Freq: None, Timezone: None
また、pd.to_datetime は文字列以外の 日付-like なものも処理できる。自分は numpy.datetime64 を datetime.datetimeに変換するのに使ったりもする。
pd.to_datetime(datetime.datetime(2014, 11, 9)) # Timestamp('2014-11-09 00:00:00') import numpy as np pd.to_datetime(np.datetime64('2014-11-09 10:00Z')) # Timestamp('2014-11-09 10:00:00')
補足 Timestamp.to_datetime() で Timestamp -> datetime.datetime へ変換できる
dt.to_datetime()
# datetime.datetime(2014, 11, 9, 10, 0)
日時の加算/減算
pd.Timestamp に対して datetime.timedelta, dateutil.relativedeltaを使って日時の加減算を行うこともできるが、
dt + datetime.timedelta(days=1) # Timestamp('2014-11-10 10:00:00') from dateutil.relativedelta import relativedelta dt + relativedelta(days=1) # Timestamp('2014-11-10 10:00:00')
一般的な時間単位は pandas が offsets として定義しているため、そちらを使ったほうが直感的だと思う。使えるクラスの一覧は こちら。
dt # Timestamp('2014-11-09 10:00:00') import pandas.tseries.offsets as offsets # 翌日 dt + offsets.Day() # Timestamp('2014-11-10 10:00:00') # 3日後 dt + offsets.Day(3) # Timestamp('2014-11-12 10:00:00')
pd.offsets を使うメリットとして、例えば dateutil.relativedelta では days と day を間違えるとまったく違う意味になってしまう。が、pd.offsets ならより明瞭に書ける。
# relativedelta の場合 days なら 3日後 dt + relativedelta(days=3) # Timestamp('2014-11-12 10:00:00') # day なら月の 3日目 dt + relativedelta(day=3) # Timestamp('2014-11-03 10:00:00') # 月の3日目を取りたいならこっちのが明瞭 dt - offsets.MonthBegin() + offsets.Day(2) # Timestamp('2014-11-03 10:00:00')
時刻部分を丸めたければ normalize=True。
dt + offsets.MonthEnd(normalize=True) # Timestamp('2014-11-30 00:00:00')
また、 pd.offsets は datetime.datetime, np.datetime64 にも適用できる。
datetime.datetime(2014, 11, 9, 10, 00) + offsets.Week(weekday=2) # Timestamp('2014-11-12 10:00:00')
※ np.datetime64 に対して適用する場合は 加算/減算オペレータではなく、pd.offsets インスタンスの .apply メソッドを使う。.apply は Timestamp / datetime にも使える。(というか加減算オペレータも裏側では .apply を使っている)。
offsets.Week(weekday=2).apply(np.datetime64('2014-11-09 10:00:00Z')) # Timestamp('2014-11-12 10:00:00') offsets.Week(weekday=2).apply(dt) # Timestamp('2014-11-12 10:00:00')
ある日付が offsets に乗っているかどうか調べる場合は .onOffset。例えば ある日付が水曜日 ( weekday=2 ) かどうか調べたければ、
# 11/09は日曜日 ( weekday=6 ) dt # Timestamp('2014-11-09 10:00:00') pd.offsets.Week(weekday=2).onOffset(dt) # False pd.offsets.Week(weekday=2).onOffset(dt + offsets.Day(3)) # True
タイムゾーン
タイムゾーン関連の処理は tz_localize と tz_convert 二つのメソッドで行う。引数は 文字列、もしくは pytz, dateutil.tz インスタンス (後述)。
上でつくった Timestamp について、現在の表示時刻 (10:00) をそのままにして タイムゾーンを付与する場合は tz_localize。
dt # Timestamp('2014-11-09 10:00:00') dt.tz_localize('Asia/Tokyo') # Timestamp('2014-11-09 10:00:00+0900', tz='Asia/Tokyo')
現在の表示時刻を世界標準時 (GMTで10:00) とみて タイムゾーン付与する場合は まず GMT に tz_localize してから tz_convert。
dt.tz_localize('GMT').tz_convert('Asia/Tokyo') # Timestamp('2014-11-09 19:00:00+0900', tz='Asia/Tokyo')
一度 タイムゾーンを付与したあとは、tz_convert でタイムゾーンを変更したり、
dttz = dt.tz_localize('Asia/Tokyo') dttz # Timestamp('2014-11-09 10:00:00+0900', tz='Asia/Tokyo') dttz.tz_convert('US/Eastern') # Timestamp('2014-11-08 20:00:00-0500', tz='US/Eastern')
タイムゾーンを削除したりできる。
tz_localize(None):Timestampのローカル時刻を引き継いで タイムゾーンを削除tz_convert(None):Timestampの GMT 時刻を引き継いで タイムゾーンを削除
dttz # Timestamp('2014-11-09 10:00:00+0900', tz='Asia/Tokyo') dttz.tz_localize(None) # Timestamp('2014-11-09 10:00:00') dttz.tz_convert(None) # Timestamp('2014-11-09 01:00:00')
また、tz_localize, tz_convert は pytz, dateutil.tz を区別なく扱える。標準の datetime.datetime では pytz と dateutil.tz でタイムゾーンの初期化方法が違ったりするのでこれはうれしい。
import pytz ptz = pytz.timezone('Asia/Tokyo') dt.tz_localize(ptz) # Timestamp('2014-11-09 10:00:00+0900', tz='Asia/Tokyo') from dateutil.tz import gettz dtz = gettz('Asia/Tokyo') dt.tz_localize(dtz) # Timestamp('2014-11-09 10:00:00+0900', tz='dateutil//usr/share/zoneinfo/Asia/Tokyo')
条件に合致する日時のリスト生成
単純な条件での生成
pd.date_range で以下で指定する引数の条件に合致した日時リストを一括生成できる。渡せるのは 下の4つのうち 3つ。例えば start, periods, freq を渡せば end は自動的に計算される。
start: 開始時刻end: 終端時刻periods: 内部に含まれる要素の数freq: 生成する要素ごとに変化させる時間単位/周期 ( 1時間ごと、1日ごとなど)。freqとして指定できる文字列は こちら。またpd.offsetsも使える (後述)。
返り値は上でも少しふれた DatetimeIndex になる。print された DatetimeIndex の読み方として、表示結果の 2行目が 生成された日時の始点 (11/01 10:00) と終端 (11/04 09:00)、3行目が生成されたリストの長さ (72コ)と周期 ( H = 1時間ごと) になる。
dtidx = pd.date_range('2014-11-01 10:00', periods=72, freq='H') dtidx # <class 'pandas.tseries.index.DatetimeIndex'> # [2014-11-01 10:00:00, ..., 2014-11-04 09:00:00] # Length: 72, Freq: H, Timezone: None
生成された DatetimeIndex から Timestamp を取得する場合、ふつうに要素選択するか、 .asobject.tolist() を使う。.asobject.tolist() の意味は、
asobject:DatetimeIndex内の要素をTimestampオブジェクトに明示的に変換 (メソッドでなくプロパティなので注意)tolist():DatetimeIndex自身と同じ中身のリストを生成
dtidx[1] # Timestamp('2014-11-01 11:00:00', offset='H') dtidx.asobject.tolist() # [Timestamp('2014-11-01 10:00:00', offset='H'), # Timestamp('2014-11-01 11:00:00', offset='H'), # ... # Timestamp('2014-11-04 08:00:00', offset='H'), # Timestamp('2014-11-04 09:00:00', offset='H')]
少し複雑 (n個おき) な条件での生成
また、freq に offsets を与えればもう少し複雑な条件でもリスト生成できる。ある日時から3時間ごとの Timestamp を生成したい場合は以下のように書ける。
pd.date_range('2014-11-01 10:00', periods=72, freq=offsets.Hour(3)) # <class 'pandas.tseries.index.DatetimeIndex'> # [2014-11-01 10:00:00, ..., 2014-11-10 07:00:00] # Length: 72, Freq: 3H, Timezone: None
さらに複雑な条件 (任意の条件) での生成
さらに複雑な条件を使いたい場合は、一度単純な条件でリスト生成し、必要な部分をスライシングして取得するのがよい。例えば 深夜 0時から5時間ごとに1時間目/4時間目の時刻をリストとして取得したい場合は以下のようにする。
これで dateutil.rrule でやるような複雑な処理も (大部分? すべて?) カバーできるはず。
dtidx = pd.date_range('2014-11-01 00:00', periods=20, freq='H') selector = [False, True, False, False, True] * 4 dtidx[selector].asobject.tolist() # [Timestamp('2014-11-01 01:00:00'), # Timestamp('2014-11-01 04:00:00'), # Timestamp('2014-11-01 06:00:00'), # Timestamp('2014-11-01 09:00:00'), # Timestamp('2014-11-01 11:00:00'), # Timestamp('2014-11-01 14:00:00'), # Timestamp('2014-11-01 16:00:00'), # Timestamp('2014-11-01 19:00:00')]
DatetimeIndex に対する一括処理
最後に、先のセクションで記載した 時刻 加減算, タイムゾーン処理の方法/メソッドは DatetimeIndex に対しても利用できる。全体について 何かまとめて処理したい場合は DatetimeIndex の時点で処理しておくと楽。
dtidx # <class 'pandas.tseries.index.DatetimeIndex'> # [2014-11-01 00:00:00, ..., 2014-11-01 19:00:00] # Length: 20, Freq: H, Timezone: None # 1日後にずらす dtidx + offsets.Day() # <class 'pandas.tseries.index.DatetimeIndex'> # [2014-11-02 00:00:00, ..., 2014-11-02 19:00:00] # Length: 20, Freq: H, Timezone: None # タイムゾーンを設定 dtidx.tz_localize('Asia/Tokyo') # <class 'pandas.tseries.index.DatetimeIndex'> # [2014-11-01 00:00:00+09:00, ..., 2014-11-01 19:00:00+09:00] # Length: 20, Freq: H, Timezone: Asia/Tokyo
まとめ
pandasなら日時/タイムスタンプ関連の操作もカンタン。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (9件) を見る
Python lifelines で生存分析
サイトのアクセス履歴をみていたら "Python", "生存分析", "Kaplan-Meier" なんかがちらほらあったので、知っている方法を書いてみる。
生存分析とは
いくつかのサンプルについて、何らかのイベントが起きるまでの時間とイベント発生率との関係をモデル化するもの。よく使われるのは
- 何かの病気になってからの死亡率
- 部品の故障率
など。
この説明がわかりやすい。
Python でやる方法
どちらか。lifelines のほうが手間はすくない。
- lifelines パッケージを使う
- rpy2 でデータを R に渡して survival パッケージ を使う
1. lifelines パッケージを使う
これを使えば Pure-Python で生存分析できる。使えるのは、以下のノンパラ、セミパラ系の手法。Weibull などのパラメトリック系はない (Scipy 使えばできる)。
- Kaplan-Meier
- Nelson-Aalen
- Cox 比例ハザード
インストール
pip install lifelines
ドキュメント
バージョン古いようなのでアップデートされてないかも。
Lifelines — lifelines 0.8.0.1 documentation
やり方
データは pandas の DataFrame で扱うのがよい。ここでは、lifelinesにもサンプルデータとして入っている lung (肺ガン) データを使う。
lungデータ
https://github.com/CamDavidsonPilon/lifelines/blob/master/datasets/lung.csv
データの意味はこちら。
ここでは、性別でグループ分けしたデータについて、イベントの発生時間 ("time" カラムの値) とそのときの状態 ("death" カラムの値, Trueなら死亡, Falseなら censored = 観測打ち切り) で Kaplan-Meierモデルを作る。
これは lifelines.KaplanMeierFitter クラスでできる。引数の詳細はドキュメント参照。
Introduction to using lifelines
import pandas as pd # 表示する行数を指定 pd.options.display.max_rows = 5 # データの読み込み df = pd.read_csv("lung.csv") # 元データには、死亡なら2, censored なら1が入っている # lifelines に渡すため、死亡なら True, censored なら False となるよう変換する df['death'] = df['status'] == 2 df # Unnamed: 0 inst time status age sex ph.ecog ph.karno pat.karno \ # 0 1 3 306 2 74 1 1 90 100 # 1 2 3 455 2 68 1 0 90 90 # .. ... ... ... ... ... ... ... ... ... # 226 227 6 174 1 66 1 1 90 100 # 227 228 22 177 1 58 2 1 80 90 # # meal.cal wt.loss death # 0 1175 NaN True # 1 1225 15 True # .. ... ... ... # 226 1075 1 False # 227 1060 0 False # # [228 rows x 12 columns] from lifelines import KaplanMeierFitter # グラフを描画する Axes ax = None # 性別でグループ分け for name, group in df.groupby('sex'): kmf = KaplanMeierFitter() kmf.fit(group['time'], event_observed=group['death'], label = 'sex=' + str(name)) # 描画する Axes を指定。None を渡すとエラーになるので場合分け if ax is None: ax = kmf.plot() else: ax = kmf.plot(ax=ax) # Kaplan-Meier曲線を描画 import matplotlib.pyplot as plt plt.show()
Kaplan-Meier曲線 (横軸が経過時間、縦軸が生存率 = イベントの未発生率) となる。性別 = 2 (女性) のほうが 同じ経過時間時点での生存率が高い傾向がある。
2. rpy2 でデータを R に渡して survival パッケージ を使う
pandas と rpy2 との連携はこちらを参照。
やり方
やりたいことにもよるが、結果を Python に戻してからの処理がちょっとめんどい。
import numpy as np import pandas as pd # 表示する行数を指定 pd.options.display.max_rows = 5 # rpy2 使って survival::lung データセットを呼び出す # 手元のデータを使うなら不要 import rpy2.robjects as robjects import pandas.rpy.common as com robjects.packages.importr('survival') lung = com.load_data('lung') # rpy2 形式に変換し、 Rへ渡す lungr = com.convert_to_r_dataframe(lung) robjects.r.assign('lungr', lungr); # R で生存分析実施 robjects.r('fitted <- survfit(Surv(time, status) ~ sex, data = lungr)'); # 結果を Python に戻し、各プロパティを numpy.arrayとして取得 fitted = robjects.r['fitted'] time = np.array(fitted.rx2('time')) surv = np.array(fitted.rx2('surv')) lower = np.array(fitted.rx2('lower')) upper = np.array(fitted.rx2('upper')) # strata のキーを順序どおりに取得するには一手間必要 robjects.r('strata_keys <- names(fitted$strata)') strata_keys = com.convert_robj(robjects.r['strata_keys']) strata = com.convert_robj(fitted.rx2('strata')) strata = reduce(lambda x, y: x + y, [ [k] * strata[k] for k in strata_keys]) # DataFrame作成 result = pd.DataFrame({'time': time, 'surv': surv, 'lower': lower, 'upper': upper, 'strata': strata}) result # lower strata surv time upper # 0 0.954230 sex=1 0.978261 11 1.000000 # 1 0.943423 sex=1 0.971014 12 0.999412 # .. ... ... ... ... ... # 204 0.001582 sex=2 0.011111 821 0.078022 # 205 NaN sex=2 0.000000 965 NaN # プロット import matplotlib.pyplot as plt # グラフを描画する Axes ax = None # 生存曲線と 信頼区間を別々に描画する必要がある # strata ごとに同じ色を使うために グループ名 : 色の辞書を作成 colors = {'sex=1': 'green', 'sex=2': 'blue'} for name, group in result.groupby('strata'): # pandas は index を x軸として plot するので、 # time の値を indexに設定 group = group.set_index('time') color = colors[name] # pandas.Series.plot で生存曲線自体を描画 ax = group['surv'].plot(ax=ax, color=colors[name]) # matplotlib.Axes.fill_between で信頼区間を描画 ax.fill_between(group.index, group['lower'], group['upper'], color=colors[name], alpha=0.3) plt.show()
