データ解析のための統計モデリング入門


自分の勉強のため．
中はRで書かれてるのでそれをPythonに書き換えて頑張る．
コード中のコメントアウトはR言語．

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

Rには種子数のカウントデータがあるらしい．
多分numpy,pandasを使って処理した方がいい気がするから，
`numpy.array`と`pandas.Series`それぞれでもデータを用意．
あとグラフ用に`pyplot`も．

```pycon
>>> data = [2, 2, 4, 6, 4, 5, 2, 3, 1, 2, 0, 4, 3, 3, 3, 3, 4, 2, 7, 2, 4, 3, 3, 3, 4, 3, 7, 5, 3, 1, 7, 6, 4, 6, 5, 2, 4, 7, 2, 2, 6, 2, 4, 5, 4, 5, 1, 3, 2, 3]
>>> import numpy
>>> import pandas
>>> matplotlib.pyplot as plt
>>> ndata = numpy.asarray(data)
>>> pdata = pandas.Series(data)
```

種子数は50個．

```pycon
# length(data)
>>> len(data)
50
>>> len(ndata)
50
>>> len(pdata)
50
```


データの標本平均や最小値・最大値・四分位数などを表示するには

```pycon
# summary(data)
>>> pdata.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
```

度数分布を取得する

```pycon
# table(data)
>>> pdata.value_counts()
3    12
2    11
4    10
5     5
7     4
6     4
1     3
0     1
dtype: int64
```

種子数５は５個体，種子数６は４個体であることを確認できる．

これをヒストグラムで表示する

```pycon
# hist(data, breaks = seq(-0.5, 9.5, 1))
>>> pdata.hist(bins=[i - 0.5 for i in xrange(11)])
<matplotlib.axes._subplots.AxesSubplot object at 0x10f1a9a10>
>>> plt.show()
```

![Screen Shot 2015-05-29 at 23.38.08.png](https://qiita-image-store.s3.amazonaws.com/0/48856/6ecaafa9-209d-a3f4-ece6-3a6355500e6a.png)

データのばらつきを表す統計量”標本分散”の算出．
numpyとpandasでデフォルトの分散の計算方法が違うらしい．（[参考](http://sinhrks.hatenablog.com/entry/2014/11/27/232150)）
ちなみにRの場合は不偏推定量っぽい．


```pycon
# var(data)
>>> numpy.var(ndata) # 標本統計量
2.9264000000000006
>>> numpy.var(ndata, ddof=True) # 不偏推定量 
2.9861224489795921
>>> pdata.var() # 不偏推定量
2.986122448979593
>>> pdata.var(ddof=False) # 標本統計量
2.926400000000001
```


標準偏差のほうも分散と同様の挙動を示す模様．
今回はpandasの例だけ．


```pycon
# sd(data)
>>> pdata.std()
1.7280400600042793

# sqrt(var(data))

```

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

種子数データには以下の特徴があることがわかる．

- １個，２個，...と数えられるカウントデータ
- １個体の種子数の標本平均は$3.56$
- 個体ごとの趣旨ううにはばらつきがあり，ヒストグラムを描くと一山の分布になる

ばらつきを表すため，**確率分布**を使う．
種子数データを統計モデルとして表現するために，
今回は**ポアソン分布**を用いるらしい．

ある植物個体$i$の種子数$y_i$を**確率変数**と呼ぶ．

> 個体１の種子数$y_1=2$となりうる確率はどれくらいか？

を表現する確率分布は比較的簡単な数式で定義され，パラメータによってその形が決定する．

ポアソン分布の場合は，パラメータは**分布の平均**のみ．

今回の例題で言えば，平均$3.56$のポアソン分布となる．
分布関係は`scipy`がいいっぽい．

```pycon
# y <- 0:9
# prob <- dpois(y, lambda = 3.56)
# plot(y, prob, type="b", lyt=2)
>>> import scipy.stats as sct
>>> y = range(10)
>>> prob = sct.poisson.pmf(y, mu=3.56)
>>> plt.plot(y, prob, 'bo-')
>>> plt.show()
```

![Screen Shot 2015-05-30 at 00.34.21.png](https://qiita-image-store.s3.amazonaws.com/0/48856/3aa5873a-1051-d31f-6602-4e82cc759ef7.png)

実際のデータと重ね合わせてみる

```pycon
>>> pdata.hist(bins=[i - 0.5 for i in range(11)])
>>> pandas.Series(prob*50).plot()
>>> plt.plot()
```

![Screen Shot 2015-05-30 at 00.45.01.png](https://qiita-image-store.s3.amazonaws.com/0/48856/e4ecd2a0-40d5-d916-c135-64be758ec022.png)

この結果ならポアソン分布で観察データを表現できていると考えられる．


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

ポアソン分布についてもう少し詳しく．

### ポアソン分布の定義

平均が$\lambda$であるときの確率変数が$y$であるときの確率

```math
p(y|\lambda) = \frac{\lambda ^{y}\exp(-\lambda)}{y!}
```

### 性質

- $\lambda$の値によって曲線の形が変わる
- $y \in${$0, 1, 2, \dots,\infty$}の値をとり，全ての$y$について

```math
\sum^{\infty}_{y=0}p(y|\lambda) = 1
```

- 確率分布の平均は$\lambda$である（$\lambda \geq 0$）
- 分散と平均は等しい


### ポアソン分布を今回選んだ理由

1. データに含まれている値$y_i$が非負の整数である（カウントデータである）
2. $y_i$に下限（０）はあるが，上限はない
3. 観測データの平均と分散がほぼ等しい


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

### 最尤推定法

**尤度**と呼ばれるモデルの「当てはまりの良さ」を表す統計量を決める手法．

ポアソン分布の場合は$\lambda$を決める．

> ある$\lambda$が決まった場合の全ての個体$i$についての確率$p(y_i|\lambda)$の積

尤度は$L(\lambda)$と表記される．
今回の場合は，

```math
\begin{eqnarray*}
L( \lambda ) &=& (y_1が2である確率) \times (y_2が2である確率) \times \dots \times (y_50が3である確率)\\
&=& p(y_1 | \lambda) \times p(y_2 | \lambda) \times p(y_3 | \lambda) \times \dots \times p(y_50 | \lambda)\\
&=& \prod_{i}p(y_i | \lambda) = \prod_{i} \frac{\lambda^{y_i} \exp (-\lambda)}{y_i !}
\end{eqnarray*}
 

```

となる．50個の個体が観測と同じになる確率（＝50個の事象が真である確率）を計算するため，積にする．

尤度関数$L(\lambda)$のままだと使いにくいため，大体は**対数尤度関数**（**log likelihood function**）を使う．

```math
\begin{eqnarray*}
\log L(\lambda) &=& \log \left( \prod_{i} \frac{\lambda^{y_i} \exp (-\lambda)}{y_i !} \right) \\
&=& \sum_i \log \left( \frac{\lambda^{y_i} \exp (-\lambda)}{y_i !} \right) \\
&=& \sum_i \left( \log(\lambda^{y_i} \exp (-\lambda)) - \log y_i !   \right)\\
&=& \sum_i \left( y_i \log \lambda - \lambda - \sum^{y_i}_{k} \log k \right)
\end{eqnarray*}
```

この対数尤度$\log L(\lambda)$と$\lambda$の関係をグラフにする．

```pycon
# logL <- function(m) sum(dpois(data, m, log=TRUE))
# lambda <- seq(2, 5, 0.1)
# plot(lambda, sapply(lambda, logL), type="l")
>>> logL = lambda m: sum(sct.poisson.logpmf(data, m))
>>> lmd = [i / 10. for i in range(20, 50)]
>>> plt.plot(lmd, [logL(m) for m in lmd])
>>> plt.show()

```

![Screen Shot 2015-06-01 at 01.34.53.png](https://qiita-image-store.s3.amazonaws.com/0/48856/a5bbea2d-4766-c561-edbe-91fddeeabb30.png)

対数尤度$\log L(\lambda)$は尤度$L(\lambda)$の単調増加関数であるから，
対数尤度が最大のとき，尤度も最大となる．

> グラフを見た感じ最も確率が高いのは，$\lambda = 3.5$付近ということがわかる．

具体的な最大値は対数尤度の極大値を求めれば良い（＝傾きが０になるとき）

つまり，

```math
\begin{eqnarray*}
\frac{\partial \log L(\lambda)}{\partial \lambda} = \sum_i \left\{ \frac{y_i}{\lambda} - 1 \right\} = \frac{1}{\lambda} \sum_i y_i - 50 &=& 0 \\
\lambda &=& \frac{1}{50}\sum_i y_i \\
&=& \frac{全部のy_iの和}{データ数}(=データの標本平均) \\
&=& 3.56
\end{eqnarray*}
```

となり，対数尤度および尤度が最大になる$lambda$を**最尤推定量**，具体的な$y_i$で評価された$\lambda$を**最尤推定値**と呼ぶ．


### 一般化

$\theta$をパラメータとする確率分布から観測データ$y_i$が発生した場合の確率を$p(y_i | \theta)$とすると，
尤度，対数尤度は

```math
\begin{eqnarray*}
L(\theta | Y) &=& \prod_i p(y_i | \theta) \\

\log L(\theta | Y) &=& \sum_i \log p(y_i | \theta)

\end{eqnarray*}
```

となる．

# 確率分布の選びかた

次の点について考える．

- 説明したい量は離散か連続か？
- 説明したい量の範囲は？ 
- 説明したい量の標本平均と標本分散の関係は？

カウントデータの統計モデルで用いる分布の例として：

- **ポアソン分布**：データが離散値，ゼロ以上の範囲，上限特になし， 平均$\approx$分散
- **二項分布**：データが離散値，ゼロ以上で有限の範囲$\{0, 1, 2, \dots, N\}$，分散は平均の関数

また連続分布については

- **正規分布**：データが連続値，範囲が$[-\infty, +\infty]$，分散は平均とは無関係に決まる
- **ガンマ分布**：データが連続値，範囲が$[0, +\infty]$，分散は平均の関数
- **一様分布**：データが連続値，有界

次は[一般化線形モデル（GLM）](http://qiita.com/makora9143/items/851b69a7f7af53d16e5f)．