確率密度関数の推定と結果の確認

# 概要

## 環境
- python3.5
- numpy==1.13.3
- scikit-learn==0.19.1
- matplotlib==2.1.0
- scipy==1.0.0


## やりたいこと

1次元、2次元のデータ分布に対して、確率密度推定を行い、その結果を以下のように図示して確認すること。

||データ分布|確率密度推定結果|
|---|---|---|
|1D|![image.png](https://qiita-image-store.s3.amazonaws.com/0/104992/02736516-4c0d-2cca-d1ea-78fe827ada06.png)|![image.png](https://qiita-image-store.s3.amazonaws.com/0/104992/a6a8f7b2-f31c-9eba-d108-0ce88c8e7104.png)|
|2D|![image.png](https://qiita-image-store.s3.amazonaws.com/0/104992/6a20f4bc-a952-2abf-28c2-09ad9c8dbe0c.png)|![image.png](https://qiita-image-store.s3.amazonaws.com/0/104992/08e408bd-c1ed-6767-1555-a4e9a551767b.png)|

## 確率密度推定のモデル
データ標本から確率密度関数を推定する方法は、一般的に2つに大別できる。

1. パラメトリックモデル: 母集団が正規分布やガンマ分布などであることを想定するモデル
2. ノンパラメトリックモデル: 母集団の分布を想定しないモデル。

## 手法

ここでは、ノンパラメトリックの代表的な手法と言えるKDE(Kernel Density Estimation)を用いて、標本の確率密度関数を推定する。

KDEは、標本に対して複数のカーネルを適用し、それらのカーネルの重ね合わせでデータ標本の確率密度関数を推定する方法である。

## KDEライブラリ
Python環境では、KDEを実装している主要なライブラリはいくつかある。

- scikit-learn
- scipy
- statsmodels

今回はscikit-learnを用いる。

# 1次元の場合

正規分布に従うランダムな値を生成する。
ここでは、2つの異なる正規分布を使っておく。

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats.distributions import norm

np.random.seed(0)
x = np.concatenate([norm(-1, 1.).rvs(200), norm(1, 0.3).rvs(100)])

# ヒストグラムで確認する
fig, ax = plt.subplots(figsize=(17, 4))
ax.hist(x, 10, fc='gray', histtype='stepfilled', alpha=0.5, normed=True)
ax.set_xlim(x.min()-1, x.max()+1)
plt.show()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/104992/02736516-4c0d-2cca-d1ea-78fe827ada06.png)

KDEを実行する。ここでは、scikit-learnにあるKDEを使う。


```python
from sklearn.neighbors.kde import KernelDensity
X = x.reshape(-1, 1)
kde = KernelDensity(kernel='gaussian', bandwidth=0.2).fit(X)
kde.score_samples(X)
```

なお、`kde.score_samplex()`が返す値は、log-likelihood（対数尤度）であり、対数となっている。expを撮って、積分すると（ほぼ）1である。
参考までに、確認しておく。

```python
from scipy import integrate
chk = integrate.quad(lambda x: np.exp(kde.score_samples(x)), -np.inf, np.inf)
assert np.allclose(chk, [1.0, 0.0], rtol=1e-05, atol=1e-07)
```

生成したランダムな値のヒストグラムとKDEの曲線を同時に図示して確認する。

```python
x_plot = np.linspace(-5, 5, 50)[:, np.newaxis]
fig, ax = plt.subplots(figsize=(17, 4))
ax.plot(x_plot, np.exp(kde.score_samples(x_plot)), linewidth=3, alpha=0.5)
ax.hist(X, 10, fc='gray', histtype='stepfilled', alpha=0.5, normed=True)
ax.set_xlim(X.min()-1, X.max()+1)
plt.show()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/104992/a6a8f7b2-f31c-9eba-d108-0ce88c8e7104.png)


# 2次元の場合

正規分布に従うランダムな値を2つ生成する。
ここでも、2つの異なる正規分布を使っておく。

```python
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats.distributions import norm

np.random.seed(0)
x = np.concatenate([norm(-1, 1.).rvs(200), norm(1, 0.3).rvs(100)])
y = np.concatenate([norm(-1, 1.).rvs(200), norm(1, 0.3).rvs(100)])

# 2Dヒストグラムで確認する
import matplotlib.cm as cm
fig = plt.figure(figsize=(8, 6))
ax = fig.gca()
ax.set_xlim(x.min(), x.max())
ax.set_ylim(y.min(), y.max())

H = ax.hist2d(x, y, bins=60, normed=True, cmap=cm.gray)
fig.colorbar(H[3],ax=ax)

plt.show()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/104992/6a20f4bc-a952-2abf-28c2-09ad9c8dbe0c.png)


KDEを実行する。

```python
from sklearn.neighbors.kde import KernelDensity
X = np.vstack((x, y)).T

kde = KernelDensity(kernel='gaussian', bandwidth=0.2).fit(X)
np.exp(kde.score_samples(X))
```

ここでも、expをとって、積分すると（ほぼ）1であることを確認しておく。少しだけ時間がかかります。


```python
from scipy import integrate
chk = integrate.dblquad(lambda _x, _y: np.exp(kde.score_samples(np.array([[_x, _y]]))), -np.inf, np.inf, gfun=lambda x: -np.inf, hfun=lambda x: np.inf)
assert np.allclose(chk, [1.0, 0.0], rtol=1e-05, atol=1e-07)
```

生成したランダムな値のヒストグラムとKDEの曲線を同時に図示して確認する。

```python
f0 = np.arange(x.min(), x.max(), 0.1)
f1 = np.arange(y.min(), y.max(), 0.1)
xx, yy = np.meshgrid(f0, f1)

positions = np.vstack([xx.ravel(), yy.ravel()]).T

scores = np.exp(kde.score_samples(positions)).T

f = np.reshape(scores, xx.shape)

fig = plt.figure(figsize=(8, 6))
ax = fig.gca()
ax.set_xlim(x.min(), x.max())
ax.set_ylim(y.min(), y.max())

cset = ax.contour(xx, yy, f, colors='r')
ax.clabel(cset, inline=1, fontsize=10)

H = ax.hist2d(x, y, bins=60, normed=True, cmap=cm.gray)
fig.colorbar(H[3],ax=ax)

plt.show()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/104992/08e408bd-c1ed-6767-1555-a4e9a551767b.png)


以上。