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

import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
from prophet import Prophet
import warnings

warnings.filterwarnings("ignore")  # 警告を無視

# 日本語フォントの設定
rcParams['font.family'] = 'Meiryo'

def prophet_forecast_future(df: pd.DataFrame, timestamp_col: str, value_col: str, forecast_steps: int = 14):
    """
    Prophetを使用して未来予測を行う。

    Parameters:
    df (pd.DataFrame): 時系列データのDataFrame。
    timestamp_col (str): タイムスタンプの列名。
    value_col (str): 対象の数値データの列名。
    forecast_steps (int): 未来の予測ステップ（日数単位）。

    Returns:
    pd.DataFrame: 予測結果（未来の予測値を含む）。
    """

    # Prophet用に列名を変更
    df_processed = df.rename(columns={timestamp_col: "ds", value_col: "y"}).copy()

    # 欠損値を線形補完
    df_processed["y"] = df_processed["y"].interpolate(method="linear").fillna(method="bfill").fillna(method="ffill")

    # モデル定義・学習
    print(f"\n--- Prophetによる未来予測 ({forecast_steps}日先まで) ---")
    model = Prophet(daily_seasonality=True)
    model.fit(df_processed)
    print("Prophetモデルの学習が完了しました。")

    # 未来の期間データ作成
    future = model.make_future_dataframe(periods=forecast_steps)
    forecast = model.predict(future)

    # フィッティング範囲と未来予測範囲に分割
    past_forecast = forecast.iloc[:len(df_processed)]
    future_forecast = forecast.iloc[len(df_processed):]

    # 結果のプロット
    plt.figure(figsize=(14, 7))

    # 過去の実測値（青実線）
    plt.plot(df_processed["ds"], df_processed["y"], label="過去の測定値 (学習用)", color="blue")

    # 学習データに対する予測（緑点線）
    plt.plot(past_forecast["ds"], past_forecast["yhat"], label="Prophet予測 (学習データ)", color="green", linestyle="--")

    # 未来予測（赤点線）
    plt.plot(future_forecast["ds"], future_forecast["yhat"], label="Prophet未来予測", color="red", linestyle="-.")

    # 予測開始点の縦点線
    forecast_start = df_processed["ds"].max()
    plt.axvline(x=forecast_start, color="gray", linestyle=":", label="予測開始点")

    # グラフ体裁
    plt.xlabel("日付")
    plt.ylabel("測定値")
    plt.title("Prophetによる未来予測")
    plt.legend()
    plt.grid(True)
    plt.show()

    # 未来予測だけを返す
    forecast_result = future_forecast[["ds", "yhat"]].set_index("ds")

    print("予測結果:")
    print(forecast_result)

    return forecast_result

if __name__ == "__main__":
    # データ読み込み
    csv_filename = "dummy_timeseries.csv"
    df = pd.read_csv(csv_filename, parse_dates=["timestamp"])

    # 未来予測の実行
    prophet_predictions = prophet_forecast_future(df, "timestamp", "value", forecast_steps=14)

# pip install Prophet

import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
from prophet import Prophet
import warnings

warnings.filterwarnings("ignore") # 警告メッセージを非表示にする

# 日本語フォントの設定 (お使いの環境に合わせて適宜変更してください)
rcParams['font.family'] = 'Meiryo' # Windowsの場合
# rcParams['font.family'] = 'AppleGothic' # Macの場合
# rcParams['font.family'] = 'DejaVu Sans' # デフォルト、日本語対応しない場合がある

def prophet_anomaly_detection(df: pd.DataFrame, time_col: str, value_col: str, 
                              interval_width: float = 0.95, 
                              yearly_seasonality: bool = False, 
                              weekly_seasonality: bool = True, 
                              daily_seasonality: bool = False):
    """
    Prophetモデルを使用して時系列データの異常検知を行います。
    予測区間からの逸脱を異常と見なします。

    Args:
        df (pd.DataFrame): 処理するデータフレーム。Prophetの要件に合わせて、
                           タイムスタンプ列は 'ds'、値の列は 'y' に内部でリネームされます。
        time_col (str): タイムスタンプの元の列名。
        value_col (str): 時系列データの値の元の列名。
        interval_width (float, optional): 予測区間の幅（信頼水準）。デフォルトは0.95 (95%)。
                                        この範囲外が異常となる。
        yearly_seasonality (bool, optional): 年次季節性を考慮するか。デフォルトはFalse。
        weekly_seasonality (bool, optional): 週次季節性を考慮するか。デフォルトはTrue。
        daily_seasonality (bool, optional): 日次季節性を考慮するか。デフォルトはFalse。

    Returns:
        pd.DataFrame: 'ds', 'y', 'yhat', 'yhat_lower', 'yhat_upper', 'is_anomaly', 'residuals' 列が追加されたデータフレーム。
    """
    
    # Prophetの要件に合わせて列名をリネーム
    df_prophet = df.rename(columns={time_col: 'ds', value_col: 'y'})
    
    # Prophetモデルの構築と学習
    print(f"--- Prophetモデルの適用 (予測区間幅: {interval_width*100}%) ---")
    try:
        model = Prophet(
            interval_width=interval_width,
            yearly_seasonality=yearly_seasonality,
            weekly_seasonality=weekly_seasonality,
            daily_seasonality=daily_seasonality
        )
        # 祝日がある場合は model.add_country_holidays(country_name='JP') などで追加可能
        
        model.fit(df_prophet)
        print("Prophetモデルの学習が完了しました。")
    except Exception as e:
        print(f"Prophetモデルの学習中にエラーが発生しました: {e}")
        return None

    # 学習データ期間内での予測（fitted valuesに相当）
    forecast = model.predict(df_prophet)

    # 予測結果を元のデータフレームに結合
    df_result = df_prophet.set_index('ds').copy()
    df_result = df_result.join(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].set_index('ds'))
    
    # 残差（予測誤差）を計算
    df_result['residuals'] = df_result['y'] - df_result['yhat']

    # 異常検知 (予測区間からの逸脱)
    # 実際の値が予測区間 (yhat_lower, yhat_upper) の外にある場合を異常とする
    df_result['is_anomaly'] = (df_result['y'] < df_result['yhat_lower']) | \
                              (df_result['y'] > df_result['yhat_upper'])

    # 検出された異常値の表示
    anomalies = df_result[df_result['is_anomaly']]
    if not anomalies.empty:
        print("\n--- 検出された異常値 ---")
        # 元の列名に戻して表示 (オプション)
        print(anomalies.rename(columns={'ds': time_col, 'y': value_col})[[value_col, 'yhat', 'yhat_lower', 'yhat_upper', 'residuals', 'is_anomaly']])
    else:
        print("\n--- 異常値は検出されませんでした ---")

    # --- 結果の可視化 (グラフを縦に並べる) ---
    fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(14, 10), sharex=True) # 2行1列でグラフを配置

    # 1つ目のサブプロット: 測定値、Prophet予測、予測区間
    axes[0].plot(df_result.index, df_result['y'], label='測定値', color='blue', alpha=0.7)
    axes[0].plot(df_result.index, df_result['yhat'], label='Prophet予測', color='green', linestyle='--')
    
    # 予測区間のプロット
    axes[0].fill_between(df_result.index, df_result['yhat_lower'], df_result['yhat_upper'], 
                         color='gray', alpha=0.2, label=f'{int(interval_width*100)}% 予測区間')
    
    if not anomalies.empty:
        axes[0].scatter(anomalies.index, anomalies['y'], color='red', s=100, zorder=5, label='異常値 (予測区間外)')

    axes[0].set_title('Prophetモデルによる時系列異常検知 (測定値と予測区間)')
    axes[0].set_ylabel('測定値')
    axes[0].legend()
    axes[0].grid(True)

    # 2つ目のサブプロット: 残差と閾値 (Prophetは残差を直接計算しないため、ここでは「予測誤差」を示す)
    # 修正点: anomalies['residuals'] ではなく df_result['residuals'] を参照する
    axes[1].plot(df_result.index, df_result['residuals'], label='予測誤差', color='purple', alpha=0.7)
    axes[1].axhline(0, color='gray', linestyle='--', label='予測誤差ゼロ')
    
    # 予測区間の上下限を誤差の視点から描画
    axes[1].fill_between(df_result.index, 
                         df_result['yhat_lower'] - df_result['yhat'], 
                         df_result['yhat_upper'] - df_result['yhat'], 
                         color='gray', alpha=0.2, label='予測誤差の許容範囲')

    if not anomalies.empty:
        # ここを修正: anomalies['residuals'] -> df_result.loc[anomalies.index, 'residuals']
        # または、anomalies データフレームに 'residuals' 列がある前提ならそのまま使えるが、
        # より安全に、anomalies['y'] - anomalies['yhat'] を使う
        axes[1].scatter(anomalies.index, anomalies['y'] - anomalies['yhat'], color='red', s=100, zorder=5, label='異常誤差')

    axes[1].set_title('Prophetモデル予測誤差と異常検知範囲')
    axes[1].set_xlabel('日付')
    axes[1].set_ylabel('予測誤差 (実測値 - 予測値)')
    axes[1].legend()
    axes[1].grid(True)

    plt.tight_layout() # サブプロット間のスペースを自動調整
    plt.show()

    return df_result

# --- メインの実行部分 ---
if __name__ == "__main__":
    # CSVを読み込んで関数を実行
    df = pd.read_csv("dummy_timeseries.csv", parse_dates=["timestamp"])

    result_df = prophet_anomaly_detection(df, 
                                          time_col="timestamp", 
                                          value_col="value", 
                                          interval_width=0.99, # 99%予測区間を閾値とする (より厳しく)
                                          yearly_seasonality=False, 
                                          weekly_seasonality=True,  
                                          daily_seasonality=False)
    
    if result_df is not None:
        print("\n--- 関数実行結果 (異常値検出フラグ含む) ---")
        print(result_df.tail())

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
from prophet import Prophet
import warnings

warnings.filterwarnings("ignore")
rcParams['font.family'] = 'Meiryo'

def prophet_anomaly_detection_and_cleaning(df: pd.DataFrame,
                                           time_col: str,
                                           value_col: str,
                                           interval_width: float = 0.95, # New parameter
                                           correction_method: str = "mean"):
    """
    Prophetを使って異常検知＆補正を行う関数

    Parameters:
        df (pd.DataFrame): 時系列データ
        time_col (str): タイムスタンプ列名
        value_col (str): 測定値の列名
        interval_width (float): 予測区間の幅（信頼水準）。0から1の間の値。デフォルトは0.95 (95%)。
        correction_method (str): "mean"（Prophet予測）または "median"（移動中央値）

    Returns:
        pd.DataFrame: 異常値を補正したデータフレーム
    """

    # Prophet 用に整形
    df_prophet = df[[time_col, value_col]].rename(columns={time_col: "ds", value_col: "y"}).copy()
    df_prophet["y"] = df_prophet["y"].interpolate(method="linear").fillna(method="bfill").fillna(method="ffill")

    # モデル構築と予測
    # interval_width を Prophet モデルに渡す
    model = Prophet(daily_seasonality=True, interval_width=interval_width)
    model.fit(df_prophet)
    forecast = model.predict(df_prophet)

    # 異常値の検出: 実際の値が予測区間外にあるかを判定
    df_prophet["yhat"] = forecast["yhat"]
    df_prophet["yhat_lower"] = forecast["yhat_lower"]
    df_prophet["yhat_upper"] = forecast["yhat_upper"]

    df_prophet["is_anomaly"] = (df_prophet["y"] < df_prophet["yhat_lower"]) | \
                               (df_prophet["y"] > df_prophet["yhat_upper"])

    # 異常値補正
    if correction_method == "mean":
        df_prophet["corrected"] = df_prophet.apply(
            lambda row: row["yhat"] if row["is_anomaly"] else row["y"], axis=1
        )
    elif correction_method == "median":
        rolling_median = df_prophet["y"].rolling(window=7, center=True, min_periods=1).median()
        df_prophet["corrected"] = df_prophet.apply(
            lambda row: rolling_median[row.name] if row["is_anomaly"] else row["y"], axis=1
        )
    else:
        raise ValueError("correction_method は 'mean' または 'median' のいずれかを指定してください。")

    # 元のタイムスタンプに戻す
    df_result = df.copy()
    df_result["補正後"] = df_prophet["corrected"].values
    df_result["異常値"] = df_prophet["is_anomaly"].values
    df_result["予測値"] = df_prophet["yhat"].values
    df_result["予測下限"] = df_prophet["yhat_lower"].values
    df_result["予測上限"] = df_prophet["yhat_upper"].values

    # 可視化
    plt.figure(figsize=(14, 7))
    plt.plot(df_prophet["ds"], df_prophet["y"], label="元データ", color="blue")
    plt.plot(df_prophet["ds"], df_prophet["yhat"], label="Prophet予測値", color="green", linestyle="--")
    plt.fill_between(df_prophet["ds"], df_prophet["yhat_lower"], df_prophet["yhat_upper"],
                     color='lightgreen', alpha=0.3, label=f'{int(interval_width*100)}% 予測区間')
    plt.scatter(df_prophet["ds"][df_prophet["is_anomaly"]],
                df_prophet["y"][df_prophet["is_anomaly"]],
                color="red", label="異常値（補正前）", edgecolors="darkred", linewidths=1.5)
    plt.scatter(df_prophet["ds"][df_prophet["is_anomaly"]],
                df_prophet["corrected"][df_prophet["is_anomaly"]],
                color="limegreen", label="補正値", edgecolors="darkgreen", linewidths=1.5)

    plt.xlabel("日付")
    plt.ylabel("測定値")
    plt.title(f"Prophetによる異常検知と補正 (信頼水準: {int(interval_width*100)}%)")
    plt.legend()
    plt.grid(True)
    plt.show()

    return df_result

# --- 実行例 ---
if __name__ == "__main__":
    # Assuming 'dummy_timeseries.csv' exists and has 'timestamp' and 'value' columns
    # For demonstration, let's create a dummy CSV if it doesn't exist
    try:
        df = pd.read_csv("dummy_timeseries.csv", parse_dates=["timestamp"])
    except FileNotFoundError:
        print("dummy_timeseries.csv not found. Creating a sample dummy_timeseries.csv.")
        dates = pd.to_datetime(pd.date_range(start="2023-01-01", periods=100, freq="D"))
        values = np.sin(np.linspace(0, 20, 100)) * 10 + 50 + np.random.randn(100) * 5
        # Introduce some anomalies
        values[10] = 100
        values[50] = 10
        values[80] = 120
        df = pd.DataFrame({"timestamp": dates, "value": values})
        df.to_csv("dummy_timeseries.csv", index=False)
        print("dummy_timeseries.csv created.")
        df = pd.read_csv("dummy_timeseries.csv", parse_dates=["timestamp"])


    cleaned_df = prophet_anomaly_detection_and_cleaning(df,
                                                        time_col="timestamp",
                                                        value_col="value",
                                                        interval_width=0.99, # 信頼水準を99%に設定
                                                        correction_method="mean")

    print("\n--- 補正後のデータ ---")
    print(cleaned_df.tail())