他にも多くの同様のブログポストがあるので何番煎じか分かりませんが，久保先生の緑本 を Python で書く自主練を始めることにしました。第2章から始めます。

GitHub リポジトリの URL は https://github.com/tnoda/midorihon-in-python です。

第2章 確率分布と統計モデルの最尤推定

Pandas 使います。

%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
from numpy.random import randn
import numpy as np
from pandas import Series, DataFrame
import pandas as pd
pd.__version__

u'0.17.0'

.RData データの読み込み

例題データの内容を Python で読み込みます。 .RData を Python で扱うには rpy2 が便利です。

from rpy2.robjects import r

r['load']('data.RData')
r['data']

<floatvector python:0x11eab8128 r:0x11ef42ea0>
[2.000000, 2.000000, 4.000000, ..., 3.000000, 2.000000, 3.000000]

2.1 例題: 種子数の統計モデリング

今回は R のデータを Pandas の Series に変換して扱います。

data = Series(np.array(r['data']))
data

0     2
1     2
2     4
3     6
4     4
5     5
6     2
7     3
8     1
9     2
10    0
11    4
12    3
13    3
14    3
15    3
16    4
17    2
18    7
19    2
20    4
21    3
22    3
23    3
24    4
25    3
26    7
27    5
28    3
29    1
30    7
31    6
32    4
33    6
34    5
35    2
36    4
37    7
38    2
39    2
40    6
41    2
42    4
43    5
44    4
45    5
46    1
47    3
48    2
49    3
dtype: float64

この data に含まれるデータの数を確認するには len 関数を使います。

len(data)

50

Series.describe()

また describe() メソッドによって、標本平均や最小値・最大値・四分位数などが分かります。

data.describe()

count    50.00000
mean      3.56000
std       1.72804
min       0.00000
25%       2.00000
50%       3.00000
75%       4.75000
max       7.00000
dtype: float64

上の describe() の読みかたを少し説明してみましょう。 min と max はそれぞれ data 中の最小値・最大値です。また 25%, 50%, 75% は、それぞれ data を小さい順に並べたときの 25%, 50%, 75% 点の値です。出力された観測点そのものではなく、たいていは補間予測された値です。 50% は標本中央値（中位値）とよばれる推定値、そして、 mean は標本平均 (sample mean) でいまの場合には 3.56 です。

Series.value_counts()

Pandas で度数分布を得る方法はいろいろありますが、ここでは value_counts() を使ってみましょう。

data.value_counts()

3    12
2    11
4    10
5     5
7     4
6     4
1     3
0     1
dtype: int64

これを見ると、一番多いのは種子数 3 で 12 個体、二番目に多いのは種子数 2 で 11 個体というのが分かるのですが、個体数でソートされているので見慣れた度数分布の表にはなっていません。個体数でのソートを避けるには sort=False オプションを value_counts() に追加します。

data.value_counts(sort=False)

0     1
1     3
2    11
3    12
4    10
5     5
6     4
7     4
dtype: int64

Series.hist()

それをヒストグラム (histogram) として図示してみましょう

plt.hist(data, bins=np.arange(-0.5, 8.5, 1.0))

ここで、

np.arange(-0.5, 8.5, 1.0)

は、 ({-0.5, 0.5, 1.5, …, 8.5}) という数列を生成します。

Series.var()

あるデータのばらつき (variability) をあらわす標本統計量の例として、標本分散 (sample variance) があげられます。

data.var()

2.986122448979592

Series.sd()

また、標本標準偏差 (sample standard deviation) とは標本分散の平方根です。 Python ではこんなふうに計算できます。

data.std()

1.728040060004279

from math import sqrt

sqrt(data.var())

1.728040060004279

2.2 データと確率分布の対応関係をながめる

poisson.pmf() (probability mass function)

「平均 3.56 のポアソン分布」なるものを， Python を使ってグラフとして図示しましょう。平均 3.56 のポアソン分布にしたがって「種子数が (y) であると観察される確率」を生成させるためには，たとえば以下のように Python に指示します。

from scipy.stats import poisson

y = np.arange(0, 10)
prob = Series(poisson.pmf(y, 3.56), index=y)
prob

0    0.028439
1    0.101242
2    0.180211
3    0.213851
4    0.190327
5    0.135513
6    0.080404
7    0.040891
8    0.018197
9    0.007198
dtype: float64

plot(), 'o–’

これでは分かりにくいので図示してみると，

plt.plot(prob, 'o--')

このように plot() 関数を呼ぶと y と prob との関係が図示されます。ここで， plot() 関数の引数 o-- によって「丸と折れ線による図示」になります。

つぎに，観測データのヒストグラムに平均 3.56 のポアソン分布を重ねてみます。

plt.plot(prob.mul(50), 'o--')
plt.hist(data, bins=np.arange(-0.5, 8.5, 1.0))

観察されたばらつきがポアソン分布で表現できているみたいだなぁと思いました。

2.3 ポアソン分布とは何か？

さまざまな平均 (($\lambda$)) のポアソン分布を描いてみます。

y = np.arange(0, 21)
lambdas = [3.5, 7.7, 15.1]
styles = ['ro--', 'go--', 'bo--']
labels = ['3.5', '7.7', '15.1']
for i in range(3):
    prob = Series(poisson.pmf(y, lambdas[i]), index=y)
    plt.plot(prob, styles[i], label=labels[i])
plt.legend(loc='best', title='lambda')
plt.xlabel('y')
plt.ylabel('prob')

2.4 ポアソン分布のパラメータの最尤推定

subplots()

平均 ($\lambda$) (lambda) を変化させていったポアソン分布と，観測データのあてはまりの良さ（対数尤度 logL ）をグラフにしてみましょう。

y = np.arange(0, 10)
lambdas = np.arange(2.0, 5.6, 0.4)

# Figure
fig, axes = plt.subplots(3, 3, sharex=True, sharey=True, figsize=(12, 8))
for i in range(9):
    mu = lambdas[i]
    prob = Series(poisson.pmf(y, mu), index=y)
    logL = sum(poisson.logpmf(data, mu))
    row, col = i / 3, i % 3
    ax = axes[row, col]
    ax.plot(prob.mul(50), 'o--', label='lambda=%.1f' % mu)
    ax.hist(data, bins=np.arange(-0.5, 8.5, 1.0), label='log L=%.1f' % logL)
    ax.legend(loc='best', fontsize=8)

FontProperties

この図を見ると，対数尤度が大きい（ゼロに近い）ほど観測データとポアソン分布が「似ている」ように見えます。さらに対数尤度 ($ log L(\lambda) $) と ($ \lambda $) の関係を調べるために，次のような図をつくってみましょう。

from subprocess import check_output

lambdas = Series(np.arange(2, 5, 0.1))
logL = lambda m: sum(poisson.logpmf(data, m))
plt.plot(lambdas, lambdas.apply(logL), '-')

2.4.1 擬似乱数と最尤推定値のばらつき

この例題の架空データはポアソン分布にしたがう乱数列の一つです。乱数発生関数が生成する乱数列は毎回異なるので，50 個のポアソン乱数の標本平均も試行ごとに異なります。乱数列を発生させるごとに，最尤推定値 ($ \hat{\lambda} $) はどのように変わるのでしょうか？ これを調べるために Python を使った乱数実験をやってみましょう。例題データを作ったときと同じように乱数発生関数でデータを生成し，そのたびに標本平均（最尤推定値）を評価して記録します。これを 3000 回繰り返した結果を次に示します。データが 50 個体ぶんしかないため，それなりのばらつきがあるでしょう。このような推定値のばらつきは標準誤差と呼ばれ，その大きさは調査個体数に依存しています。

res = []
for i in range(3000):
    res.append(np.random.poisson(3.5, 50).mean())
plt.hist(res, bins=20)

# LaTeX
from matplotlib import rc
rc('text', usetex=True)
plt.xlabel(u'\hat{\lambda}')

平均値が 2.5 であるポアソン乱数を発生させるには Python で

np.random.poisson(2.5, 50)

array([4, 5, 5, 3, 0, 6, 4, 2, 3, 5, 0, 7, 3, 3, 1, 1, 1, 1, 5, 4, 2, 1, 1,
       4, 2, 1, 5, 4, 2, 3, 0, 1, 3, 3, 1, 2, 4, 0, 2, 1, 4, 4, 4, 0, 1, 2,
       6, 3, 3, 2])

と指定します。