【Python実践】LSTMで時系列データ分析！予測・異常検知・補正を攻略

# =========================================================================
# このプログラムは、各章で紹介するプログラムの先頭に張り付けてご使用ください。
# =========================================================================

import pandas as pd
import numpy as np
from sklearn.metrics import mean_squared_error
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout
from tensorflow.keras.callbacks import EarlyStopping
from sklearn.preprocessing import MinMaxScaler
import matplotlib.pyplot as plt
from matplotlib import rcParams
# 日本語フォントの設定
rcParams['font.family'] = 'Meiryo'

def create_and_train_lstm_model(
    df: pd.DataFrame,
    time_col,
    value_col: str,
    sequence_length: int = 15,
    test_ratio: float = 0.2,
    epochs: int = 100,
    batch_size: int = 16
) :
    """
    LSTMモデルを構築・訓練し、テスト結果と学習曲線も表示。
    Args:
        df (pd.DataFrame): 時系列データ（タイムスタンプを含む）
        time_col (str): タイムスタンプ列の名前
        value_col (str): 予測対象の数値データの列名
        sequence_length (int, optional): LSTMが参照する過去のデータ長（デフォルト：15）
        test_ratio (float, optional): 学習データとテストデータの分割比率（デフォルト：0.2）
        epochs (int, optional): モデルの学習回数（デフォルト：100）
        batch_size (int, optional): 1回の学習で処理するデータ数（デフォルト：16）
    Returns:
        model: 学習済みモデル
        test_result_df: テストデータの予測結果
        scaler: MinMaxScaler（逆変換に使用）
        history: 学習履歴（学習曲線描画用）
    """

    # タイムスタンプ列をインデックスに設定して並べ替え
    df = df.copy()
    df.set_index(time_col , inplace=True, drop=True)
    df = df.sort_index()

    # スケーリング
    scaler = MinMaxScaler()
    scaled_values = scaler.fit_transform(df[[value_col]])

    # シーケンス作成
    X, y = [], []
    for i in range(len(scaled_values) - sequence_length):
        X.append(scaled_values[i : i + sequence_length])
        y.append(scaled_values[i + sequence_length])
    X, y = np.array(X), np.array(y)

    # 学習・テスト分割
    split_idx = int(len(X) * (1 - test_ratio))
    X_train, X_test = X[:split_idx], X[split_idx:]
    y_train, y_test = y[:split_idx], y[split_idx:]
    dates = df.index[sequence_length:]

    # モデル定義
    model = Sequential([
        LSTM(50, activation='relu', input_shape=(sequence_length, 1)),
        Dropout(0.2),
        Dense(1)
    ])
    model.compile(optimizer='adam', loss='mse')

    # 学習
    early_stopping = EarlyStopping(monitor='val_loss', patience=10, restore_best_weights=True)
    history = model.fit(
        X_train, y_train,
        epochs=epochs,
        batch_size=batch_size,
        validation_split=0.1,
        callbacks=[early_stopping],
        verbose=0
    )

    # 予測・逆変換
    y_pred_test = model.predict(X_test, verbose=0).flatten()
    y_test_inv = scaler.inverse_transform(y_test.reshape(-1, 1)).flatten()
    y_pred_inv = scaler.inverse_transform(y_pred_test.reshape(-1, 1)).flatten()

    # RMSE計算
    rmse = np.sqrt(mean_squared_error(y_test_inv, y_pred_inv))
    print(f"テストRMSE: {rmse:.4f}")

    # 結果DataFrame
    test_dates = dates[split_idx:]
    test_result_df = pd.DataFrame({
        'date': test_dates,
        'actual': y_test_inv,
        'predicted': y_pred_inv
    })

    return model, test_result_df, scaler, history


def forecast_future(model: Sequential, df: pd.DataFrame, time_col: str, value_col: str,
                    forecast_days: int, sequence_length: int, scaler: MinMaxScaler) -> pd.DataFrame:
    """
    指定した日付以降の未来データをLSTMモデルで予測する。

    Args:
        model (Sequential): 学習済みLSTMモデル。
        time_col (str): タイムスタンプ列の名前
        df (pd.DataFrame): 時系列データ（インデックスは日付）。
        value_col (str): 予測対象の列名。
        forecast_days (int): 予測する日数。
        sequence_length (int): LSTMのシーケンス長。
        scaler (MinMaxScaler): 正規化用スケーラー。

    Returns:
        pd.DataFrame: 予測された未来データ（日付, 予測値）。
    """
    # タイムスタンプ列をインデックスに設定して並べ替え
    df = df.copy()
    df.set_index(time_col , inplace=True, drop=True)
    df = df.sort_index()

    values = df[value_col].values.reshape(-1, 1)
    normalized = scaler.transform(values)

    # 既存データの最後を使用
    future_input = normalized[-sequence_length:].tolist() 

    future_predictions = []
    start_date = df.index.max() + pd.Timedelta(days=1)  # 最後の日付の翌日を予測開始点に設定
    future_dates = pd.date_range(start=start_date, periods=forecast_days, freq="D")

    for _ in range(forecast_days):
        x_input = np.array(future_input[-sequence_length:]).reshape(1, sequence_length, 1)
        y_pred = model.predict(x_input, verbose=0)
        future_predictions.append(y_pred[0][0])
        future_input.append([y_pred[0][0]])

    forecast_values = scaler.inverse_transform(np.array(future_predictions).reshape(-1, 1)).flatten()
    return pd.DataFrame({"date": future_dates, "predicted": forecast_values})

def predict_past(model: Sequential, df: pd.DataFrame, time_col: str, value_col: str,
                 sequence_length: int, scaler: MinMaxScaler) -> pd.DataFrame:
    """
    LSTMモデルを用いて過去データを予測する。

    Args:
        model (Sequential): 学習済みLSTMモデル。
        time_col (str): タイムスタンプ列の名前
        df (pd.DataFrame): 時系列データ。
        value_col (str): 予測対象の列名。
        sequence_length (int): シーケンス長。
        scaler (MinMaxScaler): 正規化スケーラー。

    Returns:
        pd.DataFrame: 実測値と予測値を含むDataFrame（date, actual, predicted）。
    """

    # タイムスタンプ列をインデックスに設定して並べ替え
    df = df.copy()
    df.set_index(time_col , inplace=True, drop=True)
    df = df.sort_index()

    values = df[value_col].values.reshape(-1, 1)
    normalized = scaler.transform(values)

    X = []
    for i in range(len(normalized) - sequence_length):
        X.append(normalized[i:i + sequence_length])
    X = np.array(X).reshape(-1, sequence_length, 1)

    preds = model.predict(X, verbose=0)
    preds_inv = scaler.inverse_transform(preds).flatten()
    actual_inv = values[sequence_length:].flatten()
    dates = df.index[sequence_length:]

    return pd.DataFrame({"date": dates, "actual": actual_inv, "predicted": preds_inv})

def detect_anomalies(pred_df: pd.DataFrame, threshold: float = 2.0) -> pd.DataFrame:
    """
    過去予測結果を用いて異常値を検知する（zスコアで判定）。

    Args:
        pred_df (pd.DataFrame): 実測値と予測値を含むDataFrame。
        threshold (float): zスコアの閾値（通常2.0以上）。

    Returns:
        pd.DataFrame: 異常フラグと誤差列を追加したDataFrame。
    """
    pred_df = pred_df.copy()
    pred_df["error"] = pred_df["actual"] - pred_df["predicted"]
    z_score = (pred_df["error"] - pred_df["error"].mean()) / pred_df["error"].std()
    pred_df["anomaly"] = z_score.abs() > threshold
    return pred_df

def correct_anomalies(anomalies_df: pd.DataFrame) -> pd.DataFrame:
    """
    異常検知結果を基にデータを補正し、新しい列 'corrected' を追加する。

    Args:
        anomalies_df (pd.DataFrame): 異常検知結果（異常値フラグを含む DataFrame）

    Returns:
        pd.DataFrame: 異常値を補正した DataFrame（新しい列 'corrected' を追加）
    """

    corrected_df = anomalies_df.copy()  # 元データをコピー

    # 異常値の場合は予測値 (`predicted`)、そうでない場合は実測値 (`actual`) を補正値としてセット
    corrected_df["corrected"] = np.where(corrected_df["anomaly"], corrected_df["predicted"], corrected_df["actual"])

    return corrected_df

def plot_results(test_df: pd.DataFrame, history) -> None:
    """ 
    テスト予測結果と学習曲線を1枚のシートに上下で表示する関数。

    Args:
        test_df (pd.DataFrame): テストデータ（'date', 'actual', 'predicted'を含む）。
        history: 訓練履歴（Kerasの学習履歴オブジェクト）。

    Returns:
        None: グラフを描画する。
    """
    fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 10))

    # テスト予測の可視化
    axes[0].plot(test_df['date'], test_df['actual'], label='実測値')
    axes[0].plot(test_df['date'], test_df['predicted'], label='予測値', linestyle='--')
    axes[0].set_title("LSTM テスト予測結果")
    axes[0].set_xlabel("日時")
    axes[0].set_ylabel("値")
    axes[0].legend()
    axes[0].grid(True)

    # 学習曲線の可視化
    axes[1].plot(history.history['loss'], label='訓練誤差')
    axes[1].plot(history.history['val_loss'], label='検証誤差')
    axes[1].set_title('学習曲線（損失）')
    axes[1].set_xlabel('エポック')
    axes[1].set_ylabel('MSE')
    axes[1].legend()
    axes[1].grid(True)

    plt.tight_layout()  # レイアウト調整
    plt.show()

if __name__ == "__main__":

    # ダミーデータの読込
    df = pd.read_csv("dummy_timeseries.csv", parse_dates=["timestamp"])

    # LSTMモデルの作成と学習
    model, test_df, scaler, history = create_and_train_lstm_model(
        df, time_col="timestamp",value_col="value", sequence_length=15
    )

    # 学習結果のグラフ化
    plot_results(test_df,history)

#========================================================================
# ここに、「LSTMモデルの作成」で紹介したプログラムを張り付けてください。
#========================================================================

def plot_lstm_predictions(df: pd.DataFrame, past_df: pd.DataFrame, future_df: pd.DataFrame, value_col: str) -> None:
    """ 
    LSTMの過去予測・未来予測を統合してプロットする関数。

    Args:
        df (pd.DataFrame): 元の時系列データ（実測値）
        past_df (pd.DataFrame): LSTMによる過去の予測結果
        future_df (pd.DataFrame): LSTMによる未来予測結果
        value_col (str): 予測対象の列名

    Returns:
        None: グラフを描画する。
    """
    # 日付データの型を統一
    past_df["date"] = pd.to_datetime(past_df["date"])
    future_df["date"] = pd.to_datetime(future_df["date"])

    plt.figure(figsize=(14, 7))

    # 過去の測定値（青色）
    plt.plot(df["timestamp"], df[value_col], label="過去の測定値 (学習用)", color="blue", alpha=0.9)

    # LSTM学習データの予測（緑の破線）
    plt.plot(past_df["date"], past_df["predicted"], label="LSTM予測 (学習データ)", color="green", linestyle="--", alpha=0.8)

    # 未来のLSTM予測（赤の破線）
    plt.plot(future_df["date"], future_df["predicted"], label="LSTM未来予測", color="red", linestyle="-.", alpha=0.8)

    # 予測開始点（縦の点線）
    prediction_start = past_df["date"].iloc[-1]
    plt.axvline(x=prediction_start, color="gray", linestyle=":", label="予測開始時点")

    # X軸の範囲を設定（過去の異常値を排除）
    plt.xlim(past_df["date"].min(), future_df["date"].max())

    plt.title("LSTMによる未来予測")
    plt.xlabel("日付")
    plt.ylabel("測定値")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":

    # ダミーデータの読込
    df = pd.read_csv("dummy_timeseries.csv", parse_dates=["timestamp"])

    # LSTMモデルの作成と学習
    model, test_df, scaler, history = create_and_train_lstm_model(
        df, time_col="timestamp",value_col="value", sequence_length=15
    )
    # 過去予測
    past_df = predict_past(model=model, df=df, time_col="timestamp",value_col="value", sequence_length=15, scaler=scaler)
    print(past_df)

    # 未来予測
    future_df = forecast_future(model=model, df=df, time_col="timestamp",value_col="value", forecast_days=14, sequence_length=15, scaler=scaler)
    print(future_df) 

    # 未来予測グラフの表示
    plot_lstm_predictions(df, past_df, future_df, value_col="value")

#========================================================================
# ここに、「LSTMモデルの作成」で紹介したプログラムを張り付けてください。
#========================================================================

def plot_anomalies_and_errors(df: pd.DataFrame, anomalies_df: pd.DataFrame) -> None:
    """ 
    LSTMモデルによる異常検知と予測誤差をプロットする関数。

    Args:
        df (pd.DataFrame): 実測値と予測値を含む元データ
        anomalies_df (pd.DataFrame): 異常値を含むDataFrame

    Returns:
        None: グラフを描画する。
    """
    plt.figure(figsize=(14, 12))

    # 上段：測定値・予測値・異常値のプロット
    plt.subplot(2, 1, 1)
    plt.plot(df["date"], df["actual"], label="実測値", color="blue", alpha=0.9)
    plt.plot(df["date"], df["predicted"], label="LSTM予測", color="orange", linestyle="--", alpha=0.8)
    
    # 異常値（赤点）
    anomalies = anomalies_df[anomalies_df["anomaly"]]
    plt.scatter(anomalies["date"], anomalies["actual"], label="異常値 (LSTM検知に基づく)", color="red", marker="o")

    plt.title("LSTMモデルによる時系列異常検知（測定値と予測）")
    plt.xlabel("日付")
    plt.ylabel("測定値")
    plt.legend()
    plt.grid(True)

    # 下段：予測誤差と異常検知閾値
    plt.subplot(2, 1, 2)
    plt.plot(df["date"], df["actual"] - df["predicted"], label="予測誤差", color="purple", alpha=0.8)
    
    # 誤差の平均値
    mean_error = (df["actual"] - df["predicted"]).mean()
    plt.axhline(mean_error, color="black", linestyle="--", label="誤差平均")

    # ±3標準偏差閾値
    std_error = (df["actual"] - df["predicted"]).std()
    plt.axhline(mean_error + 3 * std_error, color="orange", linestyle="--", label="+3.0標準偏差")
    plt.axhline(mean_error - 3 * std_error, color="orange", linestyle="--", label="-3.0標準偏差")

    # 異常値（赤点）
    plt.scatter(anomalies["date"], anomalies["actual"] - anomalies["predicted"], label="異常値", color="red", marker="o")

    plt.title("LSTMモデル予測誤差と異常検知閾値")
    plt.xlabel("日付")
    plt.ylabel("予測誤差")
    plt.legend()
    plt.grid(True)

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":

    # ダミーデータの読込
    df = pd.read_csv("dummy_timeseries.csv", parse_dates=["timestamp"])

    # LSTMモデルの作成と学習
    model, test_df, scaler, history = create_and_train_lstm_model(
        df, time_col="timestamp",value_col="value", sequence_length=15
    )
    # 過去予測
    past_df = predict_past(model=model, df=df, time_col="timestamp",value_col="value", sequence_length=15, scaler=scaler)

    # 異常検知
    anomalies_df = detect_anomalies(past_df,3)
    print(anomalies_df)

    # 異常検知結果と誤差のグラフ化
    plot_anomalies_and_errors(past_df,anomalies_df)

#========================================================================
# ここに、「LSTMモデルの作成」で紹介したプログラムを張り付けてください。
#========================================================================
def plot_anomalies(corrected_df: pd.DataFrame, threshold: float = None) -> None:
    """ 
    異常値を検知し、補正後の値を比較するグラフを描画する関数（閾値 or 信頼区間付き）。

    Args:
        corrected_df (pd.DataFrame): 異常検知済みデータ（補正値は描画しない）
        threshold (float, optional): 異常検知に使用するσ値（例: 2.0, 3.0, 3.5）。
                                     Noneの場合は信頼区間（±1.96σ）を表示。

    Returns:
        None: グラフを描画する。
    """

    plt.figure(figsize=(14, 7))

    # 実測値（青線）
    plt.plot(corrected_df["date"], corrected_df["actual"], label="元データ", color="blue", alpha=0.7)

    # LSTM予測値（緑の破線）
    plt.plot(corrected_df["date"], corrected_df["predicted"], label="LSTM予測値", color="green", linestyle="--", alpha=0.8)

    if threshold is not None:
        # 🔹 異常検知閾値（±thresholdσ）
        std_error = corrected_df["error"].std()
        mean_error = corrected_df["error"].mean()
        lower_bound = corrected_df["predicted"] + (mean_error - threshold * std_error)
        upper_bound = corrected_df["predicted"] + (mean_error + threshold * std_error)

        label_text = f"異常検知閾値 (±{threshold}σ)"
    else:
        # 🔹 信頼区間（±1.96σ, 95% CI）
        std_predicted = corrected_df["predicted"].std()
        lower_bound = corrected_df["predicted"] - 1.96 * std_predicted
        upper_bound = corrected_df["predicted"] + 1.96 * std_predicted

        label_text = "95% 信頼区間 (±1.96σ)"

    # 閾値 or 信頼区間をグレーの領域で表示
    plt.fill_between(corrected_df["date"], lower_bound, upper_bound, color="gray", alpha=0.2, label=label_text)

    # 異常値（赤点）
    anomalies = corrected_df[corrected_df["anomaly"]]
    plt.scatter(anomalies["date"], anomalies["actual"], label="異常値",
                color="red", marker="o", zorder=3)

    # 異常値の補正結果（緑点）
    plt.scatter(anomalies["date"], anomalies["predicted"], label="補正後の異常値",
                color="limegreen", marker="o", edgecolors="darkgreen", linewidths=1.5, zorder=3)

    plt.title(f"LSTMによる異常検知と {'閾値' if threshold else '信頼区間'} の範囲")
    plt.xlabel("日付")
    plt.ylabel("測定値")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":

    # ダミーデータの読込
    df = pd.read_csv("dummy_timeseries.csv", parse_dates=["timestamp"])

    # LSTMモデルの作成と学習
    model, test_df, scaler, history = create_and_train_lstm_model(
        df, time_col="timestamp",value_col="value", sequence_length=15
    )
    # 過去予測
    past_df = predict_past(model=model, df=df, time_col="timestamp",value_col="value", sequence_length=15, scaler=scaler)

    # 異常検知
    anomalies_df = detect_anomalies(past_df,3)

    # データ補正
    correct_df = correct_anomalies(anomalies_df)
    print(correct_df)

    # 異常検知結果と補正結果のグラフ
    plot_anomalies(correct_df,3)

   # 異常検知結果と補正結果のグラフ
   plot_anomalies(correct_df)