最尤法による２母数ワイブル分布の母数推定 [理論・実装(python)]

# はじめに
pythonの勉強中なので、少しプログラムを書いてみました。
途中、少し、計算過程を記載していますが、、温かい目でお願いします。

# 目的・テーマ
以下に示す２母数ワイブル分布の母数推定について勉強してみました。
ここで、α、βをそれぞれ尺度母数、形状母数です。
分布の形状を決めるパラメータですね。
```math
f(x;\alpha,\beta) = \frac{\beta}{\alpha} \left( \frac{x}{\alpha} \right)^{\beta - 1} \exp \left[ - \left( \frac{x}{\alpha} \right)^\beta \right]
```
今回は最尤法を用いて２母数ワイブル分布の母数推定にチャレンジしようと思います。


# 理論編
２母数ワイブル分布が上記の式で表される時、尤度関数は以下のようになります。
```math
\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}
```
対数尤度関数は以下のようになります。

```math
\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]  -②
```
まずは式②の右辺をわかりやすくしましょう。

```math
\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} \\
```
以下の式を用いて、α、βの値を求めます。
```math
\begin{aligned}
\frac{\partial \log L(\alpha,\beta)}{\partial \alpha} &= 0 -③ \\
\frac{\partial \log L(\alpha,\beta)}{\partial \beta} &= 0  -④
\end{aligned}
```

式②'を踏まえて、式③、④の左辺を計算してみましょう。
（式②'をαで偏微分する）
```math
\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} \\
```
（式②'をβで偏微分する）
```math
\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} \\
```

式③、⑤より、
```math
\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} \\
```

式④、⑥より、

```math
\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} \\
```

式⑦、⑧より、
```math
\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} \\
```
まとめると
```math
\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.推定したβを用いて、αの値を計算する

これから、ニュートン法を用いてβの値を求めるための準備をしましょう。

```math
\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(β)を以下のように定義します。
```math
\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(β)の微文形も用意します。
```math
\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} \\
```

```math
\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を作成するようなイメージです。

なお、検証用データと理論的なワイブル分布の概形を比較し、理論曲線に比べてどの程度一致するか確認します。

この節では以下のような関数を用意しました。
```python
# ワイブル分布に従う乱数データの生成
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()
```
作成した関数を実行します。
```python
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に近い値となればいいですね。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3831763/3c0738a9-8eed-c0bb-e3a9-d61511bb28e8.png)

## 母数βの推定

ニュートン法によりβの推定を行います。そのため、先ほど示したf(β)、その微分形f'(β)の定義をプログラム中に記載しています。
```python
# 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は収束するまでの途中経過を保存した配列になります。
作成した関数を実行してみます。
```python
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に近い値が出力されました。
初期値の与え方も重要そうなので、この辺りは今後勉強します。
## 母数αの推定

前節までのコードによってβの推定ができました。これからαの推定用のプログラムを作成します。
```python
# 母数αの推定用
#############################################
def estimate_alpha(beta,xi):
  n = len(xi)
  alpha = (np.sum(xi ** beta) / n) ** (1/beta)
  return alpha
#############################################
```

```python
alpha = estimate_alpha(beta,xi)
print("alpha:",alpha) #2.996716351398999
```
推定値
- α : 2.996716351398999
- β : 1.99908526070747

いい感じそう