【Python実践】複雑な相関関係は「次元削減」で解決！予兆検知のための多変量⇒単変量変換アプローチ

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from matplotlib import rcParams
rcParams['font.family'] = 'Meiryo'

# ---------------------------------------------------------
# 1. 4変数の多変量ダミーデータ生成関数（リアルな異常連鎖シナリオ）
# ---------------------------------------------------------
def generate_multivariate_dummy_data(n_samples=500):
    np.random.seed(42)
    
    # ① ベースライン（正常時の相関ルール）
    speed = np.linspace(1000, 3000, n_samples) + np.sin(np.linspace(0, 10, n_samples)) * 500 + np.random.normal(0, 50, n_samples)
    rpm = speed * 1.5 + np.random.normal(0, 50, n_samples)       
    pressure = rpm * 0.02 + np.random.normal(0, 1, n_samples)    
    temp = speed * 0.05 + 20 + np.random.normal(0, 2, n_samples) 

    # ② 【異常の混入】時間差を伴う相関の崩壊
    anomaly_start, anomaly_end = 350, 400
    temp_warning_start = 300  # 温度だけ少し前から予兆が出始める

    # 車速: 変化なし（そのまま順調に走り続けている）

    # 温度: 予兆区間(300)から徐々に上昇し、異常区間(350〜)でさらに上昇傾向になる（摩擦熱の蓄積）
    trend = np.linspace(0, 60, anomaly_end - temp_warning_start)
    temp[temp_warning_start:anomaly_end] += trend
    
    # ★追加：異常区間が終わった後(400〜最後)も、恒久的なダメージとして温度が高いままにする
    temp[anomaly_end:] += 60

    # 回転数: 異常区間で「15%のスリップ（滑り）」を発生させる
    rpm[anomaly_start:anomaly_end] = rpm[anomaly_start:anomaly_end] * 0.85 + np.random.normal(0, 30, 50)
    
    # 圧力: 異常区間(350〜)で高いまま張り付く
    pressure[anomaly_start:anomaly_end] = 90 + np.random.normal(0, 1, 50)

    return pd.DataFrame({'Speed': speed, 'RPM': rpm, 'Temp': temp, 'Pressure': pressure})

# ---------------------------------------------------------
# 2. 共通のグラフ描画関数（2軸で4変数を1つにまとめる）
# ---------------------------------------------------------
def plot_multivariate_result(df, score_col, title, color='purple'):
    fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True, gridspec_kw={'height_ratios': [1.8, 1]})
    
    ax1 = axes[0]
    ax2 = ax1.twinx() 
    
    # 【第1軸 (左)】 駆動力系：Speed と RPM
    l1, = ax1.plot(df['Speed'], label='Speed (速度)', color='gray', alpha=0.7)
    l2, = ax1.plot(df['RPM'], label='RPM (回転数)', color='blue', alpha=0.7)
    ax1.set_ylabel('Speed / RPM')

    # 【第2軸 (右)】 負荷系：Temp と Pressure
    l3, = ax2.plot(df['Temp'], label='Temp (温度)', color='red', alpha=0.7)
    l4, = ax2.plot(df['Pressure'], label='Pressure (圧力)', color='orange', alpha=0.7)
    ax2.set_ylabel('Temp / Pressure')

    lines = [l1, l2, l3, l4]
    labels = [l.get_label() for l in lines]
    ax1.legend(lines, labels, loc='upper left')

    ax1.axvspan(300, 350, color='yellow', alpha=0.1, label='温度上昇の予兆')
    ax1.axvspan(350, 400, color='pink', alpha=0.2, label='完全な異常発生')
    ax1.set_title("Raw Data: 4変数の生データ (速度は正常だが、温度が先行上昇（黄色帯）し、回転数が少し滑る（ピンク帯)")

    # 下段：変換後の異常スコア（単変量）
    axes[1].plot(df[score_col], color=color, linewidth=2)
    axes[1].set_title(title)
    axes[1].axvspan(300, 350, color='yellow', alpha=0.1)
    axes[1].axvspan(350, 400, color='pink', alpha=0.2)
    axes[1].set_xlabel('Time')

    plt.tight_layout()
    plt.show()

# ---------------------------------------------------------
# 4. マハラノビス距離の計算と実行
# ---------------------------------------------------------
def calc_mahalanobis_dist(df):
    df_calc = df.copy()
    scaler = StandardScaler()
    scaled_data = scaler.fit_transform(df_calc)
    
    # 4次元空間での共分散行列とその逆行列を計算
    cov_matrix = np.cov(scaled_data, rowvar=False)
    inv_cov_matrix = np.linalg.inv(cov_matrix)
    mean_distr = np.mean(scaled_data, axis=0)
    
    # 4次元空間の中心（正常な群れ）からのマハラノビス距離の2乗を計算
    m_dist = [mahalanobis(row, mean_distr, inv_cov_matrix)**2 for row in scaled_data]
    
    df_calc['mahal_dist'] = m_dist
    return df_calc

# 実行
df_dummy = generate_multivariate_dummy_data()
df_result2 = calc_mahalanobis_dist(df_dummy)
plot_multivariate_result(
    df_result2, 
    'mahal_dist', 
    "マハラノビス距離：正常な群れから外れる温度上昇(黄色帯)と15%の滑り(ピンク帯)を検知！", 
    color='green'
)

# ---------------------------------------------------------
# 3. PCA再構成誤差の計算と実行
# ---------------------------------------------------------
def calc_pca_recon_error(df):
    df_calc = df.copy()
    scaler = StandardScaler()
    scaled_data = scaler.fit_transform(df_calc)
    
    pca = PCA(n_components=1) 
    pca.fit(scaled_data) 
    
    compressed = pca.transform(scaled_data)
    reconstructed = pca.inverse_transform(compressed)
    
    df_calc['pca_recon_error'] = np.square(scaled_data - reconstructed).sum(axis=1)
    return df_calc

# 実行
df_dummy = generate_multivariate_dummy_data()
df_result1 = calc_pca_recon_error(df_dummy)
plot_multivariate_result(df_result1, 'pca_recon_error', "PCA再構成誤差：目視しづらい温度上昇(黄色帯)と15%の滑り(ピンク帯)の相関崩れを捕捉！")