Matlabによる機器の寿命解析

■背景
Rでタイトルのようなことをやっていたのですが、ちょこちょことアルゴリズムの部分で分からないことが増えてきたので、先ずはアルゴリズム理解として慣れたツールで練習したいと思います。

ちなみにRで躓いたのは、ワイブル分布の推定パラメータを2つでやっていたところ実は3パラメータで考える必要が出てきて、それにまつわる周辺の諸計算を上手く出来なかった点です。

Matlabなら何とかなるかもしれません…。

先ずはMatlabにStatistics and machine learning toolboxをインストールします。これで統計処理の関数がいろいろと使えるようになります。

そして問題はというと…

f:id:cocosuzu:20170602210027j:plain

この図でtが小さいほどワイブル分布の近似曲線から離れている点が問題です。これは2パラメータのワイブル分布では推定しきれない特性なので、第3のパラメータをワイブル分布関数で考慮する必要があります。したがって、今回は以下を課題とします。

【課題: 3パラメーターワイブル分布のパラメーター推定】

2パラメータワイブル分布！？3パラメータ！？となるので先にワイブル分布について書いておきます。

■ワイブル分布の基本

統計分布を使って事象を分析する際に一番重要なのは事象をどの統計で表すと良いのかを判定することです。ただし、ここではワイブル分布になると分かっているのでその点はスルーします。その他の分布には指数分布、対数分布などがあります。

ワイブル分布が生まれた経緯を確認すると分かりやすいので先ずはそのあたりから触れていきます。

n個の環によって構成された鎖において、一つの環がストレス $x$ で故障する累積確率分布関数を $G(x)$ とすると、故障しない確率は鎖全体で以下のようになります。

$\displaystyle \left(1-G(x)\right)^{n} \tag{1} \label{1}$

一方で鎖に $x$ のストレスを加えた際に、鎖の一部が壊れる累積確率分布関数を $F(x)$ とすると、壊れない確率は以下で表現出来ます。

$\displaystyle 1-F(x) \tag{2} \label{2}$

式\eqref{2}と式\eqref{3}は同じことをさしているので以下の等式が成り立ちます。

$\displaystyle 1-F(x)=\left(1-G(x)\right)^{n} \tag{3} \label{3}$

ここで $G(x)$ を $(1-\exp(-\phi(x))$ とすると

$\displaystyle F(x)=1-\exp(-n\phi(x)) \tag{4} \label{4}$

ここでワイブル分布の特徴的な部分として $\phi(x)$ を以下のように仮定します。

$\displaystyle n\phi(x)=\left(\frac{x-\gamma }{\eta }\right)^m \tag{5} \label{5}$

この仮定によって、故障の特性を広範囲で上手く推定できるそうです。この仮定の下、式\eqref{5}は以下のようになります。

$\displaystyle F(x)=1-\exp(-\left(\frac{x-\gamma }{\eta }\right)^m) \tag{6} \label{6}$

一般的に\eqref{6}のパラメータはそれぞれ $\gamma$ :位置パラメータ、 $m$ :形状パラメータ、 $\eta$ :尺度パラメータと呼ばれます。鎖の場合は1変数をストレス $x$ としましたが、現象によっては時間 $t$ とする場合もあります。

以上より実験結果などから3パラメータを推定することでワイブル分布で故障特性を表すことが出来る下準備が出来ました。

■ワイブル解析によるMTTFの推定

寿命の指標としてMTTFという数値があります。今回はこれを求めるために必要なパラメータをワイブル解析により求めましょうというのが趣旨です。

ワイブル分布に従う確率密度分布関数 $f(t)$ は式\eqref{6}の微分結果として以下の式で表されます。

$\displaystyle f(t|m,\eta)=\frac{m}{\eta}\left(\frac{t-\gamma}{\eta}\right)^{m-1}\exp\left(-\left(\frac{t-\gamma}{\eta}\right)^{m}\right) \tag{7} \label{7}$

また、MTTFの定義式は以下の通りです。

$\displaystyle {\rm MTTF} = \int_{0}^{\infty}tf(t)dt \tag{8} \label{8}$

確率密度関数を代入してみると以下のような形になります。

$\displaystyle {\rm MTTF} = \int_{0}^{\infty}t\frac{m}{\eta}\left(\frac{t-\gamma}{\eta}\right)^{m-1}\exp\left(-\left(\frac{t-\gamma}{\eta}\right)^{m}\right)dt \tag{9} \label{9}$

むむむ…ネイピア数の指数関数と変数の掛け算…どこかで見たことのある形…

そう！ガンマ関数！ということでガンマ関数の定義式を見てみましょう。

$\displaystyle \Gamma(x)=\int_{0}^{\infty}t^{x-1}\exp(-x)dt \tag{10} \label{10}$

無理やりガンマ関数の形に合わせるには、 $\exp$ の中身をすっきりさせたいので

$\displaystyle \left(\frac{t-\gamma}{\eta}\right)^m=u \tag{11} \label{11}$

と置いてみましょう。これによって $\exp$ 部分は以下のようになります。

$\displaystyle \exp\left(-\left(\frac{t-\gamma}{\eta}\right)^{m}\right)=\exp\left(-u\right) \tag{12} \label{12}$

両辺を $t$ で微分すると

$\displaystyle \frac{d}{dt}\exp\left(-\left(\frac{t-\gamma}{\eta}\right)^{m}\right)=\frac{du}{dt}\frac{d}{du}\exp\left(-u\right) \tag{13} \label{13}$

整理すると

$\displaystyle \frac{m}{\eta}\left(\frac{t-\gamma}{\eta}\right)^{m-1}\exp\left(-\left(\frac{t-\gamma}{\eta}\right)^{m}\right)dt=\exp\left(-u\right)du \tag{14} \label{14}$

また、式\eqref{11}を $t$ について解くと

$\displaystyle t=\eta u^{\frac{1}{m}}+\gamma \tag{15} \label{15}$

となり、式\eqref{9}を以下のように変換することができます。

$\displaystyle {\rm MTTF} = \int_{0}^{\infty} \left(\eta u^{\frac{1}{m}}+\gamma \right)\exp\left(-u\right)du \tag{16} \label{16}$

ガンマ関数式\eqref{10}と若干違うので、さらに無理やり微調節すると

$\displaystyle {\rm MTTF} = \eta \int_{0}^{\infty} u^{\left(\frac{1+m}{m}-1\right)}\exp\left(-u\right)du + \gamma \int_{0}^{\infty}\exp\left(-u\right) \tag{17} \label{17}$

ようやくガンマ関数の形を作ることが出来ました。したがってMTTFは

$\displaystyle {\rm MTTF} = \eta \Gamma\left(\frac{1+m}{m}\right)+\gamma \tag{18} \label{18}$

と簡単な形で表すことが出来ます。

■ワイブル解析の実際

ちょうど良さそうな例がMathworksのサイト【記事1】にあったのでこれを参考に進めます。

【記事1】
jp.mathworks.com

この例だと今まで2パラメータのワイブル分布推定しかしていなかったところを、3パラメータまで推定することが出来ます。2パラメータのワイブル推定では式\eqref{6}において $\gamma=0$ と考えているため、残りの2パラメータである $m$ 、 $\eta$ のみを推定しています。 $\gamma=0$ とはどんな条件下かというと、動作が開始した $t=0$ の直後から少なからず故障発生の確率がある場合です。

実際には、部品によってある程度の時間が経過しないと絶対に故障しない場合があったりします。そんな場合に $\gamma=0$ ではワイブル分布のプロットとると、ある時間tで累積故障率がすとんと落ちて急速にゼロに向かいます。これを2パラメータのワイブル推定でフィッティングすると、冒頭のようにtの小さい領域でフィッティングラインと実測値がどんどん離れていくことになります。

対して3パラメータのワイブル分布では $\gamma$ を調整することによって小さなtの領域での乖離を減らすことができ、測定結果に則したワイブル分布を正しく推定することができます。

ただ、まとめるのに疲れてしまったので、とりあえずここまでをアップして続きは後日追記します。

参考記事
【記事2】http://www.yukms.com/biostat/takahasi2/rec/archive/takahashi_04_2(2015_04_10).pdf