【Python実践】ワイブル分布（分析）で故障予知！（コピペで使えるサンプルコード付き）

import numpy as np
import pandas as pd

# --- サンプルデータの生成設定 ---
np.random.seed(42)  # 再現性のためのシード設定

# ワイブル分布の真のパラメータを設定 (ここでは摩耗故障傾向のデータを想定)
true_beta = 3.0  # 形状パラメータ (β) - 摩耗故障傾向 (>1)
true_eta = 1000  # 尺度パラメータ (η) - 特性寿命
num_samples = 100 # 生成するデータ数

# ワイブル分布から故障時間を生成
# 注意: np.random.weibull は尺度パラメータを1としているため、
# 生成後に true_eta を掛けてスケールを調整します。
raw_failure_times = np.random.weibull(true_beta, size=num_samples) * true_eta

# 打ち切りデータを導入 (ここでは、1200時間で試験を終了すると仮定)
# 1200時間を超えて故障したものは、1200時間で打ち切りとする
observed_times = np.where(raw_failure_times > 1200, 1200, raw_failure_times)

# 故障したかどうかを示すフラグ
# observed_times == raw_failure_times の場合、実際に故障した（イベント発生）
# それ以外は打ち切り（イベント未発生）
event_observed = (observed_times == raw_failure_times)

# データフレームの作成
df_sample_data = pd.DataFrame({
    'ID': range(1, num_samples + 1), # 製品ID
    'Time': observed_times,          # 観察された時間（故障時間または打ち切り時間）
    'Event': event_observed          # 故障したかどうか (True=故障, False=打ち切り)
})

# CSVファイルとして保存
csv_filename = 'weibull_sample_data.csv'
df_sample_data.to_csv(csv_filename, index=False)

print(f"サンプルデータを '{csv_filename}' に保存しました。")
print("\n生成されたデータの最初の5行:")
print(df_sample_data.head())
print("\n故障データの数:", df_sample_data['Event'].sum())
print("打ち切りデータの数:", (~df_sample_data['Event']).sum())

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import weibull_min
from matplotlib import rcParams
rcParams['font.family'] = 'Meiryo'

def plot_weibull_probability(df, duration_col, event_col, confidence=0.95):
    """
    Weibull確率プロットを描画する関数（理論線・信頼区間付き）
    
    Parameters:
    - df: pandas.DataFrame
        データフレーム（寿命データ or 任意の観測値）
    - duration_col: str
        寿命または観測値の列名
    - event_col: str
        イベントが発生したかどうか（True/False）の列名
    - confidence: float
        信頼区間の信頼水準（例: 0.95）

    Returns:
    - なし（グラフを表示）
    """

    # 故障 or イベント発生のみ取り出し
    data = df[df[event_col]][duration_col].values
    data_sorted = np.sort(data)
    n = len(data_sorted)

    # 経験的累積分布（Bernardの公式）
    F = np.array([(i - 0.3) / (n + 0.4) for i in range(1, n + 1)])
    ln_time = np.log(data_sorted)
    ln_ln_y = np.log(-np.log(1 - F))

    # ワイブル分布フィッティング。flocは位置パラメータ。
    # 多くの寿命データはゼロからスタートする（負にならない）ので、floc=0を指定するのが通例。
    shape, loc, scale = weibull_min.fit(data, floc=0)

    # 理論線と信頼区間の描画準備
    x_theory = np.linspace(data_sorted.min(), data_sorted.max(), 100)
    ln_x_theory = np.log(x_theory)
    F_theory = weibull_min.cdf(x_theory, shape, loc=0, scale=scale)
    ln_ln_F_theory = np.log(-np.log(1 - F_theory))

    # 信頼区間（正規近似）
    z = abs(np.percentile(np.random.normal(0, 1, 100000), [(1 - confidence) / 2 * 100, (1 + confidence) / 2 * 100]))
    stderr = 1 / (shape * np.sqrt(n))
    ln_ln_F_upper = ln_ln_F_theory + z[1] * stderr
    ln_ln_F_lower = ln_ln_F_theory - z[1] * stderr

    # プロット
    plt.figure(figsize=(9, 6))
    plt.plot(ln_time, ln_ln_y, 'o', label='データ点')
    plt.plot(ln_x_theory, ln_ln_F_theory, 'r--', label=f'フィット直線\n(β={shape:.2f}, η={scale:.0f})')
    plt.fill_between(ln_x_theory, ln_ln_F_lower, ln_ln_F_upper, color='r', alpha=0.2, label=f'{int(confidence*100)}% 信頼区間')
    plt.xlabel('ln(Duration)')
    plt.ylabel('ln(-ln(1 - F))')
    plt.title('Weibull Probability Plot')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

# --- 使用例 ---
if __name__ == "__main__":

    # データの読込
    csv_filename = 'weibull_sample_data.csv'
    df = pd.read_csv(csv_filename)

    # ワイブル確率プロット
    plot_weibull_probability(df, duration_col='Time', event_col='Event')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import weibull_min
from matplotlib import rcParams

rcParams['font.family'] = 'Meiryo'

def analyze_weibull_distribution(df, duration_col, event_col):
    """
    データフレーム df の故障時間と故障イベントを使用して、ワイブル分布を推定し、ヒストグラムと比較する

    Parameters:
    df : pandas.DataFrame
        ワイブル分析を行うデータフレーム
    duration_col : str
        故障時間のカラム名
    event_col : str
        故障イベント（True/False）のカラム名
    """
    # 故障イベントが発生したデータのみ抽出
    failure_times = df[df[event_col]][duration_col].values

    # ワイブル分布のパラメータを推定（形状パラメータ β、尺度パラメータ η）
    shape, loc, scale = weibull_min.fit(failure_times, floc=0)
    
    print(f"推定された形状パラメータ β: {shape:.3f}")
    print(f"推定された尺度パラメータ η: {scale:.3f}")

    # ヒストグラム＋フィットしたワイブル分布の描画
    x = np.linspace(0, failure_times.max(), 100)
    pdf_fitted = weibull_min.pdf(x, shape, loc=loc, scale=scale)

    plt.hist(failure_times, bins=15, density=True, alpha=0.6, label='実データ')
    plt.plot(x, pdf_fitted, 'r-', label='フィットしたワイブル分布')
    plt.xlabel(duration_col)
    plt.ylabel('Density')
    plt.title('ヒストグラムとフィットしたワイブル分布')
    plt.legend()
    plt.grid(True)
    plt.show()

# --- 使用例 ---
if __name__ == "__main__":
    # データフレームを作成（実際はCSVを読み込む）
    csv_filename = 'weibull_sample_data.csv'
    df = pd.read_csv(csv_filename)

    # ヒストグラム + PDF の描画
    analyze_weibull_distribution(df, duration_col='Time', event_col='Event')

import numpy as np
import pandas as pd
from scipy.stats import weibull_min
from scipy.special import gamma  # ガンマ関数

def calculate_weibull_parameters(df: pd.DataFrame, duration_col: str, event_col: str) -> dict:
    """
    DataFrameからワイブル分布の形状パラメータβ、尺度パラメータη、中央値、MTTFを推定。

    Parameters:
    - df: pandas.DataFrame
    - duration_col: 寿命や経過時間を表す列名
    - event_col: 故障イベント（True/False）を表す列名

    Returns:
    - dict: {'beta': β, 'eta': η, 'median': 中央値, 'mttf': 平均寿命}
    """

    # 故障したデータ（Event=True）のみを使用
    failure_data = df[df[event_col] == True][duration_col]

    # ワイブル分布にフィット
    beta, loc, eta = weibull_min.fit(failure_data, floc=0)

    # 中央値： η * (ln(2))^(1/β)
    median = eta * (np.log(2)) ** (1 / beta)

    # MTTF（平均寿命）: η * Γ(1 + 1/β)
    mttf = eta * gamma(1 + 1 / beta)

    return {
        'beta': beta,
        'eta': eta,
        'median': median,
        'mttf': mttf
    }

# --- 使用例 ---
if __name__ == "__main__":
    # データフレームを作成（実際はCSVを読み込む）
    csv_filename = 'weibull_sample_data.csv'
    df = pd.read_csv(csv_filename)

    # 
    result = calculate_weibull_parameters(df, duration_col='Time', event_col='Event')
    print(f"形状パラメータ β       : {result['beta']:.3f}")
    print(f"尺度パラメータ η       : {result['eta']:.3f}")
    print(f"中央値（50%が故障） : {result['median']:.3f}")
    print(f"平均寿命（MTTF）     : {result['mttf']:.3f}")

import pandas as pd
import numpy as np
from lifelines import WeibullFitter

def weibull_failure_probability(df, time_start, time_end):
    """
    Weibullモデルを使用して、特定の時間範囲での故障確率を計算する関数
    
    :param df: データフレーム (Time列とEvent列を含む)
    :param time_start: 計算開始時間
    :param time_end: 計算終了時間
    :return: 故障確率（%表記）
    """
    # Weibullモデルのフィッティング
    wf = WeibullFitter()
    wf.fit(durations=df['Time'], event_observed=df['Event'])

    # 推定されたパラメータの表示
    print("推定された形状パラメータ β:", round(wf.rho_, 3))
    print("推定された尺度パラメータ η:", round(wf.lambda_, 3))

    # 生存関数を用いた故障確率の計算
    S_start = wf.survival_function_at_times(time_start).values[0]
    S_end = wf.survival_function_at_times(time_end).values[0]

    failure_prob = S_start - S_end

    return f"部品が {time_start}～{time_end} 時間の間に故障する確率: {failure_prob:.2%}"


if __name__ == "__main__":

    # CSVデータの読み込み
    df = pd.read_csv('weibull_sample_data.csv')
    
    # 例: 500時間～700時間の間に故障する確率を計算
    result = weibull_failure_probability(df, 500, 700)
    print(result)

import pandas as pd
import numpy as np
from lifelines import WeibullFitter

def evaluate_reliability(df, t):
    """
    Weibullモデルを使用して、指定時刻での信頼性評価を行う関数
    
    :param df: データフレーム (Time列とEvent列を含む)
    :param t: 評価したい時刻（時間）
    :return: 信頼度（生存関数）と故障率（ハザード関数）
    """
    # Weibullモデルのフィッティング
    wf = WeibullFitter()
    wf.fit(durations=df['Time'], event_observed=df['Event'])

    # 信頼度（生存関数）： R(t) = P(T > t)
    reliability = wf.survival_function_at_times(float(t)).values[0]

    # 故障率（ハザード関数）： h(t) = f(t) / R(t)
    try:
        hazard_rate = wf.hazard_at_times(float(t)).values[0]
    except KeyError:
        hazard_rate = np.nan  # t が範囲外なら NaN を返す

    # 結果の出力
    print(f"形状パラメータ β: {wf.rho_:.3f}")
    print(f"尺度パラメータ η: {wf.lambda_:.3f}")

    print(f"信頼度 R({t})（t時間まで故障しない確率）: {reliability:.2%}")
    if np.isnan(hazard_rate):
        print(f"故障率 h({t})（t時間時点での故障の激しさ）: データ範囲外")
    else:
        print(f"故障率 h({t})（t時間時点での故障の激しさ）: {hazard_rate:.6f} / 時間")

if __name__ == "__main__":
    # CSVデータの読み込み
    df = pd.read_csv('weibull_sample_data.csv')
    
    # 評価する時刻（例: 800時間）
    evaluate_reliability(df, 800)