【ワイブル分布】確率紙による推定と最尤推定によるパラメータ推定方法の精度比較

はじめに

今回はワイブル分布のパラメータ推定において、確率紙による方法１と最尤推定を用いた方法２の比較をします。
結論としては、形状パラメータ $mm$ については確率紙による推定方法１が優れ、尺度パラメータ $η\eta$ は両方法で大きな違いはありませんが、若干最尤推定による推定方法２の方が優れているという結果になります。
ただ、両方法についてワイブル分布の形状パラメータ $mm$ が 0 に近づくにつれ、 尺度パラメータ $η\eta$ の推定量は正の無限大に発散するため、推定精度は低下します。
最尤推定による推定を改良した場合は両推定方法に勝ります。詳細はこちらの記事をご覧ください。

概要

前置き

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

$$f(x)=mη(xη)m−1exp⁡{−(xη)m}F(x)=1−exp⁡{−(xη)m}\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∼W(m,η)X \sim W(m,\eta)$ から互いに独立に得られた $nn$ 個の観測値を小さい順に $x(1), ... ,x(n)x_{(1)},\ ...\ ,x_{(n)}$ とします。

方法１ 確率紙による推定

$XX$ の累積分布関数 $F(x)F(x)$ と観測値 $xx$ について $y=ln⁡[−ln⁡{1−F(x)}]y=\ln{[-\ln{\{1-F(x)\}}]}$ と $ln⁡x\ln{x}$ は直線になります。
これを利用し、観測値から $ln⁡x,y\ln{x},y$ をプロットし、最小二乗法により未知パラメータを求める方法です。
得られた観測値 $x(i)x_{(i)}$ に対し、$yi=ln⁡[−ln⁡{1−F(x(i))}](i=1,... ,n)y_{i} = \ln{\left[- \ln{\{ 1-F(x_{(i)})\}} \right]} (i=1, ...\ ,n)$ とすると、
$m,ηm,\eta$ の推定値 $m^1,η^1\hat{m}_{1}, \hat{\eta}_{1}$ は

$${m^1=∑i=1n(ln⁡x(i)−ln⁡x‾)(yi−yˉ)∑i=1n(ln⁡x(i)−ln⁡x‾)2η^1=exp⁡[−yˉ−m^1ln⁡x‾m^1]\left\{\begin{aligned} \hat{m}_1 &= \dfrac{\sum_{i=1}^{n}(\ln{x_{(i)}} - \overline{\ln{x}})(y_{i}-\bar{y})} {\sum_{i=1}^{n}(\ln{x_{(i)}}- \overline{\ln{x}})^{2}} \\ \hat{\eta}_1 &= \exp{\left[ -\dfrac{\bar{y}-\hat{m}_1\overline{\ln{x}}}{\hat{m}_1} \right]} \end{aligned}\right.$$

$$(yˉ=1n∑i=1nyi , ln⁡x‾=1n∑i=1nln⁡x(i))\left( \begin{aligned} \bar{y}= \dfrac{1}{n}\sum_{i=1}^{n} y_{i}\ , \ \overline{\ln{x}} =\dfrac{1}{n} \sum_{i=1}^{n} \ln{x_{(i)}} \end{aligned} \right)$$

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

方法２ 最尤推定による推定

$X∼W(m,η)X \sim W(m,\eta)$ から互いに独立に得られた $nn$ 個の観測値を $x1, ... ,xnx_{1},\ ...\ ,x_{n}$ とすると、

$mm$ の最尤推定量 $m^2\hat{m}_{2}$ は

$$g(m)=1m+1n∑i=1nln⁡xi−∑i=1nximln⁡xi∑i=1nxim=0\begin{aligned} g(m) &= \dfrac{1}{m} + \dfrac{1}{n}\sum_{i=1}^{n}\ln{x_{i}} -\dfrac{\sum_{i=1}^{n} x_{i}^{m} \ln{x_{i}}}{\sum_{i=1}^{n} x_{i}^{m}} \\ &=0 \end{aligned}$$

を満たす $mm$ となります。
ニュートン法により、以下の漸化式を反復計算することにより $m^2\hat{m}_{2}$ を計算します。

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

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

$$η^2=1n∑i=1nxim^2m^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}$$

と表せます。詳細は [3] をご覧ください。
欠点として、$mm$ が大きく、$nn$ が小さいほど、$mm$ の最尤推定の計算が発散しやすくなります。参考程度に $m=5,n=4m=5,n=4$ で１万回に１回程度、$m=5,n=5m=5,n=5$ の場合で１０万回に１回程度発生します。$m=5,n=6m=5,n=6$ の場合では全く発生しません。観測値がすべて近しい値の場合は最尤推定量 $m^2\hat{m}_{2}$ の値が大きくなり、発散しやすい傾向にあります。

比較方法

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

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

結果

• $mm$ の推定：確率紙による推定方法１よい。観測値の数が約５０以上の場合は、最尤推定による推定方法２がよい。
• $η\eta$ の推定：最尤推定による推定方法２がよい。欠点として、両者ともに $mm$ が小さい場合は、$η\eta$ の推定値は正の無限大に発散するため、観測値の数を大きくする必要あり。

Python ソースコード

こちらをご覧ください。

[1] 逆関数法によるワイブル分布に従う乱数の生成

確率密度関数 $f(x)f(x)$、累積分布関数 $F(x)F(x)$ で表される分布に従う乱数の生成方法は、標準一様分布 $U(0,1)U(0,1)$ から得られる乱数 $uu$ を累積分布関数の逆変換 $(F−1)( F^{-1} )$ をすることにより可能となります。

ワイブル分布の累積分布関数 $F(x)F(x)$ の逆関数を求めると、

$$F(x)=1−exp⁡{−(xη)m}1−F(x)=exp⁡{−(xη)m}ln⁡(1−F(x))=−(xη)mxm=−ηmln⁡(1−F(x))x=η−ln⁡(1−F(x))m\begin{aligned} F(x)&=1-\exp\left\{-\left(\dfrac{x}{\eta} \right)^{m} \right\} \\ 1-F(x)&=\exp\left\{-\left(\dfrac{x}{\eta} \right)^{m} \right\} \\ \ln{(1-F(x))} &= -\left( \dfrac{x}{\eta} \right)^{m}\\ x^{m} &= -\eta^{m}\ln{(1-F(x))}\\ x &=\eta\sqrt[\footnotesize m]{-\ln{(1-F(x))}} \end{aligned}$$

よって、累積分布関数 $F(x)F(x)$ の逆関数 $F−1(x)F^{-1}(x)$ は

$$F−1(x)=η−ln⁡(1−x)m\begin{aligned} F^{-1}(x)&=\eta\sqrt[\footnotesize m]{-\ln{(1-x)}} \end{aligned}$$

となります。
以上より、$U(0,1)U(0,1)$ に従う一様乱数 $uu$ を生成後、 $x=F−1(u)x=F^{-1}(u)$ と変換すれば、ワイブル分布 $W(m,η)W(m,\eta)$ に従う乱数を生成することができます。

[2] 方法１：確率紙によるパラメータ推定

まず、$ln⁡ln⁡11−F(x)\ln{\ln{\dfrac{1}{1-F(x)}}}$ を計算すると、

$$ln⁡ln⁡11−F(x)=ln⁡[−ln⁡{1−F(x)}]=ln⁡(xη)m=m(ln⁡x−ln⁡η)\begin{aligned} \ln{\ln{\dfrac{1}{1-F(x)}}} &= \ln{[-\ln{\{1-F(x)\}}]} \\ &= \ln{\left(\dfrac{x}{\eta} \right)^m} \\ &= m( \ln{x} - \ln{\eta} ) \end{aligned}$$

と表すことができます。 $y=ln⁡[−ln⁡{1−F(x)}]y=\ln{[-\ln{\{1-F(x)\}}]}$ とすると $yy$ と $ln⁡x\ln{x}$ は線形の関係にあることが分かります。
よって、$X∼W(m,η)X \sim W(m,\eta)$ から互いに独立に得られた $nn$ 個の観測値を小さい順に並べたものを $x(1), ... ,x(n)x_{(1)},\ ...\ ,x_{(n)}$ とし、 $yi=ln⁡[−ln⁡{1−F(x(i))}]y_{i} = \ln{\left[- \ln{\{ 1-F(x_{(i)})\}} \right]}$ とした場合、$ln⁡x(i)\ln{x_{(i)}}$ と $yiy_{i}$ をプロットすると直線になります。累積分布関数 $F(x(i))F(x_{(i)})$ の推定については Benard の近似式を使用します。

$$F(x(i))≈i−0.3n+0.4\begin{aligned} F(x_{(i)}) \approx \dfrac{i-0.3}{n+0.4} \end{aligned}$$

累積分布関数の推定に関してはこちらの記事をご覧ください。

$ln⁡x(i)\ln{x_{(i)}}$ と $yiy_{i}$ について $y=aln⁡x+by=a \ln{x} + b$ と直線近似することを考えます。
最小二乗法により、係数 $a,ba,b$ については

$${a=∑i=1n(ln⁡x(i)−ln⁡x‾)(yi−yˉ)∑i=1n(ln⁡x(i)−ln⁡x‾)2b=yˉ−a ln⁡x‾\left\{\begin{aligned} a&=\dfrac{\sum_{i=1}^{n}(\ln{x_{(i)}} - \overline{\ln{x}})(y_{i}-\bar{y})} {\sum_{i=1}^{n}(\ln{x_{(i)}}- \overline{\ln{x}})^{2}} \\ b&= \bar{y} - a\ \overline{\ln{x}} \end{aligned}\right.$$

$$(yˉ=1n∑i=1nyi , ln⁡x‾=1n∑i=1nln⁡x(i))\left( \begin{aligned} \bar{y}= \dfrac{1}{n}\sum_{i=1}^{n} y_{i}\ , \ \overline{\ln{x}} =\dfrac{1}{n} \sum_{i=1}^{n} \ln{x_{(i)}} \end{aligned} \right)$$

で表されます。 $a=m,b=−mln⁡ηa=m, b=-m\ln{\eta}$ から、確率紙による $m,ηm,\eta$ の推定量 $m^1,η^1\hat{m}_{1},\hat{\eta}_{1}$ は

$${m^1=∑i=1n(ln⁡x(i)−ln⁡x‾)(yi−yˉ)∑i=1n(ln⁡x(i)−ln⁡x‾)2η^1=exp⁡[−yˉ−m^1ln⁡x‾m^1]\left\{\begin{aligned} \hat{m}_1 &= \dfrac{\sum_{i=1}^{n}(\ln{x_{(i)}} - \overline{\ln{x}})(y_{i}-\bar{y})} {\sum_{i=1}^{n}(\ln{x_{(i)}}- \overline{\ln{x}})^{2}} \\ \hat{\eta}_1 &= \exp{\left[ -\dfrac{\bar{y}-\hat{m}_1\overline{\ln{x}}}{\hat{m}_1} \right]} \end{aligned}\right.$$

となります。

[3] 方法２：ワイブル分布の最尤推定によるパラメータ推定

ワイブル分布の尤度関数を $L(m,η)L(m,\eta)$ とすると、観測値を $x1, ... ,xnx_{1},\ ...\ ,x_{n}$ として、

$$L(m,η)=∏i=1nmη(xiη)m−1exp⁡[−(xiη)m]=(mη)n{∏i=1n(xiη)m−1}exp⁡[−∑i=1n(xiη)m]\begin{aligned} L(m,\eta) &=\prod_{i=1}^{n}\dfrac{m}{\eta}\left( \dfrac{x_{i}}{\eta} \right) ^{m-1} \exp{\left[- \left( \dfrac{x_{i}}{\eta} \right)^{m} \right]}\\ &=\left( \dfrac{m}{\eta} \right)^{n} \left\{ \prod_{i=1}^{n}\left( \dfrac{x_{i}}{\eta} \right) ^{m-1} \right\} \exp{\left[ -\sum_{i=1}^{n}\left( \dfrac{x_{i}}{\eta}\right)^{m} \right]} \end{aligned}$$

対数尤度関数については

$$ln⁡L(m,η)=n(ln⁡m−ln⁡η)+(m−1)∑i=1n(ln⁡xi−ln⁡η)−∑i=1n(xiη)m\begin{aligned} \ln{L(m,\eta)} &=n(\ln{m} -\ln{\eta}) + (m-1)\sum_{i=1}^{n}(\ln{x_{i}}-\ln{\eta}) -\sum_{i=1}^{n}\left( \dfrac{x_{i}}{\eta}\right)^{m} \\ \end{aligned}$$

対数尤度関数が最大となる $m,ηm,\eta$ を求めると、以下の２式が成り立ちます。

$${∂∂mln⁡L(m,η)=0∂∂ηln⁡L(m,η)=0\left\{\begin{aligned} \dfrac{\partial}{\partial m}\ln{L(m,\eta)} = 0 \\ \dfrac{\partial}{\partial \eta}\ln{L(m,\eta)} = 0 \\ \end{aligned}\right.$$

$η\eta$ については

$$∂∂ηln⁡L(m,η)=−nη+(m−1)∑i=1n(−1η)−∑i=1n(−mximηm+1)=−mnη+mη∑i=1n(xiη)m=mη(−n+∑i=1n(xiη)m)=0\begin{aligned} \dfrac{\partial}{\partial \eta}\ln{L(m,\eta)} &=-\dfrac{n}{\eta} + (m-1) \sum_{i=1}^{n}\left( -\dfrac{1}{\eta} \right) -\sum_{i=1}^{n}\left( -m\dfrac{x_{i}^{m}}{\eta^{m+1}}\right) \\ &=-\dfrac{mn}{\eta} +\dfrac{m}{\eta}\sum_{i=1}^{n}\left( \dfrac{x_{i}}{\eta}\right)^{m} \\ &=\dfrac{m}{\eta} \left(-n + \sum_{i=1}^{n}\left( \dfrac{x_{i}}{\eta}\right)^{m} \right) =0\\ \end{aligned}$$

よって、

$$n=∑i=1n(xiη)mηm=1n∑i=1nximη=1n∑i=1nximm\begin{aligned} n &=\sum_{i=1}^{n}\left( \dfrac{x_{i}}{\eta}\right)^{m}\\ \eta^{m} &=\dfrac{1}{n} \sum_{i=1}^{n} x_{i}^{m} \\ \eta &= \sqrt[\footnotesize m]{\dfrac{1}{n} \sum_{i=1}^{n} x_{i}^{m}} \end{aligned}$$

$mm$ については

$$∂∂mln⁡L(m,η)=nm+∑i=1n(ln⁡xi−ln⁡η)−∑i=1n(xiη)m(ln⁡xi−ln⁡η)=nm+∑i=1nln⁡xi−nln⁡η+ln⁡η∑i=1n(xiη)m−∑i=1n(xiη)mln⁡xi=nm+∑i=1nln⁡xi−1ηm∑i=1nximln⁡xi(∵n=∑i=1n(xiη)m)=nm+∑i=1nln⁡xi−n∑i=1nximln⁡xi∑i=1nxim=n(1m+1n∑i=1nln⁡xi−∑i=1nximln⁡xi∑i=1nxim)=0\begin{aligned} \dfrac{\partial}{\partial m}\ln{L(m,\eta)} &=\dfrac{n}{m} + \sum_{i=1}^{n}(\ln{x_{i}}-\ln{\eta}) -\sum_{i=1}^{n} \left( \dfrac{x_{i}}{\eta}\right)^{m} (\ln{x_{i}} -\ln{\eta})\\ &=\dfrac{n}{m} + \sum_{i=1}^{n}\ln{x_{i}}-n\ln{\eta} +\ln{\eta} \sum_{i=1}^{n}\left( \dfrac{x_{i}}{\eta}\right)^{m} -\sum_{i=1}^{n} \left( \dfrac{x_{i}}{\eta}\right)^{m} \ln{x_{i}} \\ &=\dfrac{n}{m} + \sum_{i=1}^{n}\ln{x_{i}} -\dfrac{1}{\eta^{m}}\sum_{i=1}^{n} x_{i}^{m} \ln{x_{i}} \left(\because n = \sum_{i=1}^{n}\left( \dfrac{x_{i}}{\eta}\right)^{m} \right)\\ &=\dfrac{n}{m} + \sum_{i=1}^{n}\ln{x_{i}} -n\dfrac{\sum_{i=1}^{n} x_{i}^{m} \ln{x_{i}}}{\sum_{i=1}^{n} x_{i}^{m}} \\ &=n\left( \dfrac{1}{m} + \dfrac{1}{n}\sum_{i=1}^{n}\ln{x_{i}} -\dfrac{\sum_{i=1}^{n} x_{i}^{m} \ln{x_{i}}}{\sum_{i=1}^{n} x_{i}^{m}} \right) =0 \\ \end{aligned}$$

となります。 $mm$ については解くことができないため、ニュートン法を使用します。

$$g(m)=1m+1n∑i=1nln⁡xi−∑i=1nximln⁡xi∑i=1nxim\begin{aligned} g(m)= \dfrac{1}{m} + \dfrac{1}{n}\sum_{i=1}^{n}\ln{x_{i}} -\dfrac{\sum_{i=1}^{n} x_{i}^{m} \ln{x_{i}}}{\sum_{i=1}^{n} x_{i}^{m}} \end{aligned}$$

とすると

$$g′(m)=−1m2−∑i=1nxim(ln⁡xi)2∑i=1nxim−∑i=1nximln⁡xi(−∑i=1nximln⁡xi(∑i=1nxim)2)=−1m2−∑i=1nxim(ln⁡xi)2∑i=1nxim+(∑i=1nximln⁡xi∑i=1nxim)2\begin{aligned} g'(m) &= -\dfrac{1}{m^{2}} - \dfrac{\sum_{i=1}^{n} x_{i}^{m} (\ln{x_{i}})^{2}}{\sum_{i=1}^{n} x_{i}^{m}} -\sum_{i=1}^{n} x_{i}^{m} \ln{x_{i}} \left( - \dfrac{\sum_{i=1}^{n} x_{i}^{m} \ln{x_{i}}} { \left( \sum_{i=1}^{n} x_{i}^{m} \right)^{2} } \right) \\ &= -\dfrac{1}{m^{2}} - \dfrac{\sum_{i=1}^{n} x_{i}^{m} (\ln{x_{i}})^{2}}{\sum_{i=1}^{n} x_{i}^{m}} + \left( \dfrac{\sum_{i=1}^{n} x_{i}^{m} \ln{x_{i}}} { \sum_{i=1}^{n} x_{i}^{m} } \right)^{2} \\ \end{aligned}$$

となります。
初期値 $m0m_{0}$ を方法１で推定した $m^1\hat{m}_{1}$ として、以下の漸化式を反復計算することで、方法２の $mm$ の最尤推定量 $m^2\hat{m}_{2}$ を求めます。

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

$η\eta$ の最尤推定量 $η^2\hat{\eta}_{2}$ は $m^2\hat{m}_{2}$ を使用して、

$$η^2=1n∑i=1nxim^2m^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}$$

より求めます。

簡易シミュレーションによる比較結果

前置き

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

• 方法１：確率紙による推定
• 方法２：最尤推定量による推定

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

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

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

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

$$θ^mean=1N∑i=1Nθ^i, Sθ^2=1N−1∑i=1N(θ^i−θ^mean)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}$$

と表せます。

結果

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

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

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

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

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

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

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

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

まとめ

以上の結果から、$mm$ の推定については確率紙による推定方法１が勝りますが、乱数の生成数 = 観測値の数が約５０以上の場合は、最尤推定による推定方法２が優れます。乱数の生成数 $nn$ が小さい場合は相対誤差が大きいことから、両方法の推定量は $mm$ の不偏推定量ではないことが分かります。$nn$ が大きい場合は、相対誤差、分散ともに０に近づくので、一致推定量であると判断できます。

$η\eta$ による推定は乱数の生成数に関係なく、最尤推定による推定方法２が優れていますが、両者ともに $mm$ が小さい場合は、$η\eta$ の推定値は正の無限大に発散するため、推定値の分散が大きくなり、観測値の数＝サンプルサイズは大きくする必要があります。

Python ソースコード

こちらをご覧ください。