定義
指数分布はパラメータ \(\lambda\) を用いて以下のように定義される。
$$\mathrm{Exp}(x|\lambda)=\lambda e^{-\lambda x}$$
基本的な性質
平均
$$\mathbb{E}[x]=\frac{1}{\lambda}$$
分散
$$\mathrm{var}[x]=\frac{1}{\lambda^{2}}$$
最尤推定
公式
データ点の集合 \(\mathcal{D}=(x_{1},x_{2},\ldots,x_{N})^{\mathrm{T}}\) を、最尤推定により指数分布にフィッティングすることを考える。このとき、パラメータ \(\lambda\) は
$$\lambda=\frac{N}{\sum_{n=1}^{N}x_{n}}\tag{1}$$
で与えられる。
証明
データについての尤度関数は
$$p(\mathcal{D}|\lambda)=\prod_{n=1}^{N}\mathrm{Exp}(x_{n}|\lambda)$$
$$=\lambda^{N}\exp{\left(-\lambda\sum_{n=1}^{N}x_{n}\right)}$$
で与えられる。両辺の対数を取って、対数尤度関数を
$$L(\lambda|D)=\ln{p(\mathcal{D}|\lambda)}$$
$$=N\ln{\lambda}-\lambda\sum_{n=1}^{N}x_{n}$$
と置く。これを \(\lambda\) について微分したものを \(0\) と置くと
$$\frac{\partial L(\lambda|D)}{\partial \lambda}=\frac{N}{\lambda}-\sum_{n=1}^{N}x_{n}=0$$
より
$$\lambda=\frac{N}{\sum_{n=1}^{N}x_{n}}$$
パラメータ変換
公式
指数分布 \(\mathrm{Exp}(x|\lambda_{\alpha})\) に従うデータ点の集合 \(\mathcal{D}=(x_{1},x_{2},\ldots,x_{N})^{\mathrm{T}}\) に、変換 \(f\) を加えて、指数分布 \(\mathrm{Exp}(x|\lambda_{\beta})\) に従うデータ点の集合 \(f(\mathcal{D})=(f(x_{1}),f(x_{2}),\ldots,f(x_{N}))^{\mathrm{T}}\) を生成することを考える。このとき
$$f(x_{n})=\frac{\lambda_{\alpha}}{\lambda_{\beta}}x_{n}$$
である。とくに、\(\lambda_{\beta}=1\) のとき
$$f(x_{n})=\lambda_{\alpha}x_{n}$$
と変換することで、データが \(\mathrm{Exp}(x|1)=e^{-x}\) に従うようにスケーリングを行うことができる。
証明
変換後のデータ点の集合 \(f(\mathcal{D})\) を最尤推定した時、式 \((1)\) より
$$\lambda_{\beta}=\frac{N}{\sum_{n=1}^{N}f(x_{n})}\tag{2}$$
が得られるとする。同様に、変換前のデータ点の集合 \(\mathcal{D}\) について
$$\lambda_{\alpha}=\frac{N}{\sum_{n=1}^{N}x_{n}}$$
$$N=\lambda_{\alpha}\sum_{n=1}^{N}x_{n}$$
が得られるので、これを式 \((2)\) に代入して
$$\lambda_{\beta}=\frac{1}{\sum_{n=1}^{N}f(x_{n})}\lambda_{\alpha}\sum_{n=1}^{N}x_{n}$$
$$\lambda_{\beta}\sum_{n=1}^{N}f(x_{n})=\lambda_{\alpha}\sum_{n=1}^{N}x_{n}$$
$$\sum_{n=1}^{N}f(x_{n})=\sum_{n=1}^{N}\frac{\lambda_{\alpha}}{\lambda_{\beta}}x_{n}$$
となり、これより
$$f(x_{n})=\frac{\lambda_{\alpha}}{\lambda_{\beta}}x_{n}$$
が得られる。
乱数生成
公式
指数分布 \(\mathrm{Exp}(y|\lambda)\) に従う乱数 \(y\) は、標準一様分布 \(U(z|0,1)\) に従う乱数 \(z\) を
$$y=-\frac{1}{\lambda}\ln{z}$$
と変換することで生成できる。
証明
この乱数生成は、逆関数法によるサンプリングを利用している。逆関数法については
を参照。
ここではまず、指数分布 \(\mathrm{Exp}(x|\lambda)\) の累積分布関数について考える。
累積分布関数
指数分布 \(\mathrm{Exp}(x|\lambda)\) の累積分布関数は以下のように与えられる。
$$1-e^{-\lambda x}$$
累積分布関数の導出
$$\int_{0}^{x}\mathrm{Exp}(t|\lambda)dt=\int_{0}^{x}\lambda e^{-\lambda t}dt$$
$$=\left[-e^{-\lambda t}\right]_{0}^{x}=1-e^{-\lambda x}$$
逆関数法の導出
累積分布関数 \(z=h(y)=1-e^{-\lambda y}\) の逆関数を考えて
$$e^{-\lambda y}=1-z$$
$$-\lambda y=\ln{(1-z)}$$
$$y=-\frac{1}{\lambda}\ln{(1-z)}$$
を得る。ここでは \(0<z<1\) であることに注意されたい。また、\(z\) が標準一様分布に従うとき、\(1-z\) も同様に標準一様分布に従うため、\((1-z)\rightarrow z\) という置き換えを行ってもよい。したがって
$$y=-\frac{1}{\lambda}\ln{z}$$
を得る。
Pythonによる実験
最後に、以上の結果を用いてPythonによる実験を行う。下記のコードでは
- \(\lambda=2\) の指数分布に従うサンプルを生成
- サンプルから \(\lambda\) の値を推定
- \(\lambda=4\) となるようにサンプルを変換
- \(\lambda\) を再推定
を行い、その結果を図示している。サンプルのヒストグラムと確率分布を同時にプロットする方法については
を参照。また、累積分布関数の式から
$$\int_{\alpha}^{\beta}\mathrm{Exp}(x|\lambda)dx=\int_{0}^{\beta}\mathrm{Exp}(x|\lambda)dx-\int_{0}^{\alpha}\mathrm{Exp}(x|\lambda)dx$$
$$=1-e^{-\lambda \beta}-(1-e^{\lambda \alpha})$$
$$=e^{\lambda \alpha}-e^{-\lambda \beta}$$
の関係を利用した。
import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import expon
n = 100000 # N
z = np.random.random(n) # z ~ U(0,1) を生成
y = -np.log(z) / 2 # y ~ Exp(2) を生成
lamb2 = n / y.sum() # lambda を最尤推定
print("lambda_2 = {}".format(lamb2)) # lambda_2 = 1.9920861363139992
_, b, _ = plt.hist(y, bins=1000, alpha=0.5, label="data from 2exp(-2y)")
fy = y * 2 / 4 # lambda の値を 2 -> 4 に変換
lamb4 = n / fy.sum() # lambda を再推定
print("lambda_4 = {}".format(lamb4)) # lambda_4 = 3.9841722726279984
plt.hist(fy, bins=b, alpha=0.5, label="data from 4exp(-4y)")
cum_fit2 = np.exp(-lamb2*b) # 各ビンの端における、Exp(y|lamb2) の累積分布関数の第２項
cum_fit4 = np.exp(-lamb4*b) # 各ビンの端における、Exp(y|lamb4) の累積分布関数の第２項
b_mid = [(b[i+1] + b[i]) / 2 for i in range(len(b)-1)] # 各ビンの中心
cum_fit2_plot = [(cum_fit2[i] - cum_fit2[i+1]) * n for i in range(len(b_mid))] # 各ビン区間の確率値
cum_fit4_plot = [(cum_fit4[i] - cum_fit4[i+1]) * n for i in range(len(b_mid))] # 各ビン区間の確率値
plt.plot(b_mid, cum_fit2_plot, label="{:1.2f}exp(-{:1.2f}y) (expanded)".format(lamb2, lamb2))
plt.plot(b_mid, cum_fit4_plot, label="{:1.2f}exp(-{:1.2f}y) (expanded)".format(lamb4, lamb4))
plt.legend()
plt.title("Exponential distribution")
plt.xlabel("y")
plt.ylabel("number")
plt.savefig("exp_dist.png")
plt.show()
コメント