はじめに
pythonの勉強中なので、少しプログラムを書いてみました。
途中、少し、計算過程を記載していますが、、温かい目でお願いします。
目的・テーマ
以下に示す２母数ワイブル分布の母数推定について勉強してみました。
ここで、α、βをそれぞれ尺度母数、形状母数です。
分布の形状を決めるパラメータですね。
f(x;\alpha,\beta) = \frac{\beta}{\alpha} \left( \frac{x}{\alpha} \right)^{\beta - 1} \exp \left[ - \left( \frac{x}{\alpha} \right)^\beta \right]
今回は最尤法を用いて２母数ワイブル分布の母数推定にチャレンジしようと思います。
理論編
２母数ワイブル分布が上記の式で表される時、尤度関数は以下のようになります。
\begin{aligned}
L(\alpha,\beta) &= \prod_{i=1}^{n}f(x_i;\alpha,\beta) \\
&= \prod_{i=1}^{n}\frac{\beta}{\alpha} \left( \frac{x_i}{\alpha} \right)^{\beta - 1} \exp \left[ - \left( \frac{x_i}{\alpha} \right)^\beta \right] -①
\end{aligned}
対数尤度関数は以下のようになります。
\log L(\alpha,\beta) = \log \prod_{i=1}^{n}\frac{\beta}{\alpha} \left( \frac{x_i}{\alpha} \right)^{\beta - 1} \exp \left[ - \left( \frac{x_i}{\alpha} \right)^\beta \right] -②
まずは式②の右辺をわかりやすくしましょう。
\begin{aligned}
\log L(\alpha,\beta) &= \log \prod_{i=1}^{n}\frac{\beta}{\alpha} \left( \frac{x_i}{\alpha} \right)^{\beta - 1} \exp \left[ - \left( \frac{x_i}{\alpha} \right)^\beta \right] \\
&= \log (\frac{\beta}{\alpha})^n \prod_{i=1}^{n} \left( \frac{x_i}{\alpha} \right)^{\beta - 1} \exp \left[ - \left( \frac{x_i}{\alpha} \right)^\beta \right] \\
&= \log (\frac{\beta}{\alpha})^n \prod_{i=1}^{n} \left( \frac{x_i}{\alpha} \right)^{\beta - 1} \prod_{i=1}^{n}\exp \left[ - \left( \frac{x_i}{\alpha} \right)^\beta \right] \\
&= n\log (\frac{\beta}{\alpha}) + (\beta - 1)\sum_{i=1}^{n} \log(\frac{x_i}{\alpha}) - \sum_{i=1}^{n}(\frac{x_i}{\alpha})^\beta
-②'
\end{aligned} \\
以下の式を用いて、α、βの値を求めます。
\begin{aligned}
\frac{\partial \log L(\alpha,\beta)}{\partial \alpha} &= 0 -③ \\
\frac{\partial \log L(\alpha,\beta)}{\partial \beta} &= 0 -④
\end{aligned}
式②'を踏まえて、式③、④の左辺を計算してみましょう。
（式②'をαで偏微分する）
\begin{aligned}
\frac{\partial \log L(\alpha,\beta)}{\partial \alpha} &= -\frac{n}{\alpha} +(1-\beta)\sum_{i=1}^{n}(\frac{1}{\alpha})+\sum_{i=1}^{n}(\frac{\beta x_i^\beta}{\alpha^{\beta+1}}) \\
&= -\frac{n}{\alpha} + \frac{n}{\alpha} - \frac{n\beta}{\alpha} + \frac{\beta}{\alpha}\sum_{i=1}^{n}(\frac{x_i}{\alpha})^\beta \\
&= - \frac{n\beta}{\alpha} + \frac{\beta}{\alpha}\sum_{i=1}^{n}(\frac{x_i}{\alpha})^\beta -⑤\\
\end{aligned} \\
（式②'をβで偏微分する）
\begin{aligned}
\frac{\partial \log L(\alpha,\beta)}{\partial \beta} &= \frac{n}{\beta} + \sum_{i=1}^{n}\log (\frac{x_i}{\alpha})- \sum_{i=1}^{n}(\frac{x_i}{\alpha})^\beta \log(\frac{x_i}{\alpha}) \\
&= \frac{n}{\beta} + \sum_{i=1}^{n}\log(x_i) - \sum_{i=1}^{n}\log(\alpha)-\frac{1}{\alpha^\beta}\sum_{i=1}^{n}(x_i)^\beta\log x_i+\frac{\log\alpha}{\alpha^\beta} \sum_{i=1}^{n}x_i^\beta \\
&= \frac{n}{\beta} + \sum_{i=1}^{n}\log(x_i) - n\log(\alpha)-\frac{1}{\alpha^\beta}\sum_{i=1}^{n}(x_i)^\beta\log x_i+\frac{\log\alpha}{\alpha^\beta} \sum_{i=1}^{n}x_i^\beta -⑥\\
\end{aligned} \\
式③、⑤より、
\begin{aligned}
-\frac{n\beta}{\alpha} + \frac{\beta}{\alpha}\sum_{i=1}^{n}(\frac{x_i}{\alpha})^\beta &= 0 \\
\frac{\beta}{\alpha^{\beta + 1}}\sum_{i=1}^{n}(x_i^\beta) &= \frac{n\beta}{\alpha} \\
\alpha \sum_{i=1}^{n}(x_i^\beta) &= n \alpha^{\beta + 1} \\
\sum_{i=1}^{n}(x_i^\beta) &= n \alpha^{\beta} \\
\alpha^{\beta} &= \frac{\sum_{i=1}^{n}(x_i^\beta)}{n} -⑦\\
\alpha &= (\frac{\sum_{i=1}^{n}(x_i^\beta)}{n})^{\frac{1}{\beta}} \\
\end{aligned} \\
式④、⑥より、
\begin{aligned}
\frac{n}{\beta} + \sum_{i=1}^{n}\log(x_i) - n\log(\alpha)-\frac{1}{\alpha^\beta}\sum_{i=1}^{n}(x_i)^\beta\log x_i+\frac{\log\alpha}{\alpha^\beta} \sum_{i=1}^{n}x_i^\beta &= 0 \\
-\frac{1}{\beta} - \frac{1}{n}\sum_{i=1}^{n}\log(x_i) + \log(\alpha)
+\frac{1}{n\alpha^\beta}\sum_{i=1}^{n}(x_i)^\beta\log x_i-\frac{\log\alpha}{n\alpha^\beta} \sum_{i=1}^{n}x_i^\beta &= 0 \\
\frac{1}{n\alpha^\beta}\sum_{i=1}^{n}(x_i)^\beta\log x_i-\frac{\log\alpha}{n\alpha^\beta} \sum_{i=1}^{n}x_i^\beta+ \log(\alpha) - \frac{1}{n}\sum_{i=1}^{n}\log(x_i)-\frac{1}{\beta} &= 0 -⑧\\
\end{aligned} \\
式⑦、⑧より、
\begin{aligned}
\frac{1}{n}\frac{1}{\alpha^\beta}\sum_{i=1}^{n}(x_i)^\beta\log x_i-\frac{\log\alpha}{n}\frac{1}{\alpha^\beta} \sum_{i=1}^{n}x_i^\beta+ \log(\alpha) - \frac{1}{n}\sum_{i=1}^{n}\log(x_i)-\frac{1}{\beta} &= 0 \\
\frac{1}{n}\frac{n}{\sum_{i=1}^{n}(x_i^\beta)}\sum_{i=1}^{n}(x_i)^\beta\log x_i-\frac{\log\alpha}{n}\frac{n}{\sum_{i=1}^{n}(x_i^\beta)} \sum_{i=1}^{n}x_i^\beta+ \log(\alpha) - \frac{1}{n}\sum_{i=1}^{n}\log(x_i)-\frac{1}{\beta} &= 0 \\
\frac{\sum_{i=1}^{n}(x_i)^\beta\log x_i}{\sum_{i=1}^{n}(x_i^\beta)}- \frac{1}{n}\sum_{i=1}^{n}\log(x_i)-\frac{1}{\beta} &= 0 \\
\beta &= 1 / (\frac{\sum_{i=1}^{n}(x_i)^\beta\log x_i}{\sum_{i=1}^{n}(x_i^\beta)}- \frac{1}{n}\sum_{i=1}^{n}\log(x_i) ) -⑨\\
\end{aligned} \\
まとめると
\begin{aligned}
\alpha &= (\frac{\sum_{i=1}^{n}(x_i^\beta)}{n})^{\frac{1}{\beta}} \\
\beta &= 1 / (\frac{\sum_{i=1}^{n}(x_i)^\beta\log x_i}{\sum_{i=1}^{n}(x_i^\beta)}- \frac{1}{n}\sum_{i=1}^{n}\log(x_i) ) \\
\end{aligned} \\
形状母数βが既知の場合、尺度母数αは計算できそうですね。
他方、形状母数βは両辺にβがいるので、このままだと少し難しいですね。
そこで、今回は以下のような方針でα、βの値を推定しようと思います。
- 1.ニュートン法を用いてβの値を推定する
- 2.推定したβを用いて、αの値を計算する
これから、ニュートン法を用いてβの値を求めるための準備をしましょう。
\begin{aligned}
\frac{\sum_{i=1}^{n}(x_i)^\beta\log x_i}{\sum_{i=1}^{n}(x_i^\beta)}- \frac{1}{n}\sum_{i=1}^{n}\log(x_i)-\frac{1}{\beta} &= 0 \\
\frac{1}{\beta} + \frac{1}{n}\sum_{i=1}^{n}\log(x_i) - \frac{\sum_{i=1}^{n}(x_i)^\beta\log x_i}{\sum_{i=1}^{n}(x_i^\beta)} &= 0 \\
\end{aligned} \\
f(β)を以下のように定義します。
\begin{aligned}
f(\beta) = \frac{1}{\beta} + \frac{1}{n}\sum_{i=1}^{n}\log(x_i) - \frac{\sum_{i=1}^{n}(x_i)^\beta\log x_i}{\sum_{i=1}^{n}(x_i^\beta)} \\
\end{aligned} \\
さらに、f(β)の微文形も用意します。
\begin{aligned}
\frac{\partial \sum_{i=1}^{n}(x_i^\beta)}{\partial \beta} &= \sum_{i=1}^{n}(x_i^\beta\log x_i) \\
\frac{\partial \sum_{i=1}^{n}(x_i)^\beta\log x_i}{\partial \beta} &= \sum_{i=1}^{n}(x_i^\beta\log^2 x_i) \\
\end{aligned} \\
\begin{aligned}
f'(\beta) &= -\frac{1}{\beta^2}-\frac{\sum_{i=1}^{n}x_i^\beta\log^2 x_i\sum_{i=1}^{n}(x_i^\beta)-(\sum_{i=1}^{n}(x_i^\beta\log x_i)^2}{(\sum_{i=1}^{n}(x_i^\beta))^2} \\
&=-\frac{1}{\beta^2}-\frac{\sum_{i=1}^{n}x_i^\beta\log^2x_i}{\sum_{i=1}^{n}(x_i^\beta)} + (\frac{\sum_{i=1}^{n}x_i^\beta\log x_i}{\sum_{i=1}^{n}(x_i^\beta)})^2
\end{aligned} \\
これでニュートン法によりβを推定する準備ができました。
実装編
この章では、２母数ワイブルの母数推定を行うためのプログラムを作成します。
使用する言語はpythonとしました。
作成するコードのイメージは以下の通りです。
- 検証用データの作成
- βの推定用
- αの推定用
検証用データの作成
この節では、指定したα、βの値を母数とする２母数ワイブル分布に従う乱数Xiを生成します。
先ほどの最尤法によるα、βの定義中のXiを作成するようなイメージです。
なお、検証用データと理論的なワイブル分布の概形を比較し、理論曲線に比べてどの程度一致するか確認します。
この節では以下のような関数を用意しました。
# ワイブル分布に従う乱数データの生成
def generate_weibull_samples(alpha,beta,n_sample):
xi = alpha * np.random.weibull(beta,n_sample)
return xi
# 理論的なワイブル分布を計算
def compute_weibull_distribution(alpha,beta,x_stop,step):
theoretical_x = np.arange(0,x_stop,step)
theoretical_y = (beta/alpha)*(theoretical_x/alpha)**(beta-1) * np.exp(-(theoretical_x/alpha)**beta)
return theoretical_x,theoretical_y
# 生成した乱数データの頻度分布とワイブル分布曲線の可視化
def plot_weibull_distribution(xi,theoretical_x,theoretical_y):
plt.hist(xi,bins=50, density=True, color='g', label="Sample Data")
plt.plot(theoretical_x, theoretical_y, color='r', label="Theoretical data")
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Weibull Distribution vs Sample Data')
plt.legend()
plt.show()
作成した関数を実行します。
alpha = 3.0
beta = 2.0
n_sample = 20000 # 生成する乱数のサンプル数
xi = generate_weibull_samples(alpha,beta,n_sample) # ワイブル分布に従う乱数の生成
x_stop = 10
step = 0.1
# 理論的なワイブル分布を計算
theoretical_values = compute_weibull_distribution(alpha,beta,x_stop,step)
# 生成した検証用データと理論曲線の可視化
plot_weibull_distribution(xi,theoretical_values[0],theoretical_values[1])
今回は、このような概形になりました。
検証用データとワイブル分布の理論曲線は概ね合っていますね。
したがって、α、βの推定値がそれぞれ3.0、2.0に近い値となればいいですね。
母数βの推定
ニュートン法によりβの推定を行います。そのため、先ほど示したf(β)、その微分形f'(β)の定義をプログラム中に記載しています。
# f(β)の定義
#############################################
def fb(beta,xi):
n = len(xi)
term1 = 1/beta
term2 = (1/n) * np.sum(np.log(xi))
term3 = - np.sum((xi**beta) * np.log(xi)) / np.sum((xi**beta))
return term1 + term2 + term3
#############################################
# f'(β)の定義
#############################################
def dfb(beta,xi):
term1 = -1/(beta**2)
term2 = -np.sum((xi**beta) * (np.log(xi))**2) / np.sum(((xi)**beta))
term3 = (np.sum((xi**beta) * (np.log(xi))) / np.sum(((xi)**beta) ))**2
return term1 + term2 + term3
#############################################
# ニュートン法の定義(母数βの推定用)
#############################################
def estimate_beta_newton(beta0,xi):
beta = beta0
box = []
for i in range(100):
if abs(fb(beta,xi)) < 1.0e-06:
break
else:
# ニュートン法の定義
beta = beta - fb(beta,xi) / dfb(beta,xi)
box.append(beta)
return beta,box
#############################################
estimate_beta_newton関数の戻り値は２つあり、「beta」は推測されたβの値。
boxは収束するまでの途中経過を保存した配列になります。
作成した関数を実行してみます。
alpha = 3.0
beta = 2.0
n_sample = 20000 # 生成する乱数のサンプル数
xi = generate_weibull_samples(alpha,beta,n_sample) # ワイブル分布に従う乱数の生成
# ニュートン法を使ってβの値を推定
beta0 = 4 # 初期値、推定開始値
beta = estimate_beta_newton(beta0,xi)[0]
print(beta) #今回は1.99908526070747
ある程度2.0に近い値が出力されました。
初期値の与え方も重要そうなので、この辺りは今後勉強します。
母数αの推定
前節までのコードによってβの推定ができました。これからαの推定用のプログラムを作成します。
# 母数αの推定用
#############################################
def estimate_alpha(beta,xi):
n = len(xi)
alpha = (np.sum(xi ** beta) / n) ** (1/beta)
return alpha
#############################################
alpha = estimate_alpha(beta,xi)
print("alpha:",alpha) #2.996716351398999
推定値
- α : 2.996716351398999
- β : 1.99908526070747
いい感じそう