ワイブル解析のやり方｜寿命データからm・ηを求める

試作したベアリングの寿命試験が終わりました。10個すべてが壊れるまで回し続け、平均寿命は929時間。報告書にそう書いたところ、設計担当から「10個に1個が壊れるのは何時間ですか」と聞かれました。

平均では答えられません。保証期間や交換周期を決めるときに知りたいのは平均ではなく、累積で何%が壊れる時点はいつかという値だからです。

寿命データはばらつきが左右非対称になりやすく、正規分布ではうまく表せません。ここで使うのがワイブル分布です。この記事では10個の実測値から2つのパラメータを推定し、B10寿命まで出す手順をExcelで解説します。

ワイブル解析を使う場面

ワイブル解析が向いているのは次のような場面です。

• 寿命や故障までの時間を扱うとき。回転数、繰り返し数、走行距離のように「使うほど蓄積する量」でも同じ扱いができます
• データ数が少ないとき。10個前後でも形状の見当がつきます。寿命試験はサンプルを壊すため数を確保しにくく、この特性が効きます
• 故障のタイプを知りたいとき。初期故障なのか摩耗故障なのかを、パラメータの値から判断できます

逆に、故障の原因が途中で変わっている場合は1本の直線に乗りません。初期不良と摩耗が混ざったデータをまとめて解析すると、どちらでもない中途半端なパラメータが出ます。プロットが折れ線に見えたら、原因別に分けてから解析してください。原因の切り分けにはFTAやFMEAが使えます。

例題｜ベアリング10個の寿命データ

試作ベアリング10個を同じ荷重・同じ回転数で連続運転し、振動値が基準を超えた時点を故障としました。故障時間（時間）を小さい順に並べます。

順位 i	故障時間 t（h）	順位 i	故障時間 t（h）
1	320	6	940
2	480	7	1,060
3	610	8	1,210
4	720	9	1,400
5	830	10	1,720

単純な平均は929時間です。ただし最短は320時間、最長は1,720時間と5倍以上の開きがあります。この非対称なばらつきを2つのパラメータで表すのがワイブル分布です。

ワイブル分布の累積分布関数（時刻 \( t \) までに故障する確率）は次の式です。

\[ F(t) = 1 – \exp\left\{ -\left(\frac{t}{\eta}\right)^{m} \right\} \]

\( m \) は形状パラメータで分布の形を、\( \eta \) は尺度パラメータで時間の大きさを決めます。この2つを実測データから求めるのがワイブル解析です。

メジアンランクでプロット用の座標を作る

先ほどの式を2回対数を取って変形すると、直線の形が現れます。

\[ \ln\left\{ \ln\frac{1}{1-F(t)} \right\} = m \ln t – m \ln \eta \]

左辺を \( y \)、\( \ln t \) を \( x \) と置けば \( y = mx – m\ln\eta \) です。傾きが m そのものです。あとは各データ点の \( F(t) \) を決めれば、散布図に点を打てます。

この \( F(t) \) の推定に使うのがメジアンランク法です。i番目に壊れたサンプルの累積故障確率を次の近似式で与えます。

\[ F_i = \frac{i – 0.3}{n + 0.4} \]

単純に \( i/n \) としないのは、10個目のデータで \( F=1 \) になり対数が計算できなくなるためです。ExcelではA列に順位、B列に故障時間を入れて次のように組みます。

C2: =(A2-0.3)/(10+0.4)      メジアンランク
D2: =LN(B2)                  x
E2: =LN(-LN(1-C2))           y

10行分を計算した結果です。

i	t（h）	F	x = ln t	y
1	320	0.0673	5.7683	−2.6638
2	480	0.1635	6.1738	−1.7233
3	610	0.2596	6.4135	−1.2020
4	720	0.3558	6.5793	−0.8217
5	830	0.4519	6.7214	−0.5086
6	940	0.5481	6.8459	−0.2304
7	1,060	0.6442	6.9660	0.0329
8	1,210	0.7404	7.0984	0.2990
9	1,400	0.8365	7.2442	0.5940
10	1,720	0.9327	7.4501	0.9927

途中で試験を打ち切った場合

全部が壊れる前に試験を止めた場合や、故障以外の理由で外れたサンプルがある場合は、この単純なメジアンランクは使えません。打ち切りを考慮した平均順位法で順位を補正する必要があります。まずは全数が故障するまで試験した完全データで手順を身につけてください。

傾きと切片から m・η を求める

x と y の散布図を描き、最小二乗法で直線を当てはめます。傾きと切片はSLOPE関数とINTERCEPT関数で直接出ます。

=SLOPE(E2:E11, D2:D11)      → 2.1835   ← m
=INTERCEPT(E2:E11, D2:D11)  → -15.2097 ← -m·lnη
=RSQ(E2:E11, D2:D11)        → 0.9991   ← 直線への当てはまり

形状パラメータは \( m = 2.18 \) と求まりました。尺度パラメータは切片から逆算します。切片が \( -m\ln\eta \) なので、次の式です。

\[ \eta = \exp\left(-\frac{b}{m}\right) = \exp\left(\frac{15.2097}{2.1835}\right) = 1060 \]

=EXP(-INTERCEPT(E2:E11,D2:D11)/SLOPE(E2:E11,D2:D11))  → 1059.6

尺度パラメータは約1,060時間です。\( \eta \) は「累積63.2%が故障する時間」という意味を持ちます。式に \( t = \eta \) を入れると \( F = 1 – e^{-1} = 0.632 \) になるためです。平均寿命ではない点に注意してください。

決定係数は0.9991で、10点がほぼ直線に乗っています。プロットが曲がっている場合はワイブル分布が合っていないか、故障の原因が混ざっています。決定係数の読み方も参考にしてください。

m が語る故障のタイプ

形状パラメータmは、故障率が時間とともにどう動くかを表します。バスタブ曲線のどの領域にいるかが1つの数値で分かる、というのがワイブル解析の実務的な価値です。

mの値	故障率	該当する領域
m < 1	時間とともに減る	初期故障期
m = 1	一定	偶発故障期（指数分布）
m > 1	時間とともに増える	摩耗故障期

今回は \( m = 2.18 \) なので摩耗故障です。使うほど壊れやすくなる、つまり寿命で決まっている状態を示します。

この読み分けは対策の方向を変えます。摩耗故障なら定期交換や予防保全が効きます。交換周期を決めれば故障を先回りできるからです。一方でmが1を下回る初期故障なら、定期交換しても意味がありません。新品ほど壊れやすいので、交換するたびにリスクが戻ってしまいます。この場合は初期不良を選別するエージングや、製造工程の見直しが対策です。

まとめると、m > 1 なら保全計画、m < 1 なら工程改善。数値を出す目的はここにあります。

B10寿命と信頼度を計算する

冒頭の「10個に1個が壊れるのは何時間か」に答えます。累積故障確率が10%になる時間をB10寿命と呼び、ベアリングや電子部品の寿命表記に広く使われます。

累積分布関数を \( t \) について解くと次の式です。

\[ t = \eta \left\{ -\ln(1 – F) \right\}^{1/m} \]

F = 0.1 を代入します。

\[ B_{10} = 1059.6 \times \left\{ -\ln(0.9) \right\}^{1/2.1835} = 378 \]

=1059.6*(-LN(0.9))^(1/2.1835)  → 378.1

B10寿命は378時間。単純平均929時間の4割程度です。保証や交換周期を平均で決めると、10%の個体がその前に壊れます。この差が、寿命データを平均で報告してはいけない理由です。

逆に、時間を指定して生存率を求めることもできます。信頼度 \( R(t) = 1 – F(t) \) です。

経過時間	信頼度 R(t)	意味
378h（B10）	90.0%	10個に1個が故障
500h	82.4%	6個に1個が故障
1,000h	41.4%	半数以上が故障

=EXP(-(500/1059.6)^2.1835)              → 0.8237
=1-WEIBULL.DIST(500, 2.1835, 1059.6, TRUE)  → 0.8237

WEIBULL.DIST関数を使えば累積確率を直接出せます。第4引数のTRUEが累積分布関数です。

参考までに、ワイブル分布の平均寿命はガンマ関数を使って \( \eta \Gamma(1 + 1/m) \) で求まり、938時間です。10個の単純平均929時間とは一致しません。分布から求めた期待値と、手元の10個の標本平均は別物として扱ってください。

よくある質問

Q. データは何個あればワイブル解析できますか？

A. 実務では7個以上が目安です。5個程度でも計算はできますが、mの推定誤差が大きく、故障タイプの判定が揺れます。特にmの値が1をまたぐ付近では判断が変わってしまうため、少数のデータでm=1.2だから摩耗故障、と言い切るのは避けてください。

Q. ηと平均寿命はどう違いますか？

A. ηは累積63.2%が故障する時間、平均寿命は分布の期待値です。今回の例ではη=1,060時間、平均寿命=938時間で一致しません。m=1（指数分布）のときだけ両者が一致します。仕様書に書くときはどちらを指すか明記してください。

Q. プロットが直線にならないときはどうしますか？

A. 折れ線に見えるなら故障モードが2つ混ざっています。故障した部位や様子で分類し、モードごとに解析し直してください。全体が緩やかに曲がる場合は、故障が起きない期間（位置パラメータ）がある3母数ワイブル分布を検討します。

高温試験の結果を実使用の年数へ換算するには、別に加速係数が要ります。アレニウス則による換算は加速試験とは｜アレニウス則で寿命を見積もるで解説しています。

まとめ

• 寿命データは左右非対称になりやすく、平均だけでは保証期間を決められない
• メジアンランク \( F_i = (i-0.3)/(n+0.4) \) で座標を作り、\( \ln t \) と \( \ln\{\ln(1/(1-F))\} \) の散布図を描く
• 傾きがm、切片から \( \eta = \exp(-b/m) \)。ExcelのSLOPE・INTERCEPT・EXPで求まる
• m > 1 は摩耗故障で予防保全が効く。m < 1 は初期故障で定期交換は効かない
• 例題では m=2.18、η=1,060h、B10寿命は378hで単純平均929hの4割

まとめると、交換周期を決めたいならB10寿命、対策の方向を決めたいならm、という使い分けです。報告書には平均ではなくこの2つを載せてください。

この記事の前提となるバスタブ曲線や信頼度の考え方は信頼性工学の基礎で、確率分布の選び方全体は確率分布の種類と使い分けで解説しています。あわせてご確認ください。