Python でワイブル分布の計算をする方法を調べてみました。結果を記載します。

ワイブル分布とは

ワイブル分布 - Wikipedia の記載から内容を抜粋します。

  • 物体の強度を統計的に記述するためにW.ワイブル (Waloddi Weibull) によって提案された確率分布
  • 時間に対する劣化現象や寿命を統計的に記述するために利用する
  • 最弱リンクモデル

最弱リンクモデル という言葉が気になりましたが、ワイブル分布ってなに？～高校数学でわかるワイブル分布～ の説明がわかりやすかったです。

統計学入門 - 東京大学出版会 の6.13 <wbr>ワイブル分布 の気になった箇所を引用します。

故障が偶発故障なら、瞬間故障率は一定となり、これらの確率変数は指数分布に従うことが証明できる。
もし劣化が進行し故障率が増加するときはIFR(Increasing Failure Rate)といい、この分布には従わない。またいわゆる「初期故障」の時期には、故障率の減少がおこり、DFR(Decreasing Failure Rate) というが、この時も指数分布はあてはまらない。これらの現実により近いものがワイブル分布である。

故障率の計算には、ワイブル分布を使うのが良いと理解しました。

Python でのワイブル分布の計算方法

numpy.random.weibull を使う

numpy.random.weibull(a, size) の a は ワイブル係数（形状パラメータ）を表します。
a の値により、ワイブル分布の形状が変化し、a=1の場合 、指数分布になります。

import numpy as np
import matplotlib.pyplot as plt
# 
a = 1. # shape
s = np.random.weibull(a, 1000000)
plt.hist(s, bins=100)
plt.show()

201904080208_output_2_0.png - Google ドライブ

a=2 の場合、レイリー分布になります。

import numpy as np
import matplotlib.pyplot as plt
# 
a = 2. # shape
s = np.random.weibull(a, 1000000)
plt.hist(s, bins=100)
plt.show()

201904080208_output_4_0.png - Google ドライブ

形状パラメータ a=0.5、a=3.5 あたりはワイブル分布かと思いますが分布の形状は以下のようになります。

import numpy as np
import matplotlib.pyplot as plt
a = .5 # shape
s = np.random.weibull(a, 1000000)
plt.hist(s, bins=100)
plt.show()

201904080208_output_6_0.png - Google ドライブ

import numpy as np
import matplotlib.pyplot as plt
a = 3.5 # shape
s = np.random.weibull(a, 1000000)
plt.hist(s, bins=100)
plt.show()

201904080208_output_7_0.png - Google ドライブ

scipy.stats.weibull_min を使う

weibull_min.stats で モーメント を計算できます。

from scipy.stats import weibull_min
c = 1.79
mean, var, skew, kurt = weibull_min.stats(c, moments='mvsk')
print(mean, var, skew, kurt)


0.8895312083339716 0.26414885038021363 0.7869416254814524 0.5763152896528907

ppf pdf を使って、確率密度関数、パーセント点の計算が行えます。

import matplotlib.pyplot as plt
from scipy.stats import weibull_min
import numpy as np
c = 1.79
fig, ax = plt.subplots(1, 1)
# 1%点、99%点を100等分した配列を生成
x = np.linspace(weibull_min.ppf(0.01, c), weibull_min.ppf(0.99, c), 100)
# グラフを描く
ax.plot(x, weibull_min.pdf(x, c), 'r-', lw=5, alpha=0.6, label='weibull_min pdf')
rv = weibull_min(c)
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')
# ヒストグラムを描く
r = weibull_min.rvs(c, size=1000)
ax.hist(r, density=True, histtype='stepfilled', alpha=0.2)
ax.legend(loc='best', frameon=False)
plt.show()

201904080208_output_11_0.png - Google ドライブ

scipy.stats.weibull_max を使う

weibull_max.stats でモーメントの計算ができます。

from scipy.stats import weibull_max
c = 2.87
mean, var, skew, kurt = weibull_max.stats(c, moments='mvsk')
print(mean, var, skew, kurt)


-0.8862269254527579 0.21460183660255183 -0.6311106578189355 0.2450893006876287


from scipy.stats import weibull_max
import matplotlib.pyplot as plt
c = 1.79
fig, ax = plt.subplots(1, 1)
rv = weibull_max(c)
# 1%点、99%点を100等分した配列を生成
x = np.linspace(weibull_max.ppf(0.01, c), weibull_max.ppf(0.99, c), 100)
# グラフを描く
ax.plot(x, weibull_max.pdf(x, c), 'r-', lw=5, alpha=0.6, label='weibull_max pdf')
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')
r = weibull_max.rvs(c, size=1000)
ax.hist(r, density=True, histtype='stepfilled', alpha=0.2)
ax.legend(loc='best', frameon=False)
plt.show()

201904080208_output_14_0.png - Google ドライブ

scipy.stats.weibull_min と scipy.stats.weibull_max の違い

scipy.stats.weibull_min と、scipy.stats.weibull_max の違いをわかった限り記載します。

  • 計算式の違い

    • scipy.stats.weibull_min
      $f(x, c) = c x^{c-1} \exp(-x^c)$

    • scipy.stats.weibull_max
      $f(x, c) = c (-x)^{c-1} \exp(-(-x)^c)$

  • グラフの形状について
    scipy.stats.weibull_max は負の方向に、scipy.stats.weibull_min は正の方向に伸びます。

参考

以下、記事作成中に参考にしました。

コメント

カテゴリー

最近の投稿