いろいろするblog

プログラミングとか統計とかえーあい的なあれこれとか

データ解析のための統計モデリング入門をpythonでやる（3）-一般化線形モデル-

【あらすじ】
種子データを一般化線形モデルであらわしてみた。

【本文】
説明変数を組み込んだ統計モデルを扱う。
つまり、個々のデータが独自の説明変数の値を持ち、それによって従うモデルが変わるようになる。
このようなモデルを一般化線形モデルという。

今回はある植物のデータとして
種子数:y_i
体サイズ:x_i
施肥処理:C or T
であるデータを扱う。

前回のポアソン分布

\frac{\lambda^{k} exp(-\lambda)}{k!}

において、\lambdaが体サイズx_iによって変わることを考える。
今回は以下の式に従うと考える。

\lambda_i = exp(\beta_1 + \beta_x x_i)

ここで、logを取ると

\log \lambda_i = \beta_1 + \beta_x x_i

となる。
右辺のように\betaの線形結合になるものを線形予測子と呼ぶ。
このように
(\lambda_i の関数) = (線形予測子)
となる場合、左辺の\lambda_i の関数をリンク関数と呼ぶ。
今回は対数を取っているので対数リンク関数である。

平均がこの場合の対数尤度を求める。
つまり、

\log_10 L(\beta_1 , \beta_2) = \sum{(y_i \log \lambda - \lambda - \sum_{k}^{y_i} \log k)}

が最大となる \beta_1 , \beta_2 を求める。
pythonでざっくり当てはめるにはstatsmodelを使う。

[python のコード]

model = sm.formula.api.glm(formula="y ~ x", data=data, family=sm.api.families.Poisson(link=sm.api.families.links.log))

formula="y ~ x" は yの説明変数はx であることを示している。
もしy の説明変数をfと考える場合（つまり、体サイズx は関係なく、施肥の有無によってのみ種子数がかわる）は"y ~ f"とする。
また、y の説明変数をxとfの両方と考える場合（体サイズと施肥の両方が関係する）は"y ~ x + f"とする

yに対する説明変数がxのみの場合、一般化線形モデルではxの値の位置に依存した異なる分布でyの値が表現される。

データの特徴を反映した確率分布とリンク関数を使いましょう。

データ解析のための統計モデリング入門をpythonでやる（2-2）-poisson分布-

【あらすじ】
ポアソン分布の性質を確認した。
最尤推定を試した。

【本文】
ポアソン分布についてみてみる。

ポアソン分布は以下の確率質量関数に従う。
ポアソン分布は離散分布であり、パラメータは平均だけである。
分散は平均と等しい。

P(X = k) = \frac{\lambda^{k} exp(-\lambda)}{k!}

\lambda : 平均

numpyはnumpy.random.poissonでポアソン分布に従う整数の配列が作れる。
scipyはscipy.stats.poissonでポアソン分布が扱える。

{ここにpython コードを埋め込む}

ある観測データとモデルの間には尤度（あてはまりの良さ）なる統計量を定義できる。
尤度はある観測データの集合が得られたときの、各観測データが得られる確率の積である。
尤度はLであらわされる。
あるパラメータ\lambdaの時の尤度はL(\lambda)である。
平均\lambdaのポアソン分布と観測データの集合\{ y_i \}の場合の尤度を式で表すと以下のようになる。

L(\lambda) = \prod{p(y_i | \lambda)} = \prod{\frac{\lambda^{y_i} exp(-\lambda)}{y_i!}}

これだと値が小さくなりすぎて使いにくいので対数を取った対数尤度(\log_{10} L)を使う。

\log_{10} L = \sum{(y_i \log \lambda - \lambda - \sum_{k}^{y_i} \log k)}

この対数尤度(\log_{10} L)が最大となる\lambdaが最もあてはまりの良いモデルとなる。

対数尤度(\log_{10} L)が最大となる\lambdaを求めるには微分して傾きゼロの値を求めればよいので

\frac{\partial \log L(\lambda)}{\partial \lambda} = \sum \{ \frac{y_i}{\lambda} - 1\}

これがゼロの時は

\sum \{ \frac{y_i}{\lambda} - 1\} = 0
\frac{1}{\lambda} \sum y_i - (総データ数) = 0
\lambda = \frac{\sum y_i}{(総データ数)}

となり、ポアソン分布においては標本平均と等しくなる。

このようにしてパラメータを決める手法を最尤推定という。

データ解析のための統計モデリング入門をpythonでやる（2-1）-pyper導入-

【あらすじ】
pythonで.RData形式のファイルを読むためにpyperを使えるようにした。
pythonで.RData形式のファイルを読んだ。

【内容】
1)pyper導入
データ解析のための統計モデリング入門で使われるデータは著者が公開しており、以下のページにあります。
生態学データ解析 - 本/データ解析のための統計モデリング入門

これらのデータはRで使われるデータ形式である.RDataなので、そのままpythonで読めません。
というわけで、pythonからpyperを使ってRを操作して読んでいきます。

pip install pyper

pyperはpythonからRに命令を飛ばす橋渡しをするだけなので、Rをインストールしていないとpyperをimportしても使えません。
というわけでRを導入します。

普通にインストールする場合は以下のようなページを参考にしていけばいけるんじゃないですかね。
qiita.com

chocolateyを使う場合は

cinst R

でいけます。
いけますが、この方法だとPATHが通ってないっぽいのでpyper.R()をするときに呼び出し先がいないと怒られます。
その時は

pyper.R(RCMD="(Rのexeが置いてあるフォルダのパス)")

で指定してあげるとちゃんと動きます。

2).RData形式のデータを読む。
pyperで読んだ後、とりあえずpandasのSeriesに入れます。

import pyper
import pandas as pd

Rを触るためのオブジェクトを作ります。
また、このとき後でpandasとnumpyでも使えるようにしときます。
r = pyper.R(RCMD=r"(Rのパス)", use_numpy=True, use_pandas=True)

rに.RDataを読ませます。
今回はカレントディレクトリにあるdata.RDataを読ませます。
r('load("data.RData")')

data.RDataに入っている"data"という名前のデータを見たいときは
r.get("data")

とすると配列になって出てきます。
これをpandasのSeriesに入れます。
data = pd.Series(r.get("data"))

2章の一番最初のデータの要約。
In [47]: data.describe()
Out[47]:
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

histでヒストグラムの表示。
data.hist()

f:id:yshim7:20181119110759p:plain

見づらいのでビン範囲と幅の調整とグリッドを非表示に。
data.hist(bins=[x / 10 for x in range(5, 95, 10)], grid=False, rwidth=0.8)

f:id:yshim7:20181119110819p:plain

満足。

データ解析のための統計モデリング入門をpythonでやる（1）

【あらすじ】
これからデータ解析のための統計モデリング入門をpython で実装しながら読んでいく。

【本文】
データ解析のための統計モデリング入門を読みます。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

この本ではRを使っていますが、pythonの方が慣れてるのでpython(3.7.0)を使います。
あと、実行はipython（7.1.1）上でやっていきます。
他にはnumpy(1.15.4)とかpadas(0.23.4)とか使う予定。
あと、Rのデータを読むのにR(3.5.1)とpyper(1.1.2)。

その他の環境として、win10上でCmderを使ってターミナル的な環境を作っています。
python以外のあれこれの導入はChocolateyを使用。