【ワイブル分布】 未知パラメータ推定方法の改良

はじめに

今回はワイブル分布のパラメータ推定において最尤推定量による方法、確率紙による推定方法を改良した方法について説明します。
確率紙による推定方法、最尤推定による推定方法には２点問題点があります。
１点目は形状パラメータ $m$ の推定量は不偏推定量ではないため、期待値 $=$ 真値 とならない点です。
２点目は尺度パラメータ $\eta$ の推定量も不偏推定量ではなく、$m$ が小さいほど、期待値と真値との誤差が大きくなり、$m\rightarrow0$ で期待値が発散するという点です。
今回は、上記２点を改善した推定方法を発見しましたので、解説します。

概要

前置き

$X$ を形状パラメータ $m$ 、尺度パラメータ $\eta$ のワイブル分布 $W(m,\eta)$ に従う確率変数とすると、
確率密度関数 $f(x)$ 、累積分布関数 $F(x)$ は

$$\begin{aligned} f(x)&=\dfrac{m}{\eta}\left(\dfrac{x}{\eta} \right)^{m-1} \exp\left\{-\left(\dfrac{x}{\eta} \right)^{m} \right\} \\ F(x)&=1-\exp\left\{-\left(\dfrac{x}{\eta} \right)^{m} \right\} \\ \end{aligned}$$

と表されます。

$X \sim W(m,\eta)$ から互いに独立に得られた $n$ 個の観測値を $x_{1},\ ...\ ,x_{n}$ とします。

形状パラメータ m の改良版推定量

$m$ の最尤推定量を $\hat{m}_{2}$ として、改良後の推定量 $\hat{m}_{3}$ は、観測値の数 $n$ を用いて

$$\begin{aligned} \hat{m}_{3} = \dfrac{a_{m}n+b_{m}}{a_{m}n+b_{m}+1} \hat{m}_{2} \\ \left(a_{m} = 0.725703999,\ b_{m} = -1.361581962 \right) \end{aligned}$$

と表せます。$a_{m},b_{m}$ については数値解析により算出されるもので、上記の値が真の値と一致するとは断言できませんが、この数値を用いると良好な推定結果が得られます。

ここで、$m$ の最尤推定量 $\hat{m}_{2}$ は

$n$ 個の観測値 $x_{1},\ ...\ ,x_{n}$ を用いると、

$$\begin{aligned} g(\hat{m}_{2}) &= \dfrac{1}{\hat{m}_{2}} + \dfrac{1}{n}\sum_{i=1}^{n}\ln{x_{i}} -\dfrac{\sum_{i=1}^{n} x_{i}^{\hat{m}_{2}} \ln{x_{i}}}{\sum_{i=1}^{n} x_{i}^{\hat{m}_{2}}} \\ &=0 \end{aligned}$$

を満たします。解析的に解くことができないため、ニュートン法により、以下の漸化式を反復計算することにより $\hat{m}_{2}$ を計算します。

$$\begin{aligned} m_{k+1}= m_{k}- \dfrac{g(m_{k})}{g'(m_{k})} \end{aligned}$$

詳細は [1] をご覧ください

尺度パラメータ η の改良版推定量

$\eta$ の最尤推定量を $\hat{\eta}_{2}$ として、改良後の推定量 $\hat{\eta}_{3}$ は、観測値の数 $n$ 、ガンマ関数を用いて

$$\begin{aligned} \hat{\eta}_{3} &= n^{\frac{1}{\hat{m}_{3}}}\dfrac{\Gamma(n)} {\Gamma\left(n+\dfrac{1}{\hat{m}_{3}} \right)} \hat{\eta}_{2} \\ &= \dfrac{\hat{m}_{3}n^{\tfrac{1}{\hat{m}_{3}}}}{\Gamma\left( \dfrac{1}{\hat{m}_{3}} \right)} \left( \prod_{k=1}^{n-1} \dfrac{k}{ k+\dfrac{1}{\hat{m}_{3}} } \right) \hat{\eta}_{2} \end{aligned}$$

と表せます。ガンマ関数は発散しやすいため、$n$ が大きい場合は下段の式を用いた方が計算が発散しにくくなります。

$\eta$ の最尤推定量 $\hat{\eta}_{2}$ は

$$\begin{aligned} \hat{\eta}_{2} &= \sqrt[\footnotesize \hat{m}_{2}]{\dfrac{1}{n} \sum_{i=1}^{n} x_{i}^{\hat{m}_{2}}} \end{aligned}$$

と表せます。詳細は [2] をご覧ください。

各推定量の比較方法

従来と改良版の以下３方法による推定の精度を比較します。

• 方法１：確率紙を用いた推定 $\hat{m}_{1}$, $\hat{\eta}_{1}$
• 方法２：最尤推定を用いた推定 $\hat{m}_{2}, \hat{\eta}_{2}$
• 方法３：改良した最尤推定を用いた推定 $\hat{m}_{3}, \hat{\eta}_{3}$

方法の比較は以下の手順で実施します。

• 逆関数法によりワイブル分布に従う乱数を生成
• 得られた乱数（観測値）から、各方法において未知パラメータ $m,\eta$ を推定
• 上記の試行を１０００００回繰り返し、各方法について推定値の平均値と不偏分散を計算
• 平均値の相対誤差と、相対化した分散（不偏分散を真値の２乗で割った値）をグラフ化し比較

方法１、方法２の具体的な計算方法、ワイブル分布に従う乱数の生成についてはこちらの記事をご覧ください。

結果

• $m,\eta$ の両パラメータの推定において、改良版の方法３が最も相対誤差、分散が小さい
• $m$ の相対誤差は方法３ではほぼ 0 である　($\rightarrow$ 不偏推定量である）
• $\eta$ の相対誤差は 方法３でも $m$ が極端に小さく、$n$ も同時に小さい場合は大きくなるも３方法の中では最小。

もともと $m$ が小さいと観測値自体の分散が大きくなり、必然的に推定量の分散も大きくなります。いずれにしても $m$ が小さい場合は観測値の数を増やす必要があります。

詳細は [3] をご覧ください

Python ソースコード

今回の簡易シミュレーションで使用した Python のソースコードは [4] をご覧ください。

[1] 形状パラメータ m の改良版推定量の解説

まず、$m$ の最尤推定量 $\hat{m}_{2}$ の期待値 $E[\hat{m}_{2}]$ が $m,\eta,$ 乱数の数 $=$ 観測値の数 $n$ によってどのように表されるかを考えます。しかし、期待値は解析的には解くことができません。以下の手順で数値解析により算出します。

• ワイブル分布に従う乱数を生成
• 得られた乱数（観測値）から、最尤推定により $m$ を推定
• 上記の試行を１０００００回繰り返し、上で計算した推定値の平均をとり、期待値とする

$\hat{m}_{2}$ の期待値（$\approx$推定値の平均）と真値 $m$ の相対誤差 $\dfrac{E[\hat{m}_{2}]-m}{m}$ に着目すると、相対誤差と $m,\eta,n$ には以下の関係があることがわかりました

1. $\hat{m}_{2}$ の期待値の相対誤差は $m$ に依存しない
2. $\hat{m}_{2}$ の期待値の相対誤差は $\eta$ に依存しない
3. $\hat{m}_{2}$ の期待値の相対誤差の逆数と生成乱数の数をプロットすると直線になる
4. 3.で直線近似した傾き $a_{m}$ 、切片 $b_{m}$ は $m$ によって変化しない

1.については、横軸に真値 $m$ 、縦軸に $m$ の最尤推定値の平均 ($\approx E[\hat{m}_{2}]$ ) の相対誤差をとったグラフを示します。
条件：$m=0.1,0.2,\ ...\ ,5.0$ 、$\eta=100$ 、生成乱数の数 $n=10$ 、試行回数 $N=100000$

図１：$m$ を変化させたときの $m$ の最尤推定値の平均値（$\approx$ 期待値）の相対誤差の推移

上図１より、相対誤差は $m$ に関係なく一定であると考えられます。

2.については、横軸に真値 $\eta$ 、縦軸に最尤推定値の平均 ($\approx E[\hat{m}_{2}]$ ) の相対誤差をとったグラフを示します。
条件：$m=0.5$ 、$\eta=1,2,5,10,\ ...\ ,10000$ 、生成乱数の数 $n=10$ 、試行回数 $N=100000$

図２：$\eta$ を変化させたときの $m$ の最尤推定値の平均値（$\approx$ 期待値）の相対誤差の推移

上図２より、相対誤差は $\eta$ に関係なく一定であると考えられます。

3.については、横軸に生成乱数の数 $n$ 、縦軸に最尤推定値の平均 ($\approx E[\hat{m}_{2}]$ ) の相対誤差の逆数をとったグラフを示します。
条件：$m=0.5$ 、$\eta=100$ 、乱数の数 $n=4,5,\ ...\ ,50$ 、試行回数 $N=100000$

図３－１：生成乱数の数 $n$ を変化させたときの $m$ の最尤推定値の平均値（$\approx$ 期待値）の相対誤差の逆数の推移

上図３ー１より、相対誤差の逆数と乱数の数をプロットすると非常な綺麗な直線となります。 $n$ が大きくなると、近似直線からのズレが大きくなっています。以下は、生成乱数の数 $n$ を $10$ から $500$ まで変化させたときの相対誤差の逆数のグラフになります。
条件：$m=0.5$ 、$\eta=100$ 、乱数の数 $n=10,20,\ ...\ ,500$ 、試行回数 $N=100000$

図３－２：生成乱数の数 $n$ を変化させたときの $m$ の最尤推定値の平均値（$\approx$ 期待値）の相対誤差の逆数の推移 $n=10\sim500$

図３－２から、$n$ を大きくすればするほど近似直線からのズレが大きくなることが分かります。これは、$n$ が大きくなることにより、相対誤差が 0 に近づきます。従って、推定値の平均値のばらつき（ノイズ）の影響が相対的に大きくなったためと考えられます。その影響で、近似直線の計算には $n=50$ までのデータを使用しています。また、$m$ が大きく、$n$ が小さいほど、$m$ の最尤推定の計算が発散しやすくなります。$n=3$ の場合は、発散するケースが頻発するため、$n\geq4$ としています。

4.について、相対誤差と乱数の数 $n$ を直線近似し、相対誤差の逆数 $=a_{m}n +b_{m}$ とします。
このとき、$a_{m},b_{m}$ は $m$ に関係なく一定の値と考えられます。

横軸に真値 $m$ 、縦軸に $a_{m},b_{m}$ をとったグラフを以下に示します。
条件：$m=0.1,0.2,\ ...\ ,5.0$ 、$\eta=100$ 、生成乱数の数 $n=4,5,\ ...\ ,50$ 、試行回数 $N=100000$

図４：$m$ を変化させたときの傾き $a_{m}$ 切片 $b_{m}$ の推移

以上より、相対誤差の逆数が $m,\eta$ に関係なく一定で、$n$ について1次式で表せると考えると、真値 $m$ と 推定量 $\hat{m}_{2}$ の期待値の間に以下の関係があると考えられます。

$$\begin{aligned} \dfrac{m}{E[\hat{m}_{2}] -m} = a_{m} n + b_{m} \end{aligned}$$

これを、式変形することにより、

$$\begin{aligned} E[\hat{m}_{2}] -m &= \dfrac{m}{a_{m} n + b_{m}} \\ E[\hat{m}_{2}] &= \dfrac{a_{m} n + b_{m} +1}{a_{m} n + b_{m}} m \end{aligned}$$

と表すことができます。したがって、$\hat{m}_{3} = \dfrac{a_{m} n + b_{m} }{a_{m} n + b_{m}+1} \hat{m}_{2}$ とすれば、$\hat{m}_{3}$ の期待値は $m$ となり、近似的に不偏推定量として扱うことができます。

$a_{m},b_{m}$ については 4. の結果を利用し平均化すると、

$$\begin{aligned} a_{m} = 0.725703999, b_{m} = -1.361581962 \end{aligned}$$

となります。

[2] 尺度パラメータ η の改良版推定量の解説

まず、$\eta$ の最尤推定量 $\hat{\eta}_{2}$ の期待値が、$m,\eta$ の真値によってどのように変化するかを確認します。

縦軸に $\hat{\eta}_{2}$ の平均値 ($\approx$ 期待値)の相対誤差、横軸に $m$ の真値をとると、以下のようにグラフになります。
条件：$m=0.2,\ ...\ ,5.0$ 、$\eta=100$ 、生成乱数の数 $n=10$ 、試行回数 $N=100000$

図５：$m$ を変化させたときの $\eta$ の推定値の平均値（ $\approx$ 期待値）の相対誤差の推移

$m \gt 1$ の場合は、誤差は比較的少ないですが、$m$ が小さくなるほど、誤差が大きくなり、$m$ が 0 に近づくにつれ無限大に発散します。

次に、推定値 $\hat{\eta}_{2}$ が、$\eta$ の真値によってどのように変化するかを確認します。
縦軸に $\hat{\eta}_{2}$ の平均値 ($\approx$ 期待値) の相対誤差、横軸に $\eta$ の真値をとると、以下のようなグラフになります。
条件：$m=0.5$ 、$\eta=1,2,5,10,\ ...\ ,10000$ 、生成乱数の数 $n=10$ 、試行回数 $N=100000$

図６：$\eta$ を変化させたときの $\eta$ の推定値の平均値（ $\approx$ 期待値）の相対誤差の推移

上図６より、相対誤差は $\eta$ に依存しないと考えられます。

$\eta$ の最尤推定量 $\hat{\eta}_{2}$ は、不偏推定量ではなく、期待値は $m$ に依存することは分かりますが、グラフをプロットしても関係性は見えてきません。
そこで、$\hat{\eta}_{2}$ の期待値は解析的には求められませんが、形状パラメータ $m$ が既知の場合は、 期待値を求めることができます。 $\eta$ の最尤推定量 $\hat{\eta}$ は以下の式で表せます。

$$\begin{aligned} \hat{\eta} &= \sqrt[\footnotesize m]{\dfrac{1}{n} \sum_{i=1}^{n} X_{i}^{m}} \end{aligned}$$

m が既知の場合の η の最尤推定量の期待値

形状パラメータの真値 $m$ が既知である場合の $\eta$ の最尤推定量 $\hat{\eta}$ の期待値を計算します。

まず、 $X_{i} \rightarrow V_{i} =X_{i}^{m}$ に変換すると、$V_{i}$ の確率密度関数 $f_{V_{i}}(v)$ は

$$\begin{aligned} f_{V_{i}}(v) &= f_{X_{i}}(x) \left| \dfrac{dx}{dv} \right| \\ &= \dfrac{1}{\eta^{m}} mx^{m-1} \exp{\left[ -\dfrac{v}{\eta^{m}}\right]} \left| \dfrac{1}{\dfrac{dv}{dx}} \right| \\ &=\dfrac{1}{\eta^{m}} \exp{\left[ -\dfrac{v}{\eta^{m}}\right]}\\ \end{aligned}$$

となり、$V_{i} \sim Exp(\eta^{-m})$ が導かれます。

次に、 $V \rightarrow W_{n} = \sum_{i=1}^{n} V_{i}$ に変換すると、 $W_{n} \sim Gamma(n,\eta^{-m})$ となることから、
$W_{n}$ の確率密度関数 $f_{W_{n}}(w)$ は

$$\begin{aligned} f_{W_{n}}(w) &= \dfrac{w^{n-1}}{\Gamma(n) (\eta^{m})^{n}} \exp{\left[ -\dfrac{y}{\eta^{m}} \right]} \\ \end{aligned}$$

と表せます。最尤推定量 $\hat{\eta}$ の期待値を計算すると、

$$\begin{aligned} E[\hat{\eta}] &= E \left[ \left(\dfrac{1}{n}W_{n} \right)^{\tfrac{1}{m}} \right] \\ &= \int_{0}^{\infty} n^{-\tfrac{1}{m}} w^{\frac{1}{m}} \dfrac{w^{n-1}}{\Gamma(n) (\eta^{m})^{n}} \exp{\left[ -\dfrac{y}{\eta^{m}} \right]} dw\\ &= n^{-\tfrac{1}{m}} \dfrac{\Gamma\left(n+\dfrac{1}{m} \right)}{\Gamma(n)} \eta^{m \cdot \tfrac{1}{m}}\int_{0}^{\infty} \dfrac{w^{n+\frac{1}{m}}}{\Gamma \left(n+\dfrac{1}{m} \right) (\eta^{m})^{n+\frac{1}{m}}} \exp{\left[ -\dfrac{y}{\eta^{m}} \right]} dw \rightarrow Gamma\left(n+\dfrac{1}{m},\eta^{-m} \right) \\ &= n^{-\tfrac{1}{m}} \dfrac{\Gamma\left(n+\dfrac{1}{m} \right)}{\Gamma(n)} \eta \\ \end{aligned}$$

となります。以上より、補正係数を $C_{\eta}=n^{\tfrac{1}{m}} \dfrac{\Gamma(n)}{\Gamma\left(n+\dfrac{1}{m} \right)}$ とすれば、$C_{\eta} \hat{\eta}$ は $\eta$ の不偏推定量となります。また、補正係数 $C_{\eta}$ は以下のように変形できます。

$$\begin{aligned} C_{\eta} &= \dfrac{mn^{\tfrac{1}{m}}}{\Gamma\left( \dfrac{1}{m} \right)} \left( \prod_{k=1}^{n-1} \dfrac{k}{ k+\dfrac{1}{m} } \right) \end{aligned}$$

ガンマ関数は値が大きくなりやすいため、$n$ が大きい場合や $m$ が小さいは、上式の形にして計算した方が、計算が発散しにくくなります。

m が未知の場合の η の推定と各方法の比較

形状パラメータ $m$ が未知の場合は、$\hat{\eta}$ の推定量について真値 $m$ は使用できないため、最尤推定量 $\hat{m}_{2}$ または、改良した推定量 $\hat{m}_{3}$ を代用します。選択肢は、最尤推定量の計算に使う $m$ について $\hat{m}_{2},\ \hat{m}_{3}$ の２通り、補正係数に使う $m$ について $\hat{m}_{2},\ \hat{m}_{3}$ の２通りで、以下の４通りが考えられます。

1. $\hat{\eta}_{22}$ : 最尤推定の計算 $\rightarrow\hat{m}_{2}$ 、補正係数 $C_{\eta}$ の計算 $\rightarrow\hat{m}_{2}$
2. $\hat{\eta}_{33}$ : 最尤推定の計算 $\rightarrow\hat{m}_{3}$ 、補正係数 $C_{\eta}$ の計算 $\rightarrow\hat{m}_{3}$
3. $\hat{\eta}_{23}$ : 最尤推定の計算 $\rightarrow\hat{m}_{2}$ 、補正係数 $C_{\eta}$ の計算 $\rightarrow\hat{m}_{3}$
4. $\hat{\eta}_{32}$ : 最尤推定の計算 $\rightarrow\hat{m}_{3}$ 、補正係数 $C_{\eta}$ の計算 $\rightarrow\hat{m}_{2}$

上の４通りについて推定の精度を比較します。比較は以下の手順で実施します。

• ワイブル分布に従う乱数を生成
• 得られた乱数（観測値）から、４方法において尺度パラメータ $\eta$ を推定
• 上記の試行を１０００００回繰り返し、各推定値の平均値と不偏分散を計算
• 平均値の相対誤差と、相対化した分散（不偏分散を真値 $\eta$ の２乗で割った値）をグラフ化し比較

以下、上記の４通りについて、推定値の平均 ( $\approx$ 期待値) の相対誤差、分散のグラフを示します。

最初に、$m$ を $0.2$ から $5$ まで変化させたときの推定値の平均の相対誤差、相対化した分散をグラフ化したものを以下に示します。
条件：$m=0.2,\ ...\ ,5.0$ 、$\eta=100$ 、生成乱数の数 $n=10$ 、試行回数 $N=100000$

図７：$m$ を変化させたときの各方法の相対誤差と分散の推移

次に、$\eta$ を $1$ から $10000$ まで変化させたときの推定値の平均の相対誤差、相対化した分散をグラフ化したものを以下に示します。
条件：$m=0.5$ 、$\eta=1,2,5,10,\ ...\ ,10000$ 、生成乱数の数 $n=10$ 、試行回数 $N=100000$

図８：$\eta$ を変化させたときの各方法の相対誤差と分散の推移

最後に、生成乱数の数　$n$ を $5$ から $10000$ まで変化させたときの推定値の平均の相対誤差、相対化した分散をグラフ化したものを以下に示します。
条件：$m=0.5$ 、$\eta=100$ 、生成乱数の数 $n=5,10,20,50,\ ...\ ,10000$ 、試行回数 $N=100000$

図９：生成乱数の数 $n$ を変化させたときの各方法の相対誤差と分散の推移

以上より、$\hat{\eta}_{23}$ つまり、最尤推定の計算には $\hat{m}_{2}$ 、補正係数には $\hat{m}_{3}$ を使用すると、最も誤差が少なくなります。

[3] 各推定量の比較結果

前置き

本項では、pythonを使用したシミュレーションにより、以下の３方法による推定量の推定精度を比較します。

• 方法１：確率紙を用いた推定量 $\hat{m}_{1}$, $\hat{\eta}_{1}$
• 方法２：最尤推定量を用いた推定 $\hat{m}_{2}, \hat{\eta}_{2}$
• 方法３：改良した推定 $\hat{m}_{3}, \hat{\eta}_{3}$

比較方法は以下の手順で実施します。

• 逆関数法によりワイブル分布に従う乱数を生成
• 得られた乱数（観測値）から、各方法において未知パラメータ $m,\eta$ を推定
• 上記の試行を100000回繰り返し、各推定値の平均値と不偏分散を計算
• 平均値の相対誤差と、相対化した分散（不偏分散を真値の２乗で割った値）をグラフ化

推定値の評価には平均値と真値の相対誤差と分散を真値の二乗で割った、相対化した分散を使用します。
未知パラメータの真値を $\theta$、$i$ 回目の推定値を $\hat{\theta}_{i}$、推定値の平均を $\hat{\theta}_{mean}$、不偏分散を $S^{2}_{\hat{\theta}}$、試行回数を $N$ として、

• 推定値の相対誤差 $=\dfrac{\hat{\theta}_{mean}-\theta}{\theta}$
• 相対化した分散 $=\dfrac{S^{2}_{\hat{\theta}}}{\theta^{2}}$

$$\begin{aligned} \hat{\theta}_{mean} = \dfrac{1}{N}\sum_{i=1}^{N} \hat{\theta}_{i} ,\ S^{2}_{\hat{\theta}} = \dfrac{1}{N-1}\sum_{i=1}^{N} (\hat{\theta}_{i} - \hat{\theta}_{mean})^{2} \end{aligned}$$

と表せます。

結果のグラフ

最初に、$m$ を $0.1$ から $5$ まで変化させたときの推定値の平均値の相対誤差、相対化した分散をグラフ化したものを以下に示します。
条件：$m=0.1,\ ...\ ,5.0$ 、$\eta=100$ 、生成乱数の数 $n=10$ 、試行回数 $N=100000$

図１０－１：$m$ を変化させたときの各方法の相対誤差と分散の推移

$m$ が小さい場合は $\eta$ の推定値は非常に大きくなります。$m =0.1,\ ...\ ,0.5$ の場合は別途以下のグラフに示します。

図１０－２：$m$ を変化させたときの各方法の相対誤差と分散の推移 $(0 \lt m \lt 0.5)$

次に、$\eta$ を $1$ から $10000$ まで変化させたときの推定値の平均値の相対誤差、相対化した分散をグラフ化したものを以下に示します。
条件：$m=0.5$ 、$\eta=1,2,5,10,\ ...\ ,10000$ 、生成乱数の数 $n=10$ 、試行回数 $N=100000$

図１１：$\eta$ を変化させたときの各方法の相対誤差と分散の推移

最後に、生成乱数の数 $n$ を $5$ から $10000$ まで変化させたときの推定値の平均値の相対誤差、相対化した分散をグラフ化したものを以下に示します。
条件：$m=0.5$ 、$\eta=100$ 、生成乱数の数 $n=5,10,20,50,\ ...\ ,10000$ 、試行回数 $N=100000$

図１２：生成乱数の数 $n$ を変化させたときの各方法の相対誤差と分散の推移

まとめ

以上の結果から、$m,\eta$ の両推定については改良した推定方法３が最も相対誤差が少なく、分散が小さい結果となります。

[4] Python ソースコード

以下、簡易シミュレーションで使用したPythonのソースコードです。
numpy,pandas は標準ライブラリにありませんので、インストールが必要になります。
バージョン情報
python : 3.12.0
numpy : 1.26.2
pandas : 2.1.3

import numpy as np
import math
import pandas as pd
import warnings
import os


# ηの推定量の補正係数の計算
def CALC_CETA(m,num):

    c= 1.0
    for k in range(1,num):
        c *= k/(k + 1/m)

    c *= m * (num ** (1/m)) / math.gamma(1/m) 
    #c = num ** (1/m) * math.gamma(num) / math.gamma(num+1/m)

    return c

# num個の乱数によるパラメータの推定を　num_of_times 回実施した計算結果を出力する
def CALC(m,eta,num,num_of_times):

    # 計算結果の配列の定義
    m_hat1_arr = np.empty(0)    #確率紙による m
    eta_hat1_arr= np.empty(0)   #確率紙によるη

    m_hat2_arr = np.empty(0)    #最尤推定による m
    eta_hat2_arr = np.empty(0)  #最尤推定によるη

    m_hat3_arr = np.empty(0)    #最尤推定の改良版による m
    eta_hat3_arr= np.empty(0)   #最尤推定の改良版によるη

    # mの最尤推定量補正用 a_m, b_m
    a_m = 0.725703999
    b_m = -1.361581962
 
    # 累積分布関数F(x)の生成 #累積分布関数は乱数に依存しないため、ループの外で実施し、計算量を減らす
    F = np.arange(1,num+1)
    F= (F-0.3)/(num +0.4) #メジアンランク　Benardの近似式

    # y = ln {-ln(1-F(x)) }　の計算
    y=np.log(-np.log(1-F))
    y_bar = np.average(y)

    # num個の乱数によるパラメータの推定を　num_of_times 回実施
    for i in range(num_of_times):

        # 最尤推定の m が発散 -> 乱数再生成のループ
        while(1):

            # Retryフラグ
            retry_flag = 0
 
            # 乱数生成用
            rng = np.random.default_rng()

            # U(0,1)の一様乱数を生成
            u = rng.random(num)

            # 昇順に並び替え
            u = np.sort(u)

            # 逆関数法によりワイブル分布に従う乱数の生成
            x = eta * (-np.log(1-u)) **(1/m)

            logx = np.log(x)
            logx_bar = np.average(logx)

            # y = a * logx + b  と線形近似し、最小二乗法による傾きaと切片bを計算

            a = (np.dot(logx-logx_bar,y-y_bar) / np.sum((logx-logx_bar) ** 2))
            b = y_bar - logx_bar * a

            # a,bからパラメータm,ηの推定 -> 確率紙による推定
            m_hat1 = a
            eta_hat1 = np.exp(-b/a) 


            # ニュートン法による最尤推定のmの算出　L(m,η)の m 微分が 許容値以下になるまで反復計算
            error = 1 #L(m,η)の m 微分の値初期値
            rep=0 #反復回数
            error_limit = 1e-8 #許容値
            m_hat = m_hat1 # m の初期値 m0 は 確率紙による推定を使用

            
            while (abs(error) >= error_limit):

                try:
                    xm = x ** m_hat
                    xmlogx_ov_xm = np.sum(np.dot(xm,logx)) / np.sum(xm)

                    # g(m), g'(m) の計算
                    g_m = 1/m_hat + logx_bar - xmlogx_ov_xm 
                    g_m_diff = - m_hat ** (-2) - np.sum(np.dot(xm,logx ** 2))/np.sum(xm)  + xmlogx_ov_xm ** 2

                    #対数尤度関数の m 微分　d/dm lnL(m,η) = num * g(m)の計算
                    error = num * g_m

                    # mの更新
                    m_hat = m_hat - g_m / g_m_diff

                # m の発散を検知 -> retry フラグ =1
                except Exception as e:
                    #print(e)
                    retry_flag = 1
                    break

                # m の反復計算が100回を超える場合は計算方法を変更
                if 100 < rep <=200 :
                    m_hat = 1/2 * (1/( xmlogx_ov_xm  -logx_bar) + m_hat )
                    #print(f'{i}回目：反復回数{rep}、m^={m_hat}、g(m)={g_m}、g\'(m)={g_m_diff}')

                # 200回超えても計算が収束しない場合は乱数の生成からやり直し
                elif rep > 200:
                    print(f'mの最尤推定の計算が収束しません。乱数を再生成します。')
                    break
                    
                rep += 1

            # m の発散を検知 -> retryフラグ =1 -> 乱数再生成
            # m 収束　　　　 -> retryフラグ=0 ->　続行 
            if retry_flag == 1:
                print(f'{i}番目の試行:mの値{m_hat}が発散しました。乱数を再生成します。')
                continue
            else:
                break
       
        m_hat2 = m_hat

        # m の補正
        
        m_hat3 = m_hat * (a_m* num + b_m)/ (a_m* num + b_m +1)

        # ηの計算
        eta_hat2 = np.average(x ** m_hat2) ** (1/m_hat2)
        eta_hat3 = eta_hat2 * CALC_CETA(m_hat3,num)

        # 各方式の推定の計算結果を配列に追加        
        m_hat1_arr = np.append(m_hat1_arr,m_hat1)
        eta_hat1_arr = np.append(eta_hat1_arr,eta_hat1)

        m_hat2_arr = np.append(m_hat2_arr,m_hat2)
        eta_hat2_arr = np.append(eta_hat2_arr,eta_hat2)
        
        m_hat3_arr = np.append(m_hat3_arr,m_hat3)
        eta_hat3_arr = np.append(eta_hat3_arr,eta_hat3)

    # 平均、不偏分散を計算
    m_hat1_avg = np.average(m_hat1_arr)
    m_hat1_var = np.var(m_hat1_arr,ddof=1)

    eta_hat1_avg = np.average(eta_hat1_arr)
    eta_hat1_var = np.var(eta_hat1_arr,ddof=1)

    m_hat2_avg = np.average(m_hat2_arr)
    m_hat2_var = np.var(m_hat2_arr,ddof=1)

    eta_hat2_avg = np.average(eta_hat2_arr)
    eta_hat2_var = np.var(eta_hat2_arr,ddof=1)

    m_hat3_avg = np.average(m_hat3_arr)
    m_hat3_var = np.var(m_hat3_arr,ddof=1)

    eta_hat3_avg = np.average(eta_hat3_arr)
    eta_hat3_var = np.var(eta_hat3_arr,ddof=1)

    #結果出力用配列の作成
    outdata = np.array([m, eta, num, num_of_times,
                        m_hat1_avg, m_hat1_var, eta_hat1_avg, eta_hat1_var,
                        m_hat2_avg, m_hat2_var, eta_hat2_avg, eta_hat2_var, 
                        m_hat3_avg, m_hat3_var, eta_hat3_avg, eta_hat3_var
                        ])

    #結果の二次元配列化　（列を結合できるように）
    outdata = outdata.reshape(len(outdata),1)

    #結果表示
    print(f'確率紙によるmの推定値平均、分散:{m_hat1_avg}、{m_hat1_var}')
    print(f'確率紙によるηの推定値平均、分散:{eta_hat1_avg}、{eta_hat1_var}')
    print(f'最尤推定によるmの推定値平均、分散:{m_hat2_avg}、{m_hat2_var}')
    print(f'最尤推定によるηの推定値平均、分散:{eta_hat2_avg}、{eta_hat2_var}')
    print(f'最尤推定の改良版によるmの推定値平均、分散:{m_hat3_avg}、{m_hat3_var}')
    print(f'最尤推定の改良版によるηの推定値平均、分散:{eta_hat3_avg}、{eta_hat3_var}')

    return outdata
    
def main():

    ## numpy の 計算オーバーフローのwarning を error としてtry except を作動させ、乱数の生成までやり直し
    # m の発散が頻発する場合は条件を変更した方がよい。(ｍが大きい、nが小さいほど発散しやすい)
    warnings.resetwarnings()
    warnings.simplefilter('error')
    
    #出力結果格納用配列
    dataset = np.empty((16,0))

    # パラメータの設定
    m = 0.5 #形状パラメータ
    eta = 100 #尺度パラメータ
    num = 10 #生成乱数の個数 
    num_of_times = 100000 #試行回数 
    No = 1

    m_start, m_stop, m_step = 0.1, 5.0, 0.1
    num_start, num_stop,num_step = 4, 50, 1

    #for m in np.arange(m_start,m_stop+m_step,m_step):
    #for eta in np.array((1,2,5,10,20,50,100,200,500,1000,2000,5000,10000)):
    for num in np.arange(num_start,num_stop+1,num_step):
    #for num in np.array((5,10,20,50,100,200,500,1000,2000,5000,10000)):


        #計算条件の表示
        print(f'No.{No}:  m={m},  η={eta},  乱数の数={num},  試行回数={num_of_times} の条件で計算中')

        #各パラメータの計算結果の列を結合
        dataset = np.append(dataset,CALC(m,eta,num,num_of_times),axis=1)
        No +=1

    # ファイル出力
    ## 表題作成
    index1 = ['mの真値', 'ηの真値', '生成乱数の個数', '試行回数',
              '確率紙:mの推定値の平均', '確率紙:mの推定値の分散', '確率紙:ηの推定値の平均','確率紙:ηの推定値の分散',
              '最尤推定:mの推定値の平均', '最尤推定:mの推定値の分散', '最尤推定:ηの推定値の平均','最尤推定:ηの推定値の分散', 
              '最尤推定改良版:mの推定値の平均','最尤推定改良版:mの推定値の分散','最尤推定改良版:ηの推定値の平均','最尤推定改良版:ηの推定値の分散'
              ]
    
    header1 = []
    for i in range(dataset.shape[1]):
        header1.append(f'No.{i+1}')

    # データフレーム形式に変換
    pd_out = pd.DataFrame(dataset,columns=header1, index=index1)

    #ファイル名の設定
    outname0 = f'm={m}_η={eta}_乱数の数={num}_試行回数={num_of_times}'
    outname = outname0 + '.csv'


    #ファイル名が重複する場合、ファイル名を自動変更
    filenum = 1
    while( os.path.exists(outname) ):
        outname = outname0 + f'_{filenum}.csv'
        filenum += 1

    # csvデータ出力　#shift-jis でないとエクセルで文字化け
    pd.DataFrame(pd_out).to_csv(outname, encoding='shift-jis')  

    input('計算が完了しました。ENTERで最初からやり直します。')

if __name__ == "__main__":

    while(1):
        main()