公開！週末研究ノート０６　ー　任意の分布からのサンプリング

# はじめに　ー　週末研究ノートとは？

個人的に研究的な活動をやるにあたり、オープンにしてみたら面白いかもと思い、自分が興味を持っている ざっくりテーマについて、これから、ゆるい週末研究を公開していこうと思います。（有識者の方のアドバイスも、ちょっとというかかなり期待してます！笑）

どこかの権威的な学会やジャーナルなどで発表する予定はないため、万が一、私の記事を利用する際には自己責任でお願いします。そんな人はいないと思いますが、念のため。



# 今回のサマリ （TL; DR）

任意の分布関数からサンプリングする方法を実際に試しました。


# 環境

- Colaboratory
    - コードは、[こちら](https://gist.github.com/tkosht/73cd4bf45f102fa1986a2f96c1f02eab)


# 今回の週末研究ノート

今回は、任意の分布からサンプリングってできるはずだけど、どうやるんだっけ？を調査し実機で確認しました。
よく確率密度関数とか頑張って計算してるのを見るけど、そこからサンプリングできるイメージがわかないと、何やってるんだっけってなる。そう、迷子になります笑
がんばって計算してるけど・・・で？ってなる笑

## 棄却サンプリング

今回は、比較的シンプルな棄却サンプリングを調査し試しました。

棄却サンプリングは、確率密度関数$P(x)$とその上限がわかっていること（より正確には、一様に上から抑えられる関数$Q(x)$が定義できること; i.e. $\forall x, P(x) \le Q(x)$）が前提です。詳細は、[参考文献](#参考文献)を参照ください。

ざっくりいうと、乱数$rx$ に対して、 $\frac{P(rx)}{Q(rx)}$ の確率で、その乱数$rx$ を採用し、そうでない場合に棄却して、サンプリングを残す方法です。
特に、$P(x)$ は、規格化されている必要がない、という点もかなり利便性が高いです。

ここでは、扱いやすくおそらく汎用的に使えると思われる $Q(x)$ が定数関数($Q(x) = k$)の場合を前提として、実際に実機で確かめてみます。


## 実際に試してみる

早速、コードを見ていきましょう。


### 正規分布風 P(x)

```python
import numpy
from matplotlib import pyplot

rs = numpy.random.RandomState(12345)
rs
```

$P(x), Q(x)$ を以下のように定義しておきます。

```python
# P(x)
def pdf_fake_normal(x, mu, sigma):
    return numpy.exp(-0.5 * ((x - mu) / sigma) ** 2)

# Q(x)
def qdf_const(x, k=1):
    return k * numpy.ones_like(x)

```

サンプリングする関数を、以下のように定義します。最後の$px$ が求めるサンプリング結果です。（i.e. 分布$P(x)$ からのサンプリング結果です）

```python
def sampler(pdf, qdf, a=5, N=10000):
    rx = rs.uniform(-a, a, size=N)
    p = pdf(rx)
    q = qdf(rx)
    u = q * rs.uniform(0, 1, size=q.shape)
    px = rx[u <= p]
    return px
```

便利のために、密度関数の$x, y=P(x)$を生成する関数を定義しておきます。

```python
def target_maker(pdf, a):
    x = numpy.linspace(-a, a)
    _y = pdf(x) 
    y = _y / (_y.sum() * numpy.abs(x[0] - x[1]))
    return x, y
```


実際に、サンプリングするコードは、以下の通りです。pdf, qdf を定義して、sampler を実行します。

```python
%%time
N = 10000
a = 5
pdf = lambda x: pdf_fake_normal(x, mu=0, sigma=1)
qdf = lambda x: qdf_const(x, k=pdf(0))
px = sampler(pdf, qdf, a, N)
```


以下のコードで、密度関数とサンプルした乱数のヒストグラムをプロットして、正しくサンプリングできたかをグラフで確認します。

```python
x, y = target_maker(pdf, a)
pyplot.hist(px, bins=50, density=True, alpha=0.5, label='sampled')
pyplot.plot(x, y, label='target')
pyplot.legend()
pyplot.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226427/49c1f107-9f0d-94f7-56ff-cc4900039ce4.png)

すばらしい！ まさに正規分布です！！


### P(x) = x^2

今度は、分布 $P(x) := x^2$ からサンプリングしてみます。

$P(x)$ は、以下のように定義します。

```python
def pdf_x_squared(x: numpy.ndarray):
    return x ** 2
```

先ほどと同様にして、サンプリングします。pdf と、qdf を定義するパラメータ$k$ が異なっています。

```python
%%time
N = 10000
a = 5
pdf = lambda x: pdf_x_squared(x)
qdf = lambda x: qdf_const(x, k=pdf(a))
px = sampler(pdf, qdf, a, N)
```


先ほどと同様に、密度関数とサンプルした乱数のヒストグラムをプロットして、正しくサンプリングできたかを確認します。

```python
x, y = target_maker(pdf, a)
pyplot.hist(px, bins=50, density=True, alpha=0.5, label='sampled')
pyplot.plot(x, y, label='target')
pyplot.legend()
pyplot.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226427/78f2374f-c2fe-30c2-b869-ec605b39e2a7.png)

おー！、それっぽい！！

成功ですね！


# まとめ

以上のように、任意の分布からサンプリングする方法を実機で確認し、確かにサンプリングできることがわかりました。
これで、分布関数をがんばって式変形する意味がわかります！笑


# 参考文献

- [統計学&pythonレシピ | 任意の確率分布からのサンプリング2(棄却サンプリング導入編)](https://the-silkworms.com/stats_rejectsampling_intro/1050/)