pandasでクロス集計

テーブルデータを分析する時、2つの列の値の関係を調べたいことが良くあります。
対象の列がどちらも連続値をとるのであれば、散布図を書きますが、
これが連続値ではなく性別や都道府県、世代、アンケートのn段階評価などの時は、散布図はあまり適しません。

このような場合は、クロス集計を使うと便利です。
(wikipediaでは分割表と呼んでいるようです。)

pandasにクロス集計を行う専用の関数があるのでそれを紹介します。
(職場ではcsvに書き出してtableauに読み込ませてやることも多いのですが、pythonで完結できればそれはそれ便利です。)

ドキュメントはこちら。
pandas.crosstab

aggfuncに値を渡せば、数え上げ以外にも色々できることがあり、pivot_table的な使い方もできるのですが、
とりあえずダミーデータを用意してデータの数え上げでクロス表を作ってみましょう。
ABテストなどで頻繁に使いますからね。


import pandas as pd
import numpy as np
# ダミーデータ生成
df = pd.DataFrame(
        {
            "data1": np.random.randint(1, 6, 100),
            "data2": np.random.choice(["foo", "bar"], 100),
        }
    )
# サンプル(5件)
print(df.sample(5))
# 出力
'''
    data1 data2
71      2   bar
13      1   foo
7       5   bar
16      1   bar
56      2   foo
'''

# クロス集計
cross_df = pd.crosstab(df.data2, df.data1)
print(cross_df)
# 出力
'''
data1   1   2   3   4  5
data2                   
bar    10  13   8  12  8
foo    12   9  11   9  8
'''

とても楽にクロス集計ができました。
SQLでこれを実装しようとするとかなりの手間なので、地味に重宝しています。

完全に蛇足ですが、crosstabを知らなかった頃はこういうコードを書いていました。
indexとcolumnsだけ指定して値が全部0のデータフレームを作って、
元のデータの値を順番に数え上げています。
今思えばcrosstab知らなくても、もう少しまともな書き方がありそうです。


cross_df_2 = pd.DataFrame(
        0,
        index=set(df.data2),
        columns=set(df.data1),
    )
for _, row in df.iterrows():
    cross_df_2.loc[row.data2, row.data1] += 1    
print(cross_df_2)
# 出力
'''
     1   2   3   4  5
foo  12   9  11   9  8
bar  10  13   8  12  8
'''

pythonで一般化線形モデル(GLM) – ポアソン回帰

以前の記事で、statsmodelsを使って線形回帰をやったので、今回はポアソン回帰をやってみます。
参考:statsmodelsで重回帰分析

データは久保拓弥先生の、データ解析のための統計モデリング入門 (通称緑本)の第3章から拝借し、
本に載っているのと同じ結果を得ることを目指します。
ちなみに本のコードはRで実装されています。
Rは勉強中なので写経はしましたがそれはそれとして、
今回はいつも使っているpythonのstatsmodelsでやってみます。

データはこちらのサポートページからダウンロードできます。
生態学のデータ解析 – 本/データ解析のための統計モデリング入門
3章の data3a.csv を保存しておきましょう。

このデータは、架空の植物の種子の数に関するもので、
種子の数が$y_i$, 植物のサイズが$x_i$, 施肥処理を行ったかどうかが$f_i$列に格納されています。

そして、種子の個数$y_i$が期待値 $\lambda_i=\exp(\beta_1+\beta_2x_i+\beta_3f)$のポアソン分布に従うと仮定した場合の尤度を最大化するように係数を決定します。
($f$は施肥処理を行ったら1, 行ってない場合は0)

早速やってみましょう。
ドキュメントはこのあたりが参考になります。
statsmodels.genmod.generalized_linear_model.GLM
Model Family に Poisson を指定すると、リンク関数は自動的にlogが選ばれます。


import pandas as pd
import statsmodels.api as sm

# データの読み込み
# ファイルの取得元
# http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html
df = pd.read_csv("./data3a.csv")

# f 列の値を Tならば 1 , その他(Cのみ)ならば0に符号化
df["fT"] = (df.f == "T").astype(int)

y = df.y
X = df[["x", "fT"]]
# 定数列(1)を作成
X = sm.add_constant(X)

# モデル生成と学習
model = sm.GLM(y, X, family=sm.families.Poisson())
result = model.fit()

# 結果出力
print(result.summary())

# 以下出力結果
                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       97
Model Family:                 Poisson   Df Model:                            2
Link Function:                    log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -235.29
Date:                Sun, 21 Apr 2019   Deviance:                       84.808
Time:                        16:36:24   Pearson chi2:                     83.8
No. Iterations:                     4   Covariance Type:             nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.2631      0.370      3.417      0.001       0.539       1.988
x              0.0801      0.037      2.162      0.031       0.007       0.153
fT            -0.0320      0.074     -0.430      0.667      -0.178       0.114
==============================================================================

本のP58に載っているのと全く同じ係数を得ることができました。

ついでに実際のデータとこのモデルから予測される種子数をプロットして可視化してみます。


import matplotlib.pyplot as plt
import numpy as np

# プロット用のデータ作成
xx = np.linspace(7, 13, 101)
yy0 = np.exp(result.params["const"] + result.params["x"]*xx)
yy1 = np.exp(result.params["const"] + result.params["x"]*xx + result.params["fT"])

# 可視化
fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel("個体のサイズ")
ax.set_ylabel("種子の個数")
ax.scatter(X[X.fT == 0]["x"], y[X.fT == 0].values, marker="x", label="肥料無し")
ax.scatter(X[X.fT == 1]["x"], y[X.fT == 1].values, alpha=0.7, label="肥料有り")
ax.plot(xx, yy0, linestyle="--",  label="肥料無し")
ax.plot(xx, yy1, label="肥料有り")
plt.legend()
plt.show()

出力結果がこちら。

pandasで散布図行列を書く

何かしらのデータについて調べる時、各データごとの関係をまとめて可視化するために
散布図行列(matrix of scatter plots)を描きたくなることがあります。

matplotlibでゴリゴリ実装するのは手間なので、seabornが良く使われていると思うのですが、
実はpandasでも作成できるのでその紹介です。
(seabornを使えばいいじゃん、というのはごもっともなのですが、
なぜか僕の自宅環境ではwarningが発生し、その原因調査中なのであまり使いたくないのです。)

ドキュメントはこちら。
pandas.plotting.scatter_matrix

早速ですが、irisでやってみましょう。


import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

iris = load_iris()
df = pd.DataFrame(
        iris.data,
        columns=iris.feature_names,
    )
pd.plotting.scatter_matrix(df, figsize=(15, 15), c=iris.target)
plt.show()

そして表示される結果がこちらです。

オプションで、 引数cにクラス情報を渡して色もつけました。
このコーディング量でかける図としては十分なクオリティではないでしょうか。

statsmodelsで重回帰分析

前回に引き続きpythonで重回帰分析を行い、回帰係数と切片を求める方法の紹介です。
今回はstatsmodelsを使います。

ドキュメントはこちら。
statsmodels.regression.linear_model.OLS

データは同じものを使い、結果が一致することを確認したいので
保存してたものを読み込みます。


import numpy as np
import statsmodels.api as sm

# データの読み込み
npzfile = np.load("sample_data.npz")
X = npzfile["X"]
y = npzfile["y"]

つぎに、回帰分析に入りますが、statsmodelsで回帰分析する場合には、一点注意があります。
それは、上で読み込んだ Xとyをそのままモデルに食わせると、切片(定数項)を0として結果を返してくることです。
それを回避するために、 X に対して、全ての値が1の列を追加する必要があります。
それは、次の関数を使って実行しています。
statsmodels.tools.tools.add_constant


# 配列Xに1だけからなる列を追加する。
X0 = sm.add_constant(X)
# 回帰分析実行
st_model = sm.OLS(y, X0)
results = st_model.fit()

# 切片と回帰係数
print(results.params)
# [ 4.94346432  2.00044153 -2.99255801  3.98231315  0.01708309]

無事にscikit-learnで行った時と全く同じ結果を得ることができましたね。
これだけだと、add_constantをやらなくてよかった分scikit-learnの圧勝に見えますが、
statsmodels を使うと一つメリットがあります。
それは、各係数が0でないことの検定ができることです。
resultsの sumary()って関数を実行してみましょう。


print(results.summary())

# 以下出力です。
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.999
Model:                            OLS   Adj. R-squared:                  0.999
Method:                 Least Squares   F-statistic:                 2.136e+04
Date:                Mon, 15 Apr 2019   Prob (F-statistic):          2.24e-139
Time:                        00:59:22   Log-Likelihood:                -143.94
No. Observations:                 100   AIC:                             297.9
Df Residuals:                      95   BIC:                             310.9
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9435      0.108     45.726      0.000       4.729       5.158
x1             2.0004      0.020    101.534      0.000       1.961       2.040
x2            -2.9926      0.018   -166.158      0.000      -3.028      -2.957
x3             3.9823      0.017    231.575      0.000       3.948       4.016
x4             0.0171      0.019      0.914      0.363      -0.020       0.054
==============================================================================
Omnibus:                        2.170   Durbin-Watson:                   1.994
Prob(Omnibus):                  0.338   Jarque-Bera (JB):                1.862
Skew:                          -0.208   Prob(JB):                        0.394
Kurtosis:                       2.477   Cond. No.                         7.12
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

出力の真ん中にある通り、定数項および回帰係数のt値と、p値が出力されています。
そして、x4 は P値が 0.05より大きいので、x4の係数が0であるという帰無仮説を棄却できていないことがわかります。
もともと、x4の回帰係数は0としてサンプルデータを作っていたので、これは良い結果です。
そのほか、AICなどの値も出してくれています。

scikit-learnでこれらの情報を収集するのは手間なので、
目的によってはstatsmodelsが便利な場面もあります。

scikit-learnで重回帰分析

今回と次回でpythonで重回帰分析を実行する方法を二つ紹介します。
今回はscikit-learnのLinearRegressionを使う方法です。

ドキュメントはこちら。
sklearn.linear_model.LinearRegression

最初に検証用のダミーデータを作ります。
$x_{i,j}$を -10 ~ 10の一様分布からサンプリングし、次の式で$y_i$を作ります。
$x_{i,3}$の係数が0になっていることから分かる通り、$x_{i,3}$は$y_i$には何の関係もない値です。
また、ノイズとして正規分布に従う乱数加えておきます。
$$
y_i = 5 + 2x_{i,0} -3x_{i,1} + 4x_{i,2} + 0x_{i,3} + \varepsilon_i,\\
\varepsilon_i \sim N(0,1) , \ \ \ i = 0,1,\cdots, 99
$$

サンプルデータを作って保存するコードがこちら。


import numpy as np
X = np.random.uniform(-10, 10, size=(100, 4))
y = 5 + X@[2, -3, 4, 0] + np.random.normal(0, 1, 100)
np.savez("sample_data.npz", X=X, y=y)

早速、回帰分析して回帰係数と定数項ついでに決定係数を求めてみましょう。
(回帰分析の目的が予測モデルを作ることであれば、データを訓練用と評価用に分けるべきなのですが、
今回は回帰係数を求める方法の紹介が主目的なので分けていません。)


import numpy as np
from sklearn.linear_model import LinearRegression

# データの読み込み
npzfile = np.load("sample_data.npz")
X = npzfile["X"]
y = npzfile["y"]

# モデルのインスタンス生成
model = LinearRegression()
#学習
model.fit(X, y)
# LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

# 回帰係数
print(model.coef_)
# [ 2.00044153 -2.99255801  3.98231315  0.01708309]

# 切片
print(model.intercept_)
# 4.943464324307898

# 決定係数 R^2
model.score(X, y)
# 0.9988893054003364

各コードの後ろにつけているコメントが僕の環境で実行した時の結果です。
回帰係数も切片もそれぞれほぼ正解に近い値が算出されていますね。

注:サンプルデータを乱数で生成しているので、データ生成からやり直せば結果は変わります。

scikit-learnを使うと、非常に手軽に回帰分析ができることがわかりました。
次回はstatsmodelsで同じことをやってみます。

numpyの配列をファイルに保存する

日常の実際で必要になったことはないのですが、
numpyに配列(ndarray)をファイルに保存する機能があるのでその紹介です。

(実用上、ファイルに保存したくなるようなデータは pandasのデータフレームに変換したあと、
to_csvやto_excelを使ってcsvかエクセルにしてしまうことが多いです。)

ドキュメントはこの辺り
numpy.save
numpy.load
numpy.savez

簡単に言ってしまえば、
numpy.save(ファイル名, 保存したい配列) で保存して、
numpy.load(ファイル名) で読み込める。
numpy.savez(ファイル名, (名前=)保存したい配列, (名前=)保存したい配列, …) で、
1ファイルに複数配列保存できる、と言った使い方になります。(名前=)は省略可能。

実際にやってみましょう。
まず検証用にダミーデータを作ります。


import numpy as np

# 保存のテスト用データを作成
data0 = np.arange(12).reshape(3, 4)
data1 = data0 ** 2
print(data0)
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''
print(data1)
'''
[[  0   1   4   9]
 [ 16  25  36  49]
 [ 64  81 100 121]]
'''

まずは、save関数で配列を一つ保存し、それを別の変数に読み込んでみましょう。
特に難しいことはなく非常に直感的な使い方で上手く動きます。


# .npyフォーマットで保存
np.save("data0.npy", data0)

# 読み込み
data0_load = np.load("data0.npy")
print(data0_load)
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''

続いて、savez関数で複数配列を1ファイルに保存し、復元してみます。
こちらは、loadした時に返されるオブジェクトが保存した配列そのものではないので注意です。
まずは名前をつけずに保存するパターン。


# savez を使うと複数の配列を.npzフォーマットの1ファイルにまとめて保存できる
np.savez("data0_1.npz", ary0, ary1)

# .npzファイルを読み込むと、配列ではなく、NpzFile オブジェクトが返される。
npzfile = np.load("data0_1.npz")
# .filesプロパティを見ると、中に二つの配列を持っていることがわかる
print(npzfile.files)
# ['arr_0', 'arr_1']

# それぞれの配列の取り出し
print(npzfile["arr_0"])
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''
print(npzfile["arr_1"])
'''
[[  0   1   4   9]
 [ 16  25  36  49]
 [ 64  81 100 121]]
'''

また、名前をつけて保存すると次のようになります。


# 名前(例ではxとy)をつけて保存することも可能
np.savez("named_data0_1.npz", x=ary0, y=ary1)
# .npzファイルを読み込むと、配列ではなく、NpzFile オブジェクトが返される。
npzfile2 = np.load("named_data0_1.npz")
# .filesプロパティを見ると、中に二つの配列を持っていることがわかる
print(npzfile2.files)
# ['x', 'y']

# 保存時の名前で取り出すことが可能
print(npzfile2["x"])
'''
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
'''
print(npzfile2["y"])
'''
[[  0   1   4   9]
 [ 16  25  36  49]
 [ 64  81 100 121]]
'''

勝手にarr_0とか命名されるよりも、自分で名前をつけておいた方が混乱がなさそうですね。

pythonでMD5

業務でMD5を使って文字列をハッシュ化する機会があったのでそのメモです。
(最近ではMD5自体があまり安全では無く、暗号学的な強度が必要なときは使うべきでは無いのでご注意ください。)

python3では、標準ライブラリにhashlibというものがあり、これにMD5が実装されています。
ドキュメント: hashlib — セキュアハッシュおよびメッセージダイジェスト

使い方は簡単で、インポートして hashlib.md5(“text”) するだけです、
となれば楽なのですが、これでは動きません。こんなふうにエラーになります。


import hashlib
hashlib.md5("text")

# 以下出力
TypeError                                 Traceback (most recent call last)
 in ()
----> 1 hashlib.md5("text")

TypeError: Unicode-objects must be encoded before hashing

これは、bytes型でデータを渡す必要があるからです。
そのため、b"text""text".encode("utf-8")のようにする必要があります。

これで一応エラーは無くなりますが、まだ少し不十分です。
というのも、戻り値がHASH object という型で返ってきます。


import hashlib
text = "Hello, World!"
bytes_text = text.encode("utf-8")
print(hashlib.md5(bytes_text))

# 出力
<md5 HASH object @ 0x10e3f7620>

欲しいのは、 16進法で表記された文字列だと思うので、もう一手間かけて、hexdigestという関数を実行します。


import hashlib
text = "Hello, World!"
bytes_text = text.encode("utf-8")
hash_obj = hashlib.md5(bytes_text)
print(hash_obj.hexdigest())

# 出力
65a8e27d8879283831b664bd8b7f0ad4

# 1行でやる方法
print(hashlib.md5(b"Hello, World!").hexdigest())

# 出力は同じ
65a8e27d8879283831b664bd8b7f0ad4

これでできました。

sha1やsha256もhashlibに実装されているので同じように使えます。

pythonで累積和

数列に対して、最初の項からn番目の項まで足したものの数列が必要になることがあります。
日々の売上データからその日までの累計の売上データを出すような計算ですね。

イメージとしては、
1,2,3,4,5,…に対して、1,3,6,10,15,… を得るような操作です。

スライスと内包表記を組み合わせてもいいし、足し算くらいならfor文を回しても問題ないのですが、
numpyやpandas に専用の関数が用意されているのでその紹介です。

ドキュメントはこの辺
numpy.cumsum
pandas.Series.cumsum

早速実行してみました。


import numpy as np
import pandas as pd

ary = np.arange(1, 11)
print(ary)
# [ 1  2  3  4  5  6  7  8  9 10]

print(ary.cumsum())
# [ 1  3  6 10 15 21 28 36 45 55]

series = pd.Series(np.arange(1, 11))
print(series)
'''
0     1
1     2
2     3
3     4
4     5
5     6
6     7
7     8
8     9
9    10
dtype: int64
'''

print(series.cumsum())
'''
0     1
1     3
2     6
3    10
4    15
5    21
6    28
7    36
8    45
9    55
dtype: int64
'''

よく似た関数として、
numpyには累積の積を返すcumprodが用意されており、
さらにpandasの方にはcumprodに加えて、そこまでの最大値と最小値を得られるcummax, cumminが用意されています。

JSON文字列をオブジェクトに変換する

昨日の記事の逆です。(まとめて書けばよかった)

ドキュメントも同じ場所。
json — JSON エンコーダおよびデコーダ

JSON型の文字列でデータが渡されたとき、それをオブジェクトに変換するには、
json.loads という関数を使います。

早速やってみましょう。
printするだけだと型が変わっていることがわかりにくいので、
typeの結果も出しました。


import json
json_data = '{"key1": "value1", "key2": "value2", "ary1": ["element1", "element2"]}'
data = json.loads(json_data)

print(data)
print(type(data))

# 出力
# {'key1': 'value1', 'key2': 'value2', 'ary1': ['element1', 'element2']}
# <class 'dict'>

辞書型のオブジェクトをJSON文字列に変換する

仕事で必要になったのでメモです。

pythonには jsonという標準ライブラリがあり、
それを使うことで、配列や辞書型のオブジェクトをJSON文字列に簡単に変換することができます。

ドキュメントはこちら。
json — JSON エンコーダおよびデコーダ

インポートして dumps関数を使うだけなので早速やってみましょう。


import json
data = {
        "key1": "value1",
        "key2": "value2",
        "ary1": ["element1", "element2"],
    }
json_data = json.dumps(data)
print(json_data)

# 出力
# {"key1": "value1", "key2": "value2", "ary1": ["element1", "element2"]}

ちなみに、json.dumpsを使わず、str(data)で文字列に変換すると、結果は
"{'key1': 'value1', 'key2': 'value2', 'ary1': ['element1', 'element2']}"
になります。
JSONのルールでは文字列はダブルクオーテーションで囲まないといけないので、
これはJSONではありません。

厳密にJSON型を要求する関数やAPIに渡すときはこれだと受け付けられないので、
json.dumpsを使いましょう。