この記事はGunosy Advent Calendar 2014の23日目の記事です。
松永と申します。エンジニアのアルバイトをしています。
Gunosyでは隔週でGunosyデータマイニング研究会が開催されており、 現在輪読にはデータ解析のための統計モデリング入門が使用されています。 「緑本」と呼ばれているようです(以下、この本のことを緑本と呼びます)
緑本では学んだ統計知識の実践例としてRによる実行結果が掲載されています。しかし、私はRをほとんど使ったことがないうえ、普段はPythonをメインに利用しているため、今回はPythonでやってみました。今回はIpython NoteBookを利用して、実行と解説を平行して行っていきます。 今回主に利用したモジュールは以下のとおりです。
おそらく依存関係等で他にインストールが必要なモジュールはあると思いますが、適宜よろしくお願いします。
緑本の第3章のみとなります。 第3章では、用意されたサンプルデータを使ってGLM(Generalized Linear Model, 一般線形化モデル)が紹介されています。サンプルデータの利用が3章から始まるということで3章を選択しました。行った演習内容としては基本的な統計量を取得したり分布図をプロットしたりするといった基本的なものになっています。
サンプルデータは種子データで、yはある植物の種子の数、xはその種子を持っていた植物のサイズ、fは使った肥料の種類を表しています。植物のサイズや使った肥料の種類が生存した種子の数にどのような影響を与えているのかモデリングしていくという流れになっています。 x, yは数値データですが、fは肥料の種類なので名義尺度となっています。
まずはデータの読み込みから。(p42)
サンプルデータは著者である久保先生のウェブサイトからダウンロードすることができます。
import pandas as pd d = pd.read_csv('data3a.csv', names=['y', 'x', 'f'])
pandasにはDataFrameというクラスがあり、以下のようにデータを綺麗に整形してくれます。 pandasのread_csvを利用した場合、読み込まれたデータはDataFrameクラスのインスタンスとなります。
d
y | x | f | |
---|---|---|---|
0 | y | x | f |
1 | 6 | 8.31 | C |
2 | 6 | 9.44 | C |
3 | 6 | 9.5 | C |
4 | 12 | 9.07 | C |
5 | 10 | 10.16 | C |
6 | 4 | 8.32 | C |
7 | 9 | 10.61 | C |
8 | 9 | 10.06 | C |
9 | 9 | 9.93 | C |
10 | 11 | 10.43 | C |
11 | 6 | 10.36 | C |
12 | 10 | 10.15 | C |
13 | 6 | 10.92 | C |
14 | 10 | 8.85 | C |
15 | 11 | 9.42 | C |
16 | 8 | 11.11 | C |
17 | 3 | 8.02 | C |
18 | 8 | 11.93 | C |
19 | 5 | 8.55 | C |
20 | 5 | 7.19 | C |
21 | 4 | 9.83 | C |
22 | 11 | 10.79 | C |
23 | 5 | 8.89 | C |
24 | 10 | 10.09 | C |
25 | 6 | 11.63 | C |
26 | 6 | 10.21 | C |
27 | 7 | 9.45 | C |
28 | 9 | 10.44 | C |
29 | 3 | 9.44 | C |
... | ... | ... | ... |
71 | 10 | 10.54 | T |
72 | 8 | 11.3 | T |
73 | 8 | 12.4 | T |
74 | 7 | 10.18 | T |
75 | 5 | 9.53 | T |
76 | 6 | 10.24 | T |
77 | 8 | 11.76 | T |
78 | 9 | 9.52 | T |
79 | 9 | 10.4 | T |
80 | 6 | 9.96 | T |
81 | 7 | 10.3 | T |
82 | 10 | 11.54 | T |
83 | 6 | 9.42 | T |
84 | 11 | 11.28 | T |
85 | 11 | 9.73 | T |
86 | 11 | 10.78 | T |
87 | 5 | 10.21 | T |
88 | 6 | 10.51 | T |
89 | 4 | 10.73 | T |
90 | 5 | 8.85 | T |
91 | 6 | 11.2 | T |
92 | 5 | 9.86 | T |
93 | 8 | 11.54 | T |
94 | 5 | 10.03 | T |
95 | 9 | 11.88 | T |
96 | 8 | 9.15 | T |
97 | 6 | 8.52 | T |
98 | 8 | 10.24 | T |
99 | 7 | 10.86 | T |
100 | 9 | 9.97 | T |
101 rows × 3 columns
以下のようにカラムを指定すると特定のカラムのデータだけを取り出すことができます。
d['x']
0 x 1 8.31 2 9.44 3 9.5 4 9.07 5 10.16 6 8.32 7 10.61 8 10.06 9 9.93 10 10.43 11 10.36 12 10.15 13 10.92 14 8.85 ... 86 10.78 87 10.21 88 10.51 89 10.73 90 8.85 91 11.2 92 9.86 93 11.54 94 10.03 95 11.88 96 9.15 97 8.52 98 10.24 99 10.86 100 9.97 Name: x, Length: 101, dtype: object
d['y']
0 y 1 6 2 6 3 6 4 12 5 10 6 4 7 9 8 9 9 9 10 11 11 6 12 10 13 6 14 10 ... 86 11 87 5 88 6 89 4 90 5 91 6 92 5 93 8 94 5 95 9 96 8 97 6 98 8 99 7 100 9 Name: y, Length: 101, dtype: object
d['f']
0 f 1 C 2 C 3 C 4 C 5 C 6 C 7 C 8 C 9 C 10 C 11 C 12 C 13 C 14 C ... 86 T 87 T 88 T 89 T 90 T 91 T 92 T 93 T 94 T 95 T 96 T 97 T 98 T 99 T 100 T Name: f, Length: 101, dtype: object
d['x'].dtype
dtype('O')
d['y'].dtype
d['f'].dtype
dtype('O')
int64などのデータ型が表示されると思うのですが...気にせず次に行きましょうw (p43, p44)
d.describe()
y | x | f | |
---|---|---|---|
count | 101 | 101 | 101 |
unique | 15 | 92 | 3 |
top | 6 | 10.21 | T |
freq | 20 | 2 | 50 |
d['x'].describe()
count 101 unique 92 top 10.21 freq 2 Name: x, dtype: object
describe関数はRのsummaryに相当します。 あれ?なぜか表示項目が少ないですね。 一旦他の例で試してみます。
import numpy as np
series = pd.Series(np.random.random_sample(100))
series.describe()
count 100.000000 mean 0.448686 std 0.299202 min 0.005042 25% 0.184654 50% 0.416622 75% 0.707822 max 0.992557 dtype: float64
これはRのsummaryとほぼ同じですね!! なぜcsvから読み込んだデータでは頻度に関する指標しか見られなかったのでしょうか... クラスはいずれも
type(d['x'])
pandas.core.series.Series
type(series)
pandas.core.series.Series
...答えは単純でした。私のデータ読み込み方法が間違っていました。
csvファイルの一行目にデータの項目名が入っている場合 → read_csvでnamesを指定する必要はない
csvファイルの一行目からデータが始まっている場合 → read_csvでnamesを指定
今回使用しているサンプルデータの一行目にはデータ名が入力されているためnamesは必要ないです。 namesを引数に指定すると、1行目もデータの一部であると解釈され、数字と文字が混ざったデータになってしまい、思ったような結果を得ることができません。 以下のように修正して再読み込みしてみます。
d = pd.read_csv('data3a.csv')
d.describe()
y | x | |
---|---|---|
count | 100.000000 | 100.000000 |
mean | 7.830000 | 10.089100 |
std | 2.624881 | 1.008049 |
min | 2.000000 | 7.190000 |
25% | 6.000000 | 9.427500 |
50% | 8.000000 | 10.155000 |
75% | 10.000000 | 10.685000 |
max | 15.000000 | 12.400000 |
おおおおお!!Rのsummaryと同じことが出来ました!! 次に、先ほど結果が怪しかったd['x'].dtypeを実行してみましょう。
d['x'].dtype
dtype('float64')
ばっちりですね! 実際にはcsvを読み込んでそのまま解析できるというケースは少ないとは思います。 pandasはデータクリーニングが行えるとのことなので、事前に今回のようなミスは避けるようにしたいところ。
次に、データをプロットしてみます (p46~)
import matplotlib.pyplot as plt
以下の一行はIpython notebookでグラフプロットを表示するためのものです。
%matplotlib inline
Ts = d[d['f']=='T']
Cs = d[d['f']=='C']
plt.xlabel('x', size=20)
plt.ylabel('y', size=20)
T = plt.scatter(Ts['x'], Ts['y'], color='b', label='T')
C = plt.scatter(Cs['x'], Cs['y'], color='r', label='C')
plt.legend((T, C), ('T', 'C'))
<matplotlib.legend.Legend at 0x10b32dc50>
少しRで分布図を書く時よりコードが増えてしまいましたが、本文中と同じ分布図がかけたと思います 次に箱ひげ図を書いてみましょう
plt.boxplot([Cs['y'], Ts['y']])
{'boxes': [<matplotlib.lines.Line2D at 0x10b463350>, <matplotlib.lines.Line2D at 0x10b47dfd0>], 'caps': [<matplotlib.lines.Line2D at 0x10b471310>, <matplotlib.lines.Line2D at 0x10b471950>, <matplotlib.lines.Line2D at 0x10b487f10>, <matplotlib.lines.Line2D at 0x10b491590>], 'fliers': [<matplotlib.lines.Line2D at 0x10b47d610>, <matplotlib.lines.Line2D at 0x10b49e250>], 'means': [], 'medians': [<matplotlib.lines.Line2D at 0x10b471f90>, <matplotlib.lines.Line2D at 0x10b491bd0>], 'whiskers': [<matplotlib.lines.Line2D at 0x10b4635d0>, <matplotlib.lines.Line2D at 0x10b463c90>, <matplotlib.lines.Line2D at 0x10b487290>, <matplotlib.lines.Line2D at 0x10b4878d0>]}
上手くかけてるみたいです。
それでは、次に一般線形モデルの話に移ります。
一般化線形モデル(Generalized Linear Model, GLM) によって今回使用しているデータをモデル化してみます。 GLMの詳しい解説は本に任せます←
「GLMもscipyが関数を用意してくれているので、それに当てはめるだけでOK!」 と信じていたのですが、Rが用意しているGLMとは異なるようなので、今回は「statsmodels」というライブラリを使ってみます。 インストールはpipを使っている人なら pip install statsmodels でOKです。早速使ってみましょう。
実は、statsmodelは二種類のGLMを用意していて、以下の二種類があります。
from statsmodels.api import GLM
from statsmodels.formula.api import glm
違いをすべて理解できたわけではありませんが、 1つ目は目的変数と説明変数を明示的に指定して利用するものです。今回、肥料を示す変数は名義変数となっているため、利用できませんでした。 2つ目は自分で目的変数と説明変数の関係を明示して利用するものです。今回はこちらを利用します。関数が取る引数もこちらのほうがRのglm関数と似ていますね。
今回は説明変数が種子のサイズだけの場合の統計モデルを扱います。(p49)
import statsmodels.api as sm
model = glm('y~x', data=d, family=sm.families.Poisson(sm.families.links.log)).fit()
model.summary()
Dep. Variable: | y | No. Observations: | 100 |
---|---|---|---|
Model: | GLM | Df Residuals: | 98 |
Model Family: | Poisson | Df Model: | 1 |
Link Function: | log | Scale: | 1.0 |
Method: | IRLS | Log-Likelihood: | -235.39 |
Date: | Tue, 16 Dec 2014 | Deviance: | 84.993 |
Time: | 17:43:31 | Pearson chi2: | 83.8 |
No. Iterations: | 7 |
coef | std err | z | P>|z| | [95.0% Conf. Int.] | |
---|---|---|---|---|---|
Intercept | 1.2917 | 0.364 | 3.552 | 0.000 | 0.579 2.005 |
x | 0.0757 | 0.036 | 2.125 | 0.034 | 0.006 0.145 |
本で示されていたものと同じ結果を得ることが出来ました。詳しい解説は省略します← それでは、求めたモデルを使ってyを予測し、元の分布と比較してみます。
x = np.linspace(6, 13, 1000)
result = model.predict({'x': x})
plt.plot(x, result)
Ts = d[d['f']=='T']
Cs = d[d['f']=='C']
plt.xlabel('x', size=20)
plt.ylabel('y', size=20)
T = plt.scatter(Ts['x'], Ts['y'], color='b', label='T')
C = plt.scatter(Cs['x'], Cs['y'], color='r', label='C')
plt.legend((T, C), ('T', 'C'))
<matplotlib.legend.Legend at 0x10b5e2390>
とりあえず本と同じようなグラフが書けました(笑)
今回行った演習は以上です。
このモデルの場合、どの個体も同じ特性を持っていると仮定してモデル化されています。しかし、今回扱った種子のデータでそのようなことはありえません。必ず種子には個体差が含まれます。例えば、種子の個体差以外にも植物が育てられた土、与えた水の量、日射量などに違いが存在します。したがって、これらを無視して「どの種子も同じ分布に従う」という仮定の元モデル化してもモデルの正しさには疑問が残ります。 第3章以降ではこのようなばらつきを考慮したモデルについても紹介されています。
Pythonでも同じことができる。ただし、いろんなモジュールをインストールする必要があるので少し面倒。
最後のグラフを書いている最中、この記事の上位互換なるものを発見してしまったのでシェアさせていただきます...
緑本に載っていないような検証もPythonでされているようなので非常に参考になると思います(今後参考にさせていただきます!)