R言語による機器の寿命解析
■背景
最近Pythonに集中しようと言いつつ、Javascriptやら他の言語にもに浮気しまくっています(笑)
先日、ある部品の寿命がどれだけか？ということを測定結果から推定することになったので色々と調べてみました。
そこで出会ったのがWeibull分布なる統計分布です。
その他の分布の代表的な例は【記事1】に分かりやすく取り上げられています。
【記事1】
qiita.com
今回取り扱う部品がWeibull分布で表される故障の特性になっているのでこれを使いましょうということです。また、Weibull分布を取り扱うに際して、特殊なグラフなどが必要になるのでソフトの選択肢は今のところ以下の通りです。
・R言語
・Matlab + statistics machine learnig tool box
・python + scipy
pythonでやりたかったのですが、参考例が上手く見つけられませんでした。また、matlabはアカウント数が限られているので他の人にも簡単に託せるフリーウェアのR言語でやっていきます。それに、Rやってるとか統計の玄人みたいでかっこいい…笑
■Rの開発環境
初めてRに触れるのですが、インストール自体は適当で大丈夫です!
Rの他にIDE?のようなものとしてRstudioをインストールします。これも適当で大丈夫です。
基本は***.Rファイルを作ってそれをconsoleでsource("***.Rのパス")と実行すれば動きます。
また、console上で1行づつ実行することも可能です。
xに1を代入
x <- 1
■Weibull分布と戯れる
まずWeibull分布をグラフに表すときには特殊な軸のグラフを使います。
【記事2】に特性も含めて概要が分かりやすく書いてありました。
【記事2】
www.atengineer.com
Weibullプロットの特殊な軸上にデータ点P(t,F(t))をプロットし、特性曲線からtに至るまでの累積故障率が分かるという訳です。これで、例えば全体の5%が壊れるtを製品の寿命と設定することが出来るようになります。
なんとなく分かったところで、実際にRで書くことが出来るのか試してみます。そのために参考にしたのが【記事3】です。
【記事3】
R -- ワイブル確率紙
無事描けました。
ちなみに実際の使用を想定して、data.txtファイルから故障までにかかった時間のデータ列を読み込んで表示するには、xへの代入を以下のようにするといいです。
x<-read.table("C:/Users/****/Documents/R/data.txt")列にヘッダーがある場合は引数で指定してあげると読み込んでくれます。
適当に数百時間～数千時間くらいで故障する10サンプル程度のデータを用意して描画してみると【画像1】みたいになります。
【画像1】
とりあえず描けるようになりました。ちょっとWeibullプロットとRに慣れてきました。
■Weibull分布とパラメータ推定
WeibullプロットはWeibull分布に則した故障の発生結果から求められており、その結果からtにおける故障確率密度f(t)を求めてあれこれするとMTTFが求まります。ここまではWeibullプロットという結果だけ見てきたのですが、パラメータの推定が目的なのでもう少しアプローチを変えます。
寿命であるMTTFを求めるにはWeibull分布のShape(形状)mとScale(尺度)というパラメータを求める必要があります。それっぽい参考例【記事4】があったので真似してみます。
【記事4】
stats.stackexchange.com
先ずはWeibull分布に則ったランダム数を作ります。このときサンプル数N=1000, m=1, としています。
データセットが出来たところでkernel密度推定を行っていますが、kernel密度推定ってなんですか！【記事5】に分かりやすく書いてありました。
【記事5】
http://www.agu.ac.jp/~nomura/lecture/archive/kd.pdf
なるほど！発生結果の見方を調整して確率密度分布の関数を見つけましょうねという感じでしょうか。
もう一つ大変分かりやすい解説が【記事6】にありました。
【記事6】
d.hatena.ne.jp
また、さっぱりとした説明として【記事7】を見つけました。
【記事7】
[R言語][Rstan]Rstanでワイブル分布のパラメータ推定 - gepulog
という感じでワイブル分布のパラメーターを推定することが出来ます。
その他参考記事
