Claude Codeでワイブル解析 — 打ち切りありの寿命データからB10寿命を出す

ワイブル解析で寿命データからB10寿命(10%の個体が故障する時間)を出すとき、結果を最も大きく動かすのは打ち切りデータの扱いです。試験終了時に壊れていなかった個体を「故障した」ことにすると、B10寿命が約19%長く出ます。この記事では、Claude Codeにlifelinesで当てはめを書かせ、手計算に近い別の方法で検算させるところまでを、実際に動かした数値つきで進めます。

打ち切りの渡し方で結果がどれだけ変わるか

lifelinesの WeibullFitter は、時間の配列と「故障を観測したか」の配列(1なら故障、0なら打ち切り)を受け取ります。lifelinesのドキュメントの例も wbf.fit(waltons['T'], waltons['E']) の形です。

次の20個の試験データで確かめます。試験は1000時間で打ち切り、14個が故障し、6個は1000時間時点で動いていました。

全部を故障扱いにすると、1000時間まで壊れなかった事実が「1000時間ちょうどで壊れた」に化けます。早く壊れる個体の割合は変わらないのに、形状が急になり、B10が長く出ます。捨てる場合は逆で、生き残りの情報を失うので短く出ます。どちらも「それらしい数字」が出るため、コードを読まないと気づけません。

Claude Codeに渡す指示とCLAUDE.mdの規約

Claudeが実行できる確認手段を渡すと、「できたように見える」で止まらなくなります。ワイブル解析での確認手段は、別実装による検算です。プロジェクトのCLAUDE.mdには次のように書いておきます。

## 寿命データ解析の規約
- 入力CSVの列は hours(時間)と failed(故障=1、打ち切り=0)に固定する(設備データのCSVを扱う流れは[Claude CodeでPLCのログ解析](/articles/claude-code-plc-sensor-log-csv-analysis)にあります)
- 打ち切りを捨てない。全件を WeibullFitter.fit(hours, failed) に渡す
- B10 は wf.percentile(0.90) で出し、閉形式 lambda * (-ln 0.9)^(1/rho) と一致を確認する
- scipy で対数尤度を直接最大化した lambda, rho と lifelines の値を突き合わせる
- 結果は rho の95%区間とB10のブートストラップ区間を併記する

Claudeへの依頼は、データの形と合格条件を一緒に書きます。次のような書き方です。

data/life_test.csv(hours, failed)を寿命試験の結果として読み、
weibull_b10.py を作って。打ち切りは failed=0。
B10、B50、lambda、rho を出し、閉形式と scipy の検算が
小数第1位まで一致するまで直して。一致しない場合は原因を報告して。

「小数第1位まで一致」のように数値で書くと、Claudeは出力を読んで自分で反復できます。テストやスクリプトの出力など、Claudeが会話内で読める合否の信号が反復の足場になります。

実行は python weibull_b10.py の1コマンドに絞れます。settings.jsonのpermissionsでこのコマンドだけを許可しておけば、解析中の確認が増えません。

{
  "permissions": {
    "allow": ["Bash(python weibull_b10.py)"]
  }
}

当てはめと検算のスクリプト

前提として、Pythonにlifelinesとscipyを入れておきます。動かした版はlifelines 0.30.3、scipy 1.18.1、numpy 2.5.4です。

Claude Codeが書くスクリプトの中核は次の部分です。

import numpy as np
import pandas as pd
from lifelines import WeibullFitter
from scipy.optimize import minimize
 
t = np.array([212, 305, 388, 441, 497, 530, 586, 640, 702, 755,
              810, 868, 921, 977, 1000, 1000, 1000, 1000, 1000, 1000], float)
e = (t < 1000).astype(int)          # 1000h で止めた試験。生存は 0
 
wf = WeibullFitter().fit(t, e)
lam, rho = wf.lambda_, wf.rho_
b10 = wf.percentile(0.90)           # S(t)=0.90 となる時間
b10_closed = lam * (-np.log(0.90)) ** (1 / rho)
 
def nll(p):                          # 打ち切り込みの負の対数尤度
    l, r = np.exp(p)
    z = (t / l) ** r
    ll = e * (np.log(r / l) + (r - 1) * np.log(t / l)) - z
    return -ll.sum()
 
res = minimize(nll, [np.log(800), 0.0], method="Nelder-Mead")
print(np.exp(res.x))

故障した個体は密度の対数、打ち切り個体は生存関数の対数だけを尤度に足す、というのが打ち切りの扱いの本体です。上のコードでは e が1の項だけ密度の対数が入り、-z は全個体に入ります。

出力の読み方

実行すると、lifelinesは次のような要約を返しました。

fitted with 20 total observations, 6 right-censored observations
number of events observed = 14
log-likelihood = -106.605
 
           coef  se(coef)  coef lower 95%  coef upper 95%
lambda_ 918.218   105.922         710.614        1125.821
rho_      2.329     0.547           1.257           3.401

最初に見るのは、冒頭の「6 right-censored」が手元の生存個体数と一致しているかです。ここがずれていれば、failed 列の向きを逆に渡しています。

lambdaは「63.2%が壊れた時間」

lambda(尺度)は、人口の63.2%が死亡した時間を表します。S(lambda) はexp(-1)で、スクリプトの出力も0.3679でした。918時間は「約6割が壊れる目安」であって、平均寿命ではありません。

rhoは故障率の形

rho(形状)は累積ハザードが凸か凹かを決め、加速するハザードと減速するハザードを表し分けます。ハザードの式は (rho/lambda)(t/lambda)^(rho-1) なので、rhoが1を超えれば時間とともに故障率が上がります。今回の2.329は摩耗型の挙動で、区間は1.257〜3.401と1をまたぎません。要約表の hypothesis = lambda_ != 1, rho_ != 1 のp値は、rhoが1(指数分布)と違うかを見ています。rhoが1に近い部品では、時間経過で劣化しないため、予防交換の根拠になりません。

B10の検算結果

percentile(0.90) の値と閉形式、scipyの直接最大化は次のとおり一致しました。

• lifelinesの percentile(0.90): 349.4 h
• 閉形式 lambda * (-ln 0.9)^(1/rho): 349.4 h
• scipyの尤度最大化: lambda 918.218、rho 2.329、対数尤度 -106.605

3つが揃えば、failed の向きと打ち切りの扱いは正しいと判断できます。逆に一つだけ外れる場合、その実装の入力の渡し方を疑います。

B10は点推定だけで報告しない

20個で故障が14個だと、B10の不確かさは大きくなります。この例で500回のブートストラップをかけると、B10の95%区間は約234〜510時間でした。349時間という点推定に対して、下限は3分の2です。保証期間を決める用途では、この下限側を見ます。

区間の出し方は、行(個体)を重複を許して20個選び直して当てはめる操作を500回繰り返し、B10の2.5%点と97.5%点を取る方法です。乱数の種を固定しておけば、同じCSVで毎回同じ区間が再現されます。依存のバージョンを固定して再現する流れはClaude CodeでRの解析を進めるが参考になります。lambdaの区間(710.6〜1125.8)のように、lifelinesの要約表が返す区間とは別物なので、どちらの区間かを報告に書き分けます。

Claude Codeには「点推定と区間を併記する」と規約に書いておくのが確実です。書かなければ、点推定だけを表にして終わる出力が返りやすくなります。

区間の幅は故障数でどう変わるか

同じ20個でも、故障の数によってrhoの区間は大きく変わります。故障が早い順に5個・8個だけを故障として残し、残りを1000時間の打ち切りに置き換えた場合の当てはめは次のとおりです(人工的な組み替えで、実際の試験結果ではありません)。

故障5個では区間が1をまたぎ、摩耗型か偶発型かを決められません。それなのにB10の点推定は375時間と、14個のときと近い値が出ます。点推定だけを見ていると、データの薄さに気づけません。Claude Codeの出力に区間を必ず併記させる理由はここにあります。

ノンパラメトリックな曲線と重ねて仮定を確かめる

ワイブルは「この形の曲線になる」という仮定を置いたモデルです。仮定が外れていないかは、仮定を置かないKaplan-Meier推定と重ねて確かめます。lifelinesのガイドも、ワイブルをNelson-Aalen推定のようなパラメータを持たない推定と対比して説明しています。

同じ20個での生存率は次のとおりでした。

どの時間でも差は0.03前後に収まり、ワイブルの当てはめが大きく外れていないと読めます。B10そのものを見ると、Kaplan-Meierの階段では305時間、ワイブルでは349.4時間です。階段は実際に故障した時間の飛び飛びの値しか返せないため、データ点の間を補間できるワイブルとは数十時間の差が出ます。

差が0.1を超える時間帯が続くなら、ワイブルが合っていない可能性があります。その場合は、故障モードが2種類混ざっていないか(初期不良と摩耗など)を疑います。Claude Codeには、この比較表を出力に含めるよう規約に1行足しておけば、毎回自動で並びます。

取り違えたとき何が出るか

failed を逆向きに渡すと、エラーは出ませんが、同じデータでlambdaが1000.0、rhoが3.6×10の24乗というありえない値に張り付きました(lifelines 0.30.3)。故障と打ち切りを入れ替えた入力なので、モデルに合わない当てはめになったと考えられます。rhoが数十を超える値は、まず入力の向きを疑う目安になります。

また、ワイブルのパーセンタイルは percentile(p) が S(t) = p となる時間を返すので、B10は0.90、B1は0.99、中央寿命(B50)は0.50を渡します。この例では順に349.4 h、127.4 h、784.5 hでした。percentile(0.10) と書き間違えると、90%が壊れる時間(B90)を返すので、規約に「B10は0.90」と明記します。

報告に載せる項目をClaudeに固定させる

解析結果を設計や品質の会議に出すとき、数字だけだと条件が伝わりません。出力の形も規約で固定しておくと、読み手が毎回同じ順で確認できます。

• 試験条件: 個体数、故障数、打ち切り数、打ち切り時間
• 推定値: lambda、rhoと、それぞれの95%区間
• 寿命の指標: B10(必要ならB1、B50)と、ブートストラップによる区間
• 検算: lifelines、閉形式、scipyの3値が一致したか
• 仮定の確認: Kaplan-Meierとの生存率の差が最大でいくつか

打ち切りの種類も記録に残します。lifelinesのガイドは、観測できなかった側が上限になる場合を左側打ち切りとして別に扱い、fit_left_censoring を用意しています。試験終了で止めた寿命データは右側打ち切りなので fit で足ります。途中から観測を始めた個体がある場合は、WeibullFitter に開始時間を表す entry の引数があることがAPIページに記載されています。

落とし穴

• 打ち切りの混在: 試験中に外した個体(別原因の故障や取り外し)も打ち切りです。故障モードごとに分けて、対象外の故障は0で渡します。
• 時間の単位: 時間と起動回数が同じ列に混ざると、rhoの意味が崩れます。列名に単位を含めておきます。
• 2母数の仮定: WeibullFitter は2母数モデルです。WeibullFitterのページには、3母数モデルについて別の注記があります。故障が一定時間まったく起きない部品には、位置パラメータが必要かを別に検討します。
• サンプルが少ない: 故障が数個しかないと、rhoの区間が非常に広がります。区間の幅も一緒に読みます。

まとめ

解析を頼む人は、打ち切りの列を先に規約で固定し、次に3値の一致、最後に区間とKaplan-Meierの重ね合わせを見る順で進めます。点推定のB10だけが報告に載っている場合は、区間と検算の有無から確かめるのが近道です。