故障確率の推定と寿命予測

はじめに

特定の操作Oを行ったときに、確率pでサージ電圧が掛かるモジュールがある。当該のモジュールはn-1回までサージ電圧に耐えられるが、n回目のサージ電圧で故障してしまう。

当該のモジュールは繰り返し操作Oを一日あたりm回行う条件下で動作する。一度故障したモジュールは交換され、同じ条件下で動作を続ける。

というモジュールとシステムについて、故障予測を行います。

寿命予測がらみの実験結果の解析手順としてそこそこ汎用なのではないかなと思いますが、あまりググっても知見が引っかからなかったので我流でまとめました。調べ方が悪いのかもしれません。

信頼性工学は専攻したわけではなく素人です。指摘がほしいです。

以下、演算に用いてるコードはPython3.7によるものです。ライブラリのバージョンは適当なので、環境によっては合わないかもしれません。

判っているデータ

判っているデータは「故障が発生するまでに操作Oを何回試行する必要があったか」というものです。

今回はダミーデータとして下記の値を使います。この値はモジュールを20回故障させ、それぞれ何回目の操作Oで故障したかを記録したものを想定しています。（乱数シミュレーションで生成しています）

上記の値は data.csv として保存しておきます。

生成過程のパラメータ推定

何回目の操作Oで故障するかの分布は $\mu = \frac{1}{p}, n = n$ のガンマ分布に従うと仮定できます。

この仮定をもとに、上記データからpとnを推定します。

from scipy.stats import gamma
import pandas as pd

# データ読み込み
df = pd.read_csv("data.csv", index_col=None, header=None)
data = df.values.reshape([-1])

# 最尤推定によるパラメータ推定。
# 今回はlocだけは固定しておく。（分布の性質から、ずらす必要がないため）
a_hat, loc_hat, scale_hat = gamma.fit(data, floc=0)

# scipyのgammaのaはn、scaleはmuのこと。
# muからpを計算しておく。
param_n = a_hat
param_mu = scale_hat
param_p = 1 / param_mu

print(f"n: {param_n} p: {param_p}")

結果は以下になりました。

n: 2.2778957910474036 p: 0.002004748770998815

先程のデータの生成はn=3、p=0.003で行っています。推定されたnとpは少しずれていますが、これはサンプル数20で推定したため、大きめの誤差が生じたためです。*1

サンプル数10000で推定すると n=3.04、p=0.00304 とそれぞれ1%程度の誤差に落ちました。

ここまでをまとめると以下の結論となります。

このモジュールは操作Oを行ったときに約0.2%の確率で内部にサージ電圧が発生し、更に平均約2.3回のサージ電圧が発生すると故障に至る。

寿命予測

ここからは上記で推定した生成過程において、では何回の操作Oに耐えられるのかという問題に対する答えを予測します。

まず最初に、求めたガンマ分布の累積分布関数の概形を確認します。

f:id:marineotter:20200405122923p:plain — ガンマ累積分布関数

このグラフは「行う操作Oの回数x」を横軸に、「x回の操作を行った際に、x回目までで一度でも故障した確率y」を縦軸に取ったものとなっています。

ここで、この「一度でも故障した確率」が0.1を上回った時点までの操作Oの回数を「寿命が来るまでの回数」と定義することとします。*2

グラフより、操作Oを340回試行すると「一度でも故障した確率」は0.1*3より上になります。よって、340回の試行で故障するという寿命予測となります。

ここで最初に戻って、1日あたりの操作Oの試行回数mを1と仮定した場合、寿命は340日と予測できます。*4

ここで使用した計算スクリプトは以下のとおりです。

import matplotlib.pyplot as plt
import numpy as np

n = 2.774251340060173
p = 0.0028411606739312542

x = np.arange(0, 3000, 0.1)
y = gamma.cdf(x, a=n, loc=0, scale=1/p)

plt.plot(x, y)
plt.savefig("cdf.png")
plt.clf()

p_thr = 0.1
predicted_lifespan = x[np.argmin(np.abs(y - p_thr))]
print(f"lifespn_predicted: {int(np.round(predicted_lifespan))}")

おわりに

以上のような感じでふんわり寿命予測してみましたが、どこかで何かを大きく間違えている気もします。

手順を進めればすすめるほどに自信がなくなっているので、指摘とかいただけると非常に助かります。

ちなみにおまけですが、今回の計算の原理は以下の予測と同一です。

HPがXの敵キャラに対して、期待ダメージX/n、命中率pの攻撃を行って、倒すまでに必要な攻撃回数の予測。

*1:十分なn数が得られない場合、このフェーズで推定されるパラメータが大きな誤差を含んでいる点が妥当性への脅威となりえます。

*2:この定義が信頼性工学の観点から正しいのか判っていません。

*3:「0.1」という値は「宣言している寿命より短い寿命で故障してしまう確率」であり、勝手に決めました。

*4:この値はデータの平均からは大きく離れており、単純な平均値を予測として採用できない理由の一つである…でいいのですかね

川獺の外部記憶

なんでも残しておく闇鍋みたいな備忘録