statsmodels で カーネル密度推定

今回は statsmodelsを使ったカーネル密度推定を紹介します。

ドキュメントはこちらです: statsmodels.nonparametric.kde.KDEUnivariate

statsmodelsのいつもの作法で、インスタンスを作る時にデータを渡し、fitメソッドを呼び出せば完了です。
(SciPyではfitは必要なかったので注意)

fit する時に、カーネル関数の種類やバンド幅の決定方法を指定できます。
詳しくは fit メソッドのドキュメントを確認してください。

カーネル関数が 7種類指定できたり、バンド幅の指定方法が少し違ったりします。

ではやっていきましょう。 比較のため、 SciPyで推定した結果も載せていきます。
(SciPyのバンド幅指定のルールがおかしいと言う話もあるので。参考記事はこちら。)
まずはデータ生成です。


import numpy as np
import statsmodels.api  as sm
import matplotlib.pyplot as plt
from scipy.stats import norm, gaussian_kde

# 推定したい分布の真の確率密度関数
def p_pdf(x):
    return (norm.pdf(x, loc=-4, scale=2)+norm.pdf(x, loc=2, scale=1))/2


# その分布からのサンプリング
def p_rvs(size=1):
    result = []
    for _ in range(size):
        if norm.rvs() < 0:
            result.append(norm.rvs(loc=-4, scale=2))
        else:
            result.append(norm.rvs(loc=2, scale=1))
    return np.array(result)


# グラフ描写のため真の確率密度関数の値を取得しておく
xticks = np.linspace(-10, 5, 151)
pdf_values = p_pdf(xticks)

# 100件のデータをサンプリング
data = p_rvs(100)

さて、ここから SciPyと statsmodelsでそれぞれのScottのルールとSilvermanのルールでカーネル密度推定して結果を可視化します。


# SciPyのカーネル密度推定
kde_scipy_scott = gaussian_kde(data, bw_method="scott")
kde_scipy_silverman = gaussian_kde(data, bw_method="silverman")

# statusmodels のカーネル密度推定
kde_sm_scott = sm.nonparametric.KDEUnivariate(data)
kde_sm_scott.fit(kernel="gau", bw="scott")
kde_sm_silverman = sm.nonparametric.KDEUnivariate(data)
kde_sm_silverman.fit(kernel="gau", bw="silverman")

fig = plt.figure(figsize=(7, 7), facecolor="w")
ax = fig.add_subplot(1, 1, 1, title="カーネル密度推定の結果")
ax.plot(xticks, pdf_values, label="真の分布")
ax.plot(xticks, kde_scipy_scott.evaluate(xticks), label="SciPy - Scott")
ax.plot(xticks, kde_scipy_silverman.evaluate(xticks), alpha=0.5, label="SciPy - Silverman")
ax.plot(xticks, kde_sm_scott.evaluate(xticks), linestyle=":", label="statsmodels - Scott")
ax.plot(xticks, kde_sm_silverman.evaluate(xticks), label="statsmodels - Silverman")
ax.hist(data, bins=20, alpha=0.1, density=True, label="標本")
ax.legend()
plt.show()

出力がこちらです。

statsmodels で Silverman のルールを使ったのが一番良さそうですね。
ただ、これは常にそうなるわkではなく、元の分布の形や標本サイズによって何が最適かは変わってきます。

SciPyのSilvermanのルールと、statsmodelsのScottのルールが重なってしまうのは、先日の記事で書いた通り、SciPyの不備と思われます。

さて、こうしてみると、 statsmodels が一番良さそうに見えますが、当然これにも欠点があります。
SciPyやscikit-learnのカーネル密度推定では多次元のデータも扱えるのですが、
statsmodelsでは1次元のデータしか扱えません。

カーネル密度推定に関しては SciPy, scikit-learn statusmodels それぞれにメリットデメリットあるので、
目的に応じて使い分けていくのが良さそうです。
(これにそこまで高い精度が求められることも少ないので一番手軽に書ける SciPyで十分なことが多いと思います。)

scikit-learnでカーネル密度推定

今回もカーネル密度推定です。前回の記事の予告通り、今回はscikit-learnを使います。

ドキュメントはこちらです: sklearn.neighbors.KernelDensity

scikit-learnの他のモデルと同様に、import してインスタンス作って、fitしたら学習完了、となるのですが、その後の推定に注意することがあります。
version 0.23.1 時点では、他のモデルにある、predict やtransform がありません。
これは代わりに、score_samples を使います。そして戻り値は、
Evaluate the log density model on the data.
とある様に、その点での密度の”対数”です。なので、SciPyの時みたいに密度関数を得たかったら指数関数で変換する必要があります。

さて、scikit-learn でカーネル密度推定をやるメリットですが、ガウスカーネルだけでなく、全部で6種類のカーネル関数が使えます。

1点だけのデータにたいしてカーネル密度推定を実施すると、そのカーネルの形がそのまま残るのでそれを利用してカーネル関数を可視化しました。


import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity

kernel_names = [
        "gaussian",
        "tophat",
        "epanechnikov",
        "exponential",
        "linear",
        "cosine",
    ]

# x軸のメモリようの配列
xticks = np.linspace(-3, 3, 61)

fig = plt.figure(figsize=(10, 7), facecolor="w")
for i, kernel in enumerate(kernel_names, start=1):
    # 指定したカーネルでモデル作成
    kde_model = KernelDensity(bandwidth=1, kernel=kernel)
    # 1点だけのデータで学習
    kde_model.fit([[0]])
    ax = fig.add_subplot(2, 3, i, title=kernel)
    ax.set_ylim([0, 1.1])
    ax.plot(xticks, np.exp(kde_model.score_samples(xticks.reshape(-1, 1))))

plt.show()

結果がこちら。 それぞれ個性ありますね。

さて、カーネル密度推定はscikit-learnで行う場合もバンド幅の決定の問題があります。
そして、SciPyの様に自動的には推定してくれません。

そこで、scikit-learnでは、他のモデルと同じ様に、最適なパラメーターを探索することになります。
Score というメソッドで、与えられたデータに対する対数尤度の合計値が算出できるので、これを使ってグリッドサーチすると良さそうです。
他のどんなモデルもそうなのですが、訓練データと評価データはかならず分ける必要があります。
特にカーネル密度推定では、バンド幅が狭ければ狭いほど訓練データに対するスコアは高くなり、極端に過学習します。
今回は久々にグリッドサーチでもやりましょうか。
また、正解の分布はベータ分布を使います。 (正規分布だとガウスカーネルが有利すぎる気がしたので)

先ほどの 6種類のカーネル + SciPyでやってみました。


from scipy.stats import beta
from sklearn.model_selection import GridSearchCV

# グリッドサーチでバンド幅決めてベストのモデルを返す関数
def search_bandwidth(data, kernel="gaussian", cv=5):
    bandwidth = np.linspace(0.1, 2, 190)
    gs_model = GridSearchCV(
                KernelDensity(kernel=kernel),
                {'bandwidth': bandwidth},
                cv=cv
            )
    gs_model.fit(data.reshape(-1, 1))
    return gs_model.best_estimator_

# 正解の分布
frozen = beta(a=3, b=2, loc=-1, scale=2)
# データ生成
data = frozen.rvs(200)

# x軸のメモリようの配列
xticks = np.linspace(-1.1, 1.1, 221)


fig = plt.figure(figsize=(10, 11), facecolor="w")
for i, kernel in enumerate(kernel_names, start=1):
    # 指定したカーネルでモデル作成
    kde_model = search_bandwidth(data, kernel=kernel)
    bw_str = str(kde_model.bandwidth.round(2))
    ax = fig.add_subplot(3, 3, i, title=kernel + " h=" + bw_str)
    # ax.set_ylim([0, 1.1])
    ax.plot(xticks, np.exp(kde_model.score_samples(xticks.reshape(-1, 1))))
    ax.plot(xticks, frozen.pdf(xticks), linestyle="--")

# おまけ。scipy
kde = gaussian_kde(data)
ax = fig.add_subplot(3, 3, 7, title="SciPi")
ax.plot(xticks,  kde.evaluate(xticks))
ax.plot(xticks, frozen.pdf(xticks), linestyle="--")

plt.show()

結果がこちら。 一番手軽に使える SciPyが頑張ってますね。
全データを学習に使えるのも大きいかもしれません。

SciPyのガウシアンカーネル密度推定に実装されているバンド幅の算出方法について

僕の環境の SciPyは version 1.4.1 なのでこの記事ではそれを想定して読んでください。(将来のバージョンでは修正されるかも。)

さて、カーネル密度推定においてバンド幅(パラメーターh)を決めるのに使われる、
Scottのルールと、Silvermanのルールについて紹介しようと思っていたのですが、
SciPyの実装を精査していく中でちょっと挙動がおかしいことに気付きました。

それぞれのルールの原典がどっちも無料では見れないのでそれを確認できないのですが、
SciPyで実装されているSilvermanのルールが、どうやら他で言うScottのルールのようです。
SciPyで実装されているScottのルールは、これ一体何なんでしょう。
他で言われているSilvermanのルールはもっと狭いバンド幅で推定するはずです。

Scottのルールの原典:
D.W. Scott, “Multivariate Density Estimation: Theory, Practice, and Visualization”, John Wiley & Sons, New York, Chicester, 1992.
Silvermanのルールの原典:
B.W. Silverman, “Density Estimation for Statistics and Data Analysis”, Vol. 26, Monographs on Statistics and Applied Probability, Chapman and Hall, London, 1986.

本当は正しい(と言うか、この二つのルールも経験則的なルールであって、それを正しいと呼ぶかは疑問符が付きますが)それぞれの発案者が提案したルールを紹介したいのですが、
普段の利用を考えると、scipy.stats.gaussian_kde で、bw_method に指定した文字列によって挙動がどう変わるのかを正確に理解している方が役に立つと思います。
そのため、この記事のタイトル通り、発案者が何を提唱したかは一旦脇に置いといて、文字列として、”scott”と”silverman”を指定したらプログラムがどう計算するかを紹介します。

ドキュメントはこちらです: scipy.stats.gaussian_kde
GitHubのコードもかなり参考になりました。

さて、推定したバンド幅はその2乗の値が covariance という属性に格納されています。(1次元データの場合。)
これは、 学習したデータの不変分散に、 factor という属性の”2乗を”掛けることで計算される実装になっています。
バンド幅hで考えると、不変分散の平方根に factor を掛けたものになりますね。

で、その factor の計算方法が、指定した “scott”と”silverman”で変わります。
データの件数を$n$、データの次元を$d$とすると、
“scott”を指定したときは、
$$
factor = n^{\frac{-1}{d+4}}
$$
です。1次元の場合は$d=1$なので$n^{-0.2}$ですね。

“silverman”を指定したときは、
$$
factor = \left(\frac{d+2}{4}*n\right)^{\frac{-1}{d+4}}
$$
になります。
$d=1$の時は、$(3/4)^{-0.2}=1.05922\cdots$なので、だいたい、$1.06n^{-0.2}$くらいになります。
他のサイトの情報を見ていると、 0.9とか1.06とかの定数が出てくるのですが、そのうち1.06がどこから登場するのかはこれで分かりますね。
ちなみに statsmodelsのカーネル密度推定では 1.059がハードコーディングされています。

で、この 1.06は本当は Scott のルールで使われる定数のはずなのですが、なぜかSciPyでは”silverman”を指定したときに登場します。
1.06の代わりに0.9くらいの値を使うのがSilvermanのルールのはずなのですが。

ただ、以上で、SciPyがどの様にバンド幅を決定しているかはわかったので、その通りの挙動をしていることをコードで確認しておきましょう。
とりあえずデータを準備します。


from scipy.stats import gaussian_kde
from scipy.stats import norm
import numpy as np

# 50件のデータをサンプリング
data = norm(loc=2, scale=2).rvs(50)

print("標本の不偏分散: ", np.var(data, ddof=1))
# 標本の不偏分散:  4.2218467749903645

まずは “scott” の確認です。(初期値なので、コード中に”scott”の文字列は登場しません。)


# Scottのルールでカーネル密度推定
kde = gaussian_kde(data)

# バンド幅 hの2乗を確認
print("バンド幅 h^2: ", kde.covariance[0, 0])
# バンド幅 h^2:  0.8829059945819668

# Scottの因子が、データ件数の -0.2乗に等しいことを確認
print(kde.scotts_factor())
# 0.45730505192732634
print(len(data) ** (-0.2))
# 0.45730505192732634

# データの普遍分散に Scottの因子の2乗を掛けてみると、先ほどのバンド幅の2乗に等しい
print((kde.scotts_factor()**2) * np.var(data, ddof=1))
# 0.8829059945819668

以上の通り、ドキュメント通りの実装になっていました。

次が bw_method=”silverman” の場合です。ごくわずか謎の端数が出てますが、こちらもドキュメント通りの実装であることが確認できます。


# カーネル密度推定 (silvermanのルール)
kde = gaussian_kde(data, bw_method="silverman")

# バンド幅 hの2乗を確認
print("バンド幅 h^2: ", kde.covariance[0, 0])
# バンド幅 h^2:  0.9905809235665322

# Silvermanの因子が、定義通りの計算結果に等しいことを確認
print(kde.silverman_factor())
# 0.48438841363348917
print((kde.n * (kde.d + 2) / 4.)**(-1. / (kde.d + 4)))
# 0.4843884136334891

# データの普遍分散に Silvermanの因子の2乗を掛けてみると、先ほどのバンド幅の2乗に等しい
print((kde.silverman_factor()**2) * np.var(data, ddof=1))
# 0.9905809235665322

ちなみに、 statsmodelsのドキュメントを見てみましょう。
statsmodels.nonparametric.kde.KDEUnivariate.fit

The bandwidth to use. Choices are:
“scott” – 1.059 * A * nobs ** (-1/5.), where A is min(std(X),IQR/1.34)
“silverman” – .9 * A * nobs ** (-1/5.), where A is min(std(X),IQR/1.34)
“normal_reference” – C * A * nobs ** (-1/5.), where C is calculated from the kernel. Equivalent (up to 2 dp) to the “scott” bandwidth for gaussian kernels. See bandwidths.py
If a float is given, it is the bandwidth.

この様に、 “scott” のほうが、 $1.059*n^{(-0.2)}$ を掛けています。SciPyでは”silverman”の時に登場するものですね。
また、データの標準偏差よりも、IQR (これは 第3四分位数から第1四分位数を引いた値です)を1.34で割ったものが小さかったらそっちを使うと言うルールも実装されています。
SciPyのコードを見るとこれも入ってないです。

ちょっと色々疑わしいところはあるのですが、ぶっちゃけた話、ヒストグラムなどを使って可視化する時におまけにつけることがあるくらいのものなので、
あんまり気にしないで使って行こうかなと思っています。

ただ、statsmodelsの方が正しい値出してくれそうなのでこれもそのうちコードを載せて紹介しましょう。
次の記事ではガウスカーネル以外を使ったカーネル密度推定を先にやりたいので、Scikit-learnを使った実装を紹介する予定です。

SciPyによるカーネル密度推定

前回の記事に続いてカーネル密度推定の話です。
参考: カーネル密度推定法の紹介

前の記事では数式に沿って実装するため、numpyで書きましたが、普段の利用では、SciPyにクラスが実装されているのでこれを使えます。

参考: scipy.stats.gaussian_kde
クラス名から分かる通り、カーネルはガウスカーネルがあらかじめ指定されており、他の関数は使えません。
実際のところ、ガウスカーネル以外が使えなくて困ることは滅多にないでしょう。

使い方ですが、 インスタンスを作るときに、データセットを渡すだけで完了します。
とてもありがたいことに、バンド幅の指定が要りません。
デフォルトでは Scottのルールと呼ばれる方法で決定されます。
bw_methodという引数で指定することで、もっと直接的な指定や、Silvermanのルールを使うことも可能です。
ScottのルールとSilvermanのルールがそれぞれ具体的にどの様にバンド幅を決めているのかは後日のブログで紹介するとして、
とりあえず、正規分布から生成したデータで使ってみましょう。


from scipy.stats import gaussian_kde
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt

# 50件のデータをサンプリング
data = norm(loc=2, scale=2).rvs(50)

# グラフ描写のため真の確率密度関数の値を取得しておく
xticks = np.linspace(-5, 9, 141)
pdf_values = norm(loc=2, scale=2).pdf(xticks)

# カーネル密度推定
kde = gaussian_kde(data)
# evaluate には pdfというエイリアスがあるのでそちらでも可
estimated_value = kde.evaluate(xticks)

# 可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ax.plot(xticks, pdf_values, label="真の分布")
ax.plot(xticks, estimated_value, label="推定値",  linestyle="--")
# 標本データ
for x in data:
    ax.vlines(x, ymin=0, ymax=0.01, alpha=0.5)
ax.set_ylim((0, 0.25))
ax.legend()
plt.show()

出力結果がこちらです。

とても手軽に良い推定が得られていますね。

ちなみに、推定してくれたバンド幅ですが、covariance という属性にその情報を持っています。


# バンド幅の確認 (hの2乗なので注意)
print(kde.covariance)

# [[0.79905207]]

多次元にも対応するため、分散共分散行列の形でデータを持てる様になっています。
今回は1次元データなので、ちょっと面倒ですが、無駄に2次元配列になってる中から要素を持ってくる必要があります。
varianceなので、この値は分散、前回の記事のバンド幅は標準偏差なので、平方根を取る必要があることも注意です。

せっかく推定してくれたバンド幅が取れるので、前回スクラッチで描いたコードの推定値と結果が一致することも見ておきましょう。
前回の記事のコピーだと、関数名が重複してしまって面倒なので適当に、my_kdeにリネームしました。


# ガウスカーネル
def gaussian_kernel(u):
    return(np.exp(-(u**2)/2)/np.sqrt(2*np.pi))


# ガウスカーネルを用いたカーネル密度推定
def my_kde(x, h, dataset):
    return sum([gaussian_kernel((x-d)/h) for d in dataset])/(len(dataset)*h)


# バンド幅の取得
h = np.sqrt(kde.covariance[0][0])
# 可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ax.plot(xticks, pdf_values, label="真の分布")
ax.plot(xticks, kde.evaluate(xticks), label="scipy", linestyle="--")
ax.plot(xticks, my_kde(xticks, h, data), label="スクラッチ", linestyle=":")
ax.legend()
plt.show()

SciPy で出した結果と、スクラッチで書いてるやつの結果が(一致してることの確認なので)重なってしまっていて、見にくいですが、
一応線の種類を変えてあるので、よく見るとしっかり同じところに線が引かれているのが分かります。

カーネル密度推定法の紹介

今回からしばらくカーネル密度推定法について書いていきます。

参考文献ですが「パターン認識と機械学習 上 (PRML 上)」の 119ページ、 2.5.1 カーネル密度推定法 あたりがわかりやすいのではないでしょうか。
ただし、この本だと最後にパラメーター$h$の決め方が重要だ、ってところまで書いて、$h$の決め方を紹介せずに次の話題に移っているので、
これからの数回の記事の中でそこまで説明できたらいいなと思っています。

さて、カーネル密度推定法を大雑把に言うと、
「確率密度が全然わからない確率分布から標本が得られたときに、その標本からできるだけいい感じに元の確率密度を推定する方法」です。
元の確率密度が正規分布だとかポアソン分布だとかだけでもわかっていたら、その期待値や分散などのパラメーターだけ推定すればいい(これをパラメトリックという)のですが、
それさえもわからないときに使われる手法であり、ノンパラメトリックなアプローチと呼ばれます。

できるだけ一般的な形でカーネル密度推定法の方法を書くと次のような形になります。

ある$D$次元のユークリッド空間の未知の確率密度$p(\mathbf{x})$から、観測値の集合$\{\mathbf{x}_n\}\ \ (n=1,2, \cdots, N)$が得られたとします。

そして、次の性質を満たすカーネル関数$k(\mathbf{u})$を選びます。
$$
\begin{align}
k(\mathbf{u}) &\geq 0,\\
\int k(\mathbf{u}) d\mathbf{u} &= 1
\end{align}
$$
このとき、$p(\mathbf{x})$を次で推定するのがカーネル密度推定法です。($h$はパラメーター。)
$$
p(\mathbf{x}) = \frac{1}{N}\sum_{n=1}^{N}\frac{1}{h^D}k(\frac{\mathbf{x}-\mathbf{x}_n}{h}).
$$

カーネル密度関数は要は確率密度間数ですね。理論上は何でもいいそうですが、現実的には原点周りで値が大きく、そこから外れると0に近く(もしくは等しく)なる関数が好都合です。
正規分布の確率密度関数(ガウスカーネル)が頻繁に用いられます。

上記の式が導出されるまでの話はPRMLに書いてあります。
その部分を省略してしまったので代わりに上の定義通り計算して、未知の確率分布が近似できている様子を確認しておきましょう。
できるだけライブラリ使わないでやります。(ただ、正解データの分布作成でscipyは使っていますし、numpyやmatplotlibは例外とします。)
カーネル関数にはガウスカーネルと、単位超立方体内で1、その外で0を取る関数(パルツェン窓というらしい)の2種類でやってみます。

元となる分布は正規分布2個を足した、2つ山のある分布で実験します。($D$は1です)

最初に正解になる分布関数の定義と、そこからのデータサンプリングを済ませておきます。


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


# 推定したい分布の真の確率密度関数
def p_pdf(x):
    return (norm.pdf(x, loc=-4, scale=2)+norm.pdf(x, loc=2, scale=1))/2


# その分布からのサンプリング
def p_rvs(size=1):
    result = []
    for _ in range(size):
        if norm.rvs() < 0:
            result.append(norm.rvs(loc=-4, scale=2))
        else:
            result.append(norm.rvs(loc=2, scale=1))
    return np.array(result)


# グラフ描写のため真の確率密度関数の値を取得しておく
xticks = np.linspace(-10, 5, 151)
pdf_values = p_pdf(xticks)

# 50件のデータをサンプリング
X = p_rvs(50)

そして、それぞれのカーネル関数と、それを使ったカーネル密度関数を定義しておきます。


# パルツェン窓
def parzen_window(u):
    if np.abs(u) <= 0.5:
        return 1
    else:
        return 0


# ガウスカーネル
def gaussian_kernel(u):
    return(np.exp(-(u**2)/2)/np.sqrt(2*np.pi))


# パルツェン窓を用いたカーネル密度推定
def parzen_kde(x, h, dataset):
    return sum([parzen_window((x-d)/h) for d in dataset])/(len(dataset)*h)


# ガウスカーネルを用いたカーネル密度推定
def gaussian_kde(x, h, dataset):
    return sum([gaussian_kernel((x-d)/h) for d in dataset])/(len(dataset)*h)

あとは可視化です。 今回は パラメーター$h$は$0.5, 1, 2$の3つ試しました。


fig = plt.figure(figsize=(12, 12), facecolor="w")

# 3種類のhで実験
for i, h in enumerate([0.5, 1, 2]):
    parzen_kde_values = [parzen_kde(x, h, X) for x in xticks]
    gaussian_kde_values = [gaussian_kde(x, h, X) for x in xticks]

    ax = fig.add_subplot(3, 2, i*2+1, facecolor="w")
    ax.set_title(f"パルツェン窓, h={h}")
    ax.plot(xticks, pdf_values, label="真の分布")
    ax.plot(xticks, parzen_kde_values, label="パルツェン窓")
    # 標本データ
    for x in X:
        ax.vlines(x, ymin=0, ymax=0.01)
    ax.set_ylim((0, 0.37))
    ax.legend()

    ax = fig.add_subplot(3, 2, i*2+2, facecolor="w")
    ax.set_title(f"ガウスカーネル, h={h}")
    ax.plot(xticks, pdf_values, label="真の分布")
    ax.plot(xticks, gaussian_kde_values, label="ガウスカーネル")
    # 標本データ
    for x in X:
        ax.vlines(x, ymin=0, ymax=0.01)
    ax.set_ylim((0, 0.37))
    ax.legend()

plt.show()

出力結果がこれです。

青い線が真の分布、黒くて短い縦線たちが推定に使った標本データ、
オレンジの線が推定した分布です。

「元の確率密度関数の形が全くわからない」という前提から出発している割にはうまく推定できているのもあると言えるのではないでしょうか。
また、カーネル密度関数の選び方、$h$の決め方がそれぞれ重要であることも分かりますね。

やってることといえば、得られた標本の付近で大きめの値を取る関数を全ての標本について計算してその平均を取っているだけなのですが、シンプルな割に強力です。

棄却法の例

前回の記事で紹介した棄却法を実際にやってみましょう。

今回の例はベータ分布です。(とりうる値の範囲が有限のものの方が適用しやすいので)
とりあえず、$B(2,3)$でやってみましょう。
確率密度関数は区間$x\in[0,1]$の範囲では次の式で表されます。(それ以外の$x$に対しては$0$です)。
$$
f(x) = \frac{x(1-x)^2}{B(2, 3)} = 12x(1-x)^2.
$$
この関数は$x=1/3$で最大値$f(1/3)=16/9$をとります。
(単純な式なので微分してすぐに確認できます。)

さて、実際にプログラムで実行してみたのがこちらです。
念のためですが、今回の目的は前回の記事で紹介したアルゴリズムで目的とする乱数が得られることを確認することです。
単にベータ分布に従う乱数が必要な場合は、scipyのrvs関数を使いましょう。


# ベータ分布の確率密度関数
f = beta.freeze(a=2, b=3)


def beta_rejection_sampling():
    while True:
        u, v = uniform().rvs(2)
        x = u  # 今回は 0<= x <= 1 なのでu をそのまま使用
        y = 16/9 * v
        if y <= f.pdf(x):
            return x


# 棄却法で10000データ生成する
data = [beta_rejection_sampling() for _ in range(10000)]

# 生成された乱数のヒストグラムと、確率密度関数を可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ax.hist(data, bins=100, density=True, label="棄却法により生成")
ax.set_xlim([0, 1])
x = np.linspace(0, 1, 100)
ax.plot(x, f.pdf(x), label="確率密度関数")
ax.legend()
plt.show()

出力がこちら。

しっかり機能していますね。

フォンノイマンの棄却法

今回も乱数を生成するお話。
累積分布関数の逆関数が求まるなら逆関数法、正規分布の場合はボックス=ミュラー法が使えるという話を書きましたが、
もっと一般の分布で使える方法として、フォンノイマンの棄却法というのがあることを最近知りました。

Wikipediaでは英語版のみページがあるようです:Rejection sampling
(自然科学の統計学に紹介されている 別名法もこれは離散版ですがアイデアが似てるので参考になるかも)

これは次のステップで行います。
まず、生成したい分布の確率密度関数$f(x)$の最大値$M$を求めておきます。
また、取得する乱数の区間$[x_{min}, x_{max}]$をきめます。
(ベータ分布や2項分布のような有限区間の値しか取らない乱数なら容易ですが、そうでない場合は十分大きな範囲をとって適当なところで区切るしかないですね)

そして、次の手順で乱数を生成します。
1. 標準一様分布$U(0, 1)$から二つの乱数$u, v$を生成する。
2. $x = x_{min}+(x_{max}-x_{min})u$ を計算し、区間$[x_{min}, x_{max}]$の乱数を得る。
3. $y = M*v$ を計算し、 $y$と$f(x)$を比較する。
4. 結果が$y<=f(x)$であれば、乱数として$x$を採用し、そうでない場合は、二つの乱数生成に戻ります。 $x$が乱数として採用される確率が$f(x)$の値に比例するため、 結果として確率密度関数$f$に従う乱数を得ることができます。

ボックス=ミュラー法

Scipyが使える今となっては使う機会はほぼありませんが、
一様分布から正規分布に従う乱数を作成できる方法である、ボックス=ミュラー法(Box–Muller’s method)を紹介します。

正規分布は累積分布関数やその逆関数が初等関数では表現できず、
最近紹介した逆関数法で乱数を生成するのは少々困難です。

そこでこのボックス=ミュラー法が使われます。
参考:自然科学の統計学(東京大学出版会)の 11.3 正規乱数の発生法

まず、確率変数$X$,$Y$が互いに独立で、共に$(0,1)$上の一様分布に従うとします。
この時、
$$
\begin{align}
Z_1 & = \sqrt{-2\log{X}}\cos{2\pi Y},\\
Z_2 & = \sqrt{-2\log{X}}\sin{2\pi Y}
\end{align}
$$
とすると、$Z_1$と$Z_2$は標準正規分布$N(0,1)$に従う互いに独立な確率変数になります。

厳密な証明は今回は省略します。
ただ、数式をみれば、正規分布の確率密度関数が指数関数の形をしているので$\log$が出てくるのも、
円周率も出てくるなど円に関係しそうな気配があるので三角関数が出てくるのもなんとなく納得性があります。

ただ、「互いに独立な」ってのは正直驚きます。
$Z_1^2+Z_2^2=-2\log{X}$って関係式が成り立つのに独立ってことはないんじゃないかと思えますね。 
(この$X$が定数ではないのがキモで、$0<X<1$から、$-2\log{X}$が
非常に大きな値も含めて実に自由に動くので独立性が生まれるようです。)

ここだけ実験して相関係数が0に近いことを確認しておきましょう。


import numpy as np

# 一様分布に従うX, Yをそれぞれ10000個生成
X = np.random.rand(10000)
Y = np.random.rand(10000)

# Z_1, Z_2 を計算
Z_1 = np.sqrt(-2*np.log(X))*np.cos(2*np.pi*Y)
Z_2 = np.sqrt(-2*np.log(X))*np.sin(2*np.pi*Y)

# 相関係数を算出
print(np.corrcoef(Z_1, Z_2)[0, 1])
# -0.0043281961119066865

確かに独立っぽいですね。
散布図等も出してみましたが何か特別な関係性は見当たらないようです。

逆関数法で指数分布に従う乱数を生成する

前回の記事で逆関数法を紹介したので具体的な分布で試そうという趣旨の記事です。
今回は指数分布を使います。(累積分布関数もその逆関数も容易に計算できるから)

まず、指数分布の確率密度関数から。指数分布はパラメータ$\lambda$を一つ持ちます。
$$
f(x) = \left\{
\begin{matrix}
\lambda e^{-\lambda x} & (x \leq 0)\\
0 & (x < 0) \end{matrix} \right. $$ そして、累積分布関数は次のようになります。 $$ F(x) = \left\{ \begin{matrix} 1-e^{-\lambda x} & (x \leq 0)\\ 0 & (x < 0) \end{matrix} \right. $$ さらにその逆関数は、こうなります。 $$ F^{-1}(u) = - \frac{1}{\lambda}\log(1-u) \quad (0\leq u <1). $$ それではPythonで、$[0,1)$の乱数を10000個くらい生成して、 逆関数法で$f(x)$に従う乱数が得られるのか実験してみましょう。 今回は単純のため、$\lambda=1$とします。 Pythonでやってみたコードと結果がこちら。


import numpy as np
import matplotlib.pyplot as plt

# 標準一様分布 に従う乱数を10000個生成する
U_data = np.random.rand(10000)
# 逆関数方法で指数分布に従うデータ生成
result_data = -np.log(1 - U_data)
# ヒストグラムと確率密度関数をプロットする
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ax.hist(result_data, bins=100, density=True)
ax.plot(np.arange(0, 8, 0.1), np.exp(-1*np.arange(0, 8, 0.1)))
plt.show()

出力。

データ件数が多いので非常に綺麗に出ましたね。

任意の確率分布に従う乱数を逆関数法で一様分布から生成する

学生時代にExcelで色々シミュレーションしていたときによく使った方法です。
エクセル関数で生成できる乱数がrand()で得られる一様分布の乱数くらいしかなく、
他の分布に従う乱数が欲しい時はそこから変換して作っていました。
(Pythonを使うようになってから、Scipyに実装されている確率分布を使えば良いので、滅多に行うことがなくなりました。)

ただ、逆関数法という名前がついていることは最近知りました。
逆関数法 出典: フリー百科事典『ウィキペディア(Wikipedia)』
英語では inverse transform sampling というようです。

前回の記事同様、少々厳密性は犠牲にして、登場する各関数等は(連続性や滑らかさなどの)適切な仮定を満たしているものとします。

まず前提として、区間$[0,1]$上の一様分布に従う乱数$U$は取得できるものとします。
$$
\begin{align}
f_U(u) = \left\{\begin{matrix}
1& u \in[0,1]\\
0& u \notin[0,1]
\end{matrix}\right.
\end{align}
$$

そして欲しい乱数$X$が従う確率密度関数を$f(x)$とし、累積分布関数を$F(x)$とします。
$$
F(x) = \int_{-\infty}^x f(y) dy.
$$

そして$F$の逆関数を$F^{-1}$とした時、
$X=F^{-1}(U)$は、確率密度関数$f(x)$に従う確率変数となります。

証明は前回の記事で紹介した、確率変数の変換を使います。
この記事中の、$g$に相当するのが$F^{-1}$であることに注意して計算します。
まず前回の記事で紹介した返還式に、この記事中で使っている記号に注意しながら代入します。
$$
f_X(x) = f_U(F(x))\frac{du}{dx}.
$$
ここで、累積分布関数の性質より$0\leq F(x)\leq 1$より、$f_U(F(x))=1$です。
さらに、
$$
\frac{du}{dx} = \frac{d}{dx}F(x) = f(x)
$$
となることから、次の式が示されました。
$$
f_X(x) = f(x).
$$

累積分布関数(と、その逆関数)が具体的に計算できる確率分布についてはこれは非常に便利な式です。
学生時代は個人的な調べ物のため、指数分布を使うシミュレーションをよくやっていたので、よく使っていました。