【Python】緑本こと『データ解析のための統計モデリング入門』を実装していく【第2章その2】

間隔が空いてしまいましたが、緑本の第2章の続きです。

ちなみに、手元では最終章以外は既に実装まで終わっています（最後の最後に出てくる条件付き自己回帰の話が理解できずに詰んでます…）。

植物の種子数の統計モデリングということで、前回は手元にあるデータを分析し、

データはカウントデータ（非負の整数）
種子数の標本平均は3.56、分散2.99
個体ごとの種子数のヒストグラムはひと山の分布で表せそう

ということが分かったのでした。今回はそれを踏まえてモデリングしていきます。
ただ写経しても意味がないので、実装中に自分の頭の中で考えていることも書き下していきます（間違いもあるかと思うので、是非ご指摘お願い致します…）。

データがどのような確率分布から発生したか考える
観測データに基づいて、仮定した確率分布のパラメータを最尤推定（標本に基づく母集団分布のパラメータ推定）

データがどのような確率分布から発生したか考える

種子数データが従う真の確率分布（母集団分布）がどうなっているかを考える（真の分布を実際に見ることはできない）。

確率分布を使うことで、任意の個体の種子数がN(N=0,1...)個となる確率の一覧が表せる。これによって、個体ごとの種子数のばらつきを比較的簡単な数式(確率密度・質量関数)で表現できる。

今回はデータについて、以下のことが分かっている。

データが非負の整数
下限は0だが、上限はよくわからない（7ぐらい？でも7以上の可能性→二項分布は使えない）
標本の平均が3.56、分散2.99でわりと近い値

上記の特徴からポアソン分布が使えそうかな？と考える（ここは知識と経験が必要そう）。

→ここでは、全ての個体の種子数が平均 $\lambda$ (未知)の同一のポアソン分布に従っていると仮定。

観測データに基づいて、仮定した確率分布のパラメータを最尤推定（標本に基づく母集団分布のパラメータ推定）

尤度(関数) $L(\lambda)$ を導入。尤度の実態は、パラメータがある値のとき（パラメータの真の値は不明だから、とりあえず任意の値を仮定する）、それぞれの標本について確率関数の値の総積をとったもの。尤度は「手元にある観測データへの当てはまりのよさ」だと解釈できる。

では尤度 $L(\lambda)$ を最大にする $\lambda$ をどうやって決めるか？ $\lambda$ に片っ端から実数を当てはめていたら日が暮れてしまう。

→尤度関数（実際には対数変換した対数尤度関数）は上に凸だから、微分して傾きが0となる $\lambda$ を求める（これを最尤推定と呼ぶ）。解析的に $\frac{\partial logL(\lambda)}{\partial \lambda}=0$ を解くと $\hat{\lambda}=3.56$ が求まる。

（ちなみに、実装の際には負の対数尤度最小化問題と考えて以下のようにしても解を求められる。）

from scipy import optimize
from scipy.stats import poisson
def f(lambda_):
return -sum(poisson.logpmf(df.x, lambda_))
optimize.minimize(f, 1.0)

          fun: 97.24400294080678
hess_inv: array([[0.07117249]])
jac: array([-5.7220459e-06])
message: 'Optimization terminated successfully.'
nfev: 27
nit: 8
njev: 9
status: 0
success: True
x: array([3.55999988])

→3.56で尤度が最大化（負の尤度最小化）されると分かる。

$\lambda=3.56$のポアソン分布と実際の観測データの分布の比較

# 平均3.56のポアソン分布（青）と観測データの分布（グレー）を重ねる
# 参考：スケールの異なる2つのグラフを重ねて描画する方法（https://qiita.com/terakawa@github/items/24f31ddca644a7d4bb49）
fig, ax1 = plt.subplots()
plt.hist(df.x, bins=10, color=['gray'])
poisson_values = poisson.rvs(3.56, size=10000)
prod = (pd.value_counts(poisson_values) / 10000).sort_index()
ax2 = ax1.twinx()
ax2.plot(prod)
plt.show()

f:id:mimosom:20190512110809p:plain

上図を見る限り、うまく当てはまってそう。従って、今回の架空植物の種子数は $\lambda=3.56$ のポアソン分布に従っていると言えそう。

（本章では推定したモデルによる予測については議論しない、3章以降に議論）

余談：λの値を変化させていったポアソン分布と観測データ

f:id:mimosom:20190430192828p:plain

この中だと $\lambda=3.5$ のポアソン分布が観測データに一番フィットしてそう。従って、上で求めた推定値( $\lambda=3.56$ )は外れてなさそう。

今回の要点（まとめ）

話の流れとしては、

観測データをよく見て、そのデータを生み出した真の確率分布（パラメータ未知）を仮定する。
そうすると真の分布のパラメータを変数とする尤度関数が求まるので、それを最大化するようなパラメータを何かしらの手段で求める。
求められたパラメータの値（最尤推定値）を1.で仮定した確率分布に代入すると、その観測データを生み出した真の分布がどういう形になっているか推定できる。

というものでした。

以上で、第2章はおしまいです。お疲れ様でした。

データ解析のための統計モデリング入門一般化線形モデル・階層ベイズモデル・MCMC

現象を数理モデルで表現・説明するのに慣れていない人のために、章ごとに異なる例題を解決していく過程を通して、統計モデルの基本となる考えかたを説明する。前半では、応用範囲のひろい統計モデルのひとつである一般化線形モデルの基礎を、後半では、実際のデータ解析に使えるように、それらをベイズ統計モデル化する方法を説明する。

Amazon