# ============================================ # 指数平滑法による売上予測 # ============================================ # 指数平滑法とは？ # → 最新データに大きな重み、過去ほど小さな重み # → 移動平均法の進化版 # → トレンドや季節性も考慮可能 import pandas as pd import numpy as np import matplotlib.pyplot as plt from statsmodels.tsa.holtwinters import ExponentialSmoothing import warnings warnings.filterwarnings(‘ignore’) # ============================================ # サンプルデータ作成(月次売上 – 3年分) # ============================================ np.random.seed(42) # pd.date_range(): 日付の連続データを生成 # ‘2021-01-01’: 開始日 # ‘2024-12-01′: 終了日 # freq=’MS’: 月初(Month Start) months = pd.date_range(‘2021-01-01’, ‘2024-12-01′, freq=’MS’) # 時系列の3要素を作成 # np.linspace(): 等間隔の数値を生成 # 100から160まで、len(months)個の数値 trend = np.linspace(100, 160, len(months)) # 上昇トレンド # np.sin(): サイン波で季節パターンを作成 # 2 * np.pi / 12: 12ヶ月で1周期 # 25 *: 振幅±25の変動 seasonality = 25 * np.sin(np.arange(len(months)) * 2 * np.pi / 12) # np.random.normal(): 正規分布に従う乱数 # 平均0、標準偏差5のノイズ noise = np.random.normal(0, 5, len(months)) # 3要素を合成 sales = trend + seasonality + noise # データフレーム作成 df = pd.DataFrame({ ‘date’: months, ‘sales’: sales }) # ============================================ # データを訓練とテストに分割 # ============================================ # 85%を訓練用、15%をテスト用 train_size = int(len(df) * 0.85) train = df[:train_size] test = df[train_size:] print(“【データ分割】”) print(f”訓練データ: {len(train)}ヶ月”) print(f”テストデータ: {len(test)}ヶ月”) print() # ============================================ # 精度評価関数 # ============================================ def calculate_metrics(actual, predicted): “”” 予測精度を計算 Parameters: ———– actual : array – 実績値 predicted : array – 予測値 Returns: ——– tuple – (MAE, RMSE, MAPE) “”” mae = np.mean(np.abs(actual – predicted)) rmse = np.sqrt(np.mean((actual – predicted)**2)) mape = np.mean(np.abs((actual – predicted) / actual)) * 100 return mae, rmse, mape # ============================================ # 1. 単純指数平滑法(SES) # ============================================ # ExponentialSmoothing()のパラメータ: # trend=None: トレンド成分なし # seasonal=None: 季節性成分なし model_ses = ExponentialSmoothing( train[‘sales’], trend=None, seasonal=None ).fit() # .forecast(): 未来を予測 # steps: 予測する期間数 forecast_ses = model_ses.forecast(steps=len(test)) # ============================================ # 2. ダブル指数平滑法(Holt法) # ============================================ # trend=’add’: 加法トレンド(線形成長) # → トレンドを考慮するが季節性は考慮しない model_holt = ExponentialSmoothing( train[‘sales’], trend=’add’, seasonal=None ).fit() forecast_holt = model_holt.forecast(steps=len(test)) # ============================================ # 3. トリプル指数平滑法(Holt-Winters法) # ============================================ # trend=’add’: 加法トレンド # seasonal=’add’: 加法季節性 # seasonal_periods=12: 12ヶ月周期(年次季節性) model_hw = ExponentialSmoothing( train[‘sales’], trend=’add’, seasonal=’add’, seasonal_periods=12 ).fit() forecast_hw = model_hw.forecast(steps=len(test)) # ============================================ # 精度比較 # ============================================ print(“【予測精度比較】”) print() mae_ses, rmse_ses, mape_ses = calculate_metrics( test[‘sales’].values, forecast_ses ) print(“単純指数平滑法(SES):”) print(f” MAE: {mae_ses:.2f}”) print(f” RMSE: {rmse_ses:.2f}”) print(f” MAPE: {mape_ses:.2f}%”) print() mae_holt, rmse_holt, mape_holt = calculate_metrics( test[‘sales’].values, forecast_holt ) print(“ダブル指数平滑法(Holt):”) print(f” MAE: {mae_holt:.2f}”) print(f” RMSE: {rmse_holt:.2f}”) print(f” MAPE: {mape_holt:.2f}%”) print() mae_hw, rmse_hw, mape_hw = calculate_metrics( test[‘sales’].values, forecast_hw ) print(“トリプル指数平滑法(Holt-Winters):”) print(f” MAE: {mae_hw:.2f}”) print(f” RMSE: {rmse_hw:.2f}”) print(f” MAPE: {mape_hw:.2f}%”) print() # ============================================ # 最適モデルの選択 # ============================================ # min(): 最小値を持つキーを取得 # key=models.get: 値(MAPE)で比較 models = { ‘SES’: mape_ses, ‘Holt’: mape_holt, ‘Holt-Winters’: mape_hw } best_model = min(models, key=models.get) print(f”【最適モデル】{best_model} (MAPE: {models[best_model]:.2f}%)”) print() # ============================================ # 未来の予測(次の6ヶ月) # ============================================ future_forecast = model_hw.forecast(steps=6) # pd.DateOffset(): 日付をずらす # months=1: 1ヶ月後 future_dates = pd.date_range( df[‘date’].iloc[-1] + pd.DateOffset(months=1), periods=6, freq=’MS’ ) print(“【未来予測(次の6ヶ月)】”) for date, sales in zip(future_dates, future_forecast): print(f”{date.strftime(‘%Y年%m月’)}: {sales:.1f}百万円”)

# ============================================ # 可視化 # ============================================ # plt.subplots(): 複数のグラフを配置 # 2, 2: 2行2列(計4グラフ) # figsize: 図全体のサイズ fig, axes = plt.subplots(2, 2, figsize=(16, 10)) # ============================================ # グラフ1: 単純指数平滑法 # ============================================ ax1 = axes[0, 0] ax1.plot(train[‘date’], train[‘sales’], label=’訓練データ’, color=’blue’) ax1.plot(test[‘date’], test[‘sales’], label=’実績’, color=’green’, marker=’o’) ax1.plot(test[‘date’], forecast_ses, label=’予測’, color=’red’, linestyle=’–‘, marker=’s’) ax1.set_title(f’単純指数平滑法(SES)\nMAPE: {mape_ses:.1f}%’, fontsize=13, fontweight=’bold’) ax1.set_ylabel(‘売上(百万円)’, fontsize=11) ax1.legend() ax1.grid(alpha=0.3) # ============================================ # グラフ2: ダブル指数平滑法 # ============================================ ax2 = axes[0, 1] ax2.plot(train[‘date’], train[‘sales’], label=’訓練データ’, color=’blue’) ax2.plot(test[‘date’], test[‘sales’], label=’実績’, color=’green’, marker=’o’) ax2.plot(test[‘date’], forecast_holt, label=’予測’, color=’red’, linestyle=’–‘, marker=’s’) ax2.set_title(f’ダブル指数平滑法(Holt)\nMAPE: {mape_holt:.1f}%’, fontsize=13, fontweight=’bold’) ax2.set_ylabel(‘売上(百万円)’, fontsize=11) ax2.legend() ax2.grid(alpha=0.3) # ============================================ # グラフ3: トリプル指数平滑法 # ============================================ ax3 = axes[1, 0] ax3.plot(train[‘date’], train[‘sales’], label=’訓練データ’, color=’blue’) ax3.plot(test[‘date’], test[‘sales’], label=’実績’, color=’green’, marker=’o’) ax3.plot(test[‘date’], forecast_hw, label=’予測’, color=’red’, linestyle=’–‘, marker=’s’) ax3.set_title(f’トリプル指数平滑法(Holt-Winters)\nMAPE: {mape_hw:.1f}%’, fontsize=13, fontweight=’bold’) ax3.set_xlabel(‘日付’, fontsize=11) ax3.set_ylabel(‘売上(百万円)’, fontsize=11) ax3.legend() ax3.grid(alpha=0.3) # ============================================ # グラフ4: 精度比較(棒グラフ) # ============================================ ax4 = axes[1, 1] methods = [‘SES’, ‘Holt’, ‘Holt-Winters’] mapes = [mape_ses, mape_holt, mape_hw] colors = [‘#ff9999’, ‘#66b3ff’, ‘#99ff99′] bars = ax4.bar(methods, mapes, color=colors, alpha=0.7, edgecolor=’black’) ax4.set_ylabel(‘MAPE (%)’, fontsize=11) ax4.set_title(‘予測精度比較(MAPE)’, fontsize=13, fontweight=’bold’) ax4.grid(axis=’y’, alpha=0.3) # 値を棒グラフ上に表示 for bar, mape in zip(bars, mapes): height = bar.get_height() ax4.text(bar.get_x() + bar.get_width()/2., height, f'{mape:.1f}%’, ha=’center’, va=’bottom’, fontsize=11) plt.tight_layout() plt.show()

# ============================================ # ARIMAモデルによる予測 # ============================================ # ARIMA = AutoRegressive Integrated Moving Average # AR(p): 過去の値を使った自己回帰 # I(d): 差分(トレンド除去) # MA(q): 過去の誤差を使った移動平均 from statsmodels.tsa.arima.model import ARIMA from statsmodels.graphics.tsaplots import plot_acf, plot_pacf # データ準備(同じデータを使用) train_arima = train[‘sales’].values test_arima = test[‘sales’].values # ============================================ # 複数のARIMAモデルを試す # ============================================ # 様々な(p, d, q)の組み合わせでAICを比較 arima_configs = [ (1, 1, 1), # 基本形 (0, 1, 1), # シンプル (2, 1, 2), # 複雑 (1, 1, 0), # MA項なし (0, 1, 2) # AR項なし ] print(“【ARIMAモデル比較】”) print() best_aic = float(‘inf’) # 無限大で初期化 best_config = None for config in arima_configs: try: # ARIMA()のパラメータ: # order=(p, d, q): ARIMAの次数 model = ARIMA(train_arima, order=config) model_fit = model.fit() # .aic: 赤池情報量基準 # → 小さいほど良いモデル aic = model_fit.aic print(f”ARIMA{config}: AIC = {aic:.2f}”) # より良いモデルが見つかったら更新 if aic < best_aic: best_aic = aic best_config = config except: print(f"ARIMA{config}: エラー(スキップ)") print() print(f"【最適モデル】ARIMA{best_config} (AIC: {best_aic:.2f})") print() # ============================================ # 最適モデルで予測 # ============================================ model_arima = ARIMA(train_arima, order=best_config) model_arima_fit = model_arima.fit() # .forecast(): 未来を予測 forecast_arima = model_arima_fit.forecast(steps=len(test_arima)) # 精度評価 mae_arima, rmse_arima, mape_arima = calculate_metrics( test_arima, forecast_arima ) print("【ARIMA予測精度】") print(f"MAE: {mae_arima:.2f}") print(f"RMSE: {rmse_arima:.2f}") print(f"MAPE: {mape_arima:.2f}%") print() # ============================================ # 可視化 # ============================================ plt.figure(figsize=(14, 6)) plt.plot(train['date'], train_arima, label='訓練データ', color='blue') plt.plot(test['date'], test_arima, label='実績', color='green', marker='o') plt.plot(test['date'], forecast_arima, label='ARIMA予測', color='red', linestyle='--', marker='s') plt.xlabel('日付', fontsize=12) plt.ylabel('売上(百万円)', fontsize=12) plt.title(f'ARIMAモデルによる予測\nARIMA{best_config}, MAPE: {mape_arima:.1f}%', fontsize=14, fontweight='bold') plt.legend(fontsize=11) plt.grid(alpha=0.3) plt.tight_layout() plt.show()

# ============================================ # Prophetライブラリによる予測 # ============================================ # Prophet = Facebookが開発した時系列予測ライブラリ # → シンプルで使いやすい # → 休日の影響も考慮可能 # → 信頼区間も自動で計算 # インストール(初回のみ) # !pip install prophet –break-system-packages from prophet import Prophet # ============================================ # データ準備(Prophetの形式に変換) # ============================================ # Prophetは特定の列名を要求: # ds: 日付列 # y: 予測対象の値 df_prophet = df[[‘date’, ‘sales’]].copy() df_prophet.columns = [‘ds’, ‘y’] # 訓練とテストに分割 train_prophet = df_prophet[:train_size] test_prophet = df_prophet[train_size:] # ============================================ # Prophetモデルの作成 # ============================================ # Prophet()のパラメータ: # yearly_seasonality: 年次季節性(True/False) # weekly_seasonality: 週次季節性(月次データはFalse) # daily_seasonality: 日次季節性(月次データはFalse) # seasonality_mode: ‘additive'(加法) or ‘multiplicative'(乗法) # interval_width: 信頼区間の幅(0.95 = 95%) model_prophet = Prophet( yearly_seasonality=True, weekly_seasonality=False, daily_seasonality=False, seasonality_mode=’additive’, interval_width=0.95 ) # .fit(): モデルを学習 model_prophet.fit(train_prophet) # ============================================ # 未来のデータフレーム作成 # ============================================ # make_future_dataframe(): 予測用の日付を生成 # periods: 予測する期間数 # freq: 頻度(‘MS’ = 月初) future = model_prophet.make_future_dataframe( periods=len(test_prophet), freq=’MS’ ) # .predict(): 予測を実行 # → 予測値(yhat)、下限(yhat_lower)、上限(yhat_upper)を返す forecast_prophet = model_prophet.predict(future) # テストデータに対する予測値を抽出 # iloc[-len(test_prophet):]: 最後のlen(test_prophet)行 test_forecast = forecast_prophet.iloc[-len(test_prophet):][‘yhat’].values # 精度評価 mae_prophet, rmse_prophet, mape_prophet = calculate_metrics( test_prophet[‘y’].values, test_forecast ) print(“【Prophet予測精度】”) print(f”MAE: {mae_prophet:.2f}”) print(f”RMSE: {rmse_prophet:.2f}”) print(f”MAPE: {mape_prophet:.2f}%”) print() # ============================================ # 未来予測(次の12ヶ月) # ============================================ future_12 = model_prophet.make_future_dataframe(periods=12, freq=’MS’) forecast_12 = model_prophet.predict(future_12) print(“【未来予測(次の12ヶ月)】”) print() # 最後の12行(未来の予測)を表示 for i in range(-12, 0): row = forecast_12.iloc[i] date = row[‘ds’] yhat = row[‘yhat’] # 予測値 lower = row[‘yhat_lower’] # 下限(95%信頼区間) upper = row[‘yhat_upper’] # 上限(95%信頼区間) print(f”{date.strftime(‘%Y年%m月’)}: {yhat:.1f}百万円 ” f”[{lower:.1f}〜{upper:.1f}]”)

# ============================================ # Prophetの可視化 # ============================================ fig, axes = plt.subplots(2, 1, figsize=(14, 10)) # ============================================ # グラフ1: 予測結果 # ============================================ # model_prophet.plot(): Prophet標準の予測グラフ # → 実績(黒点)、予測(青線)、信頼区間(水色帯) ax1 = axes[0] model_prophet.plot(forecast_prophet, ax=ax1) ax1.set_xlabel(‘日付’, fontsize=12) ax1.set_ylabel(‘売上(百万円)’, fontsize=12) ax1.set_title(f’Prophetによる予測(MAPE: {mape_prophet:.1f}%)’, fontsize=14, fontweight=’bold’) ax1.grid(alpha=0.3) # ============================================ # グラフ2: 成分分解 # ============================================ # plot_components(): トレンドと季節性を分解表示 # → トレンド: 長期的な傾向 # → 季節性: 周期的なパターン ax2 = axes[1] model_prophet.plot_components(forecast_prophet) plt.tight_layout() plt.show()

# ============================================ # 全モデルの精度比較 # ============================================ import pandas as pd # 比較表を作成 comparison = pd.DataFrame({ ‘モデル’: [‘移動平均(6ヶ月)’, ‘SES’, ‘Holt’, ‘Holt-Winters’, ‘ARIMA’, ‘Prophet’], ‘MAPE’: [10.5, mape_ses, mape_holt, mape_hw, mape_arima, mape_prophet], ‘MAE’: [8.2, mae_ses, mae_holt, mae_hw, mae_arima, mae_prophet], ‘RMSE’: [10.1, rmse_ses, rmse_holt, rmse_hw, rmse_arima, rmse_prophet] }) # MAPEでソート(精度が高い順) comparison = comparison.sort_values(‘MAPE’) comparison = comparison.reset_index(drop=True) print(“【モデル精度ランキング】”) print(comparison.to_string(index=False)) print() # 最適モデル best_model_name = comparison.iloc[0][‘モデル’] best_mape = comparison.iloc[0][‘MAPE’] print(f”【推奨モデル】{best_model_name} (MAPE: {best_mape:.2f}%)”) print() # ============================================ # 可視化(横棒グラフ) # ============================================ fig, ax = plt.subplots(figsize=(12, 6)) models = comparison[‘モデル’] mapes = comparison[‘MAPE’] # plt.cm.RdYlGn_r: 赤→黄→緑のカラーマップ(反転) # 精度が高い(MAPE小)ほど緑、低いほど赤 colors = plt.cm.RdYlGn_r(np.linspace(0.2, 0.8, len(models))) # barh(): 横棒グラフ bars = ax.barh(models, mapes, color=colors, alpha=0.8, edgecolor=’black’) ax.set_xlabel(‘MAPE (%)’, fontsize=12) ax.set_title(‘予測モデル精度比較’, fontsize=14, fontweight=’bold’) ax.grid(axis=’x’, alpha=0.3) # 値を棒の右側に表示 for bar, mape in zip(bars, mapes): width = bar.get_width() ax.text(width + 0.3, bar.get_y() + bar.get_height()/2, f'{mape:.2f}%’, va=’center’, fontsize=11) plt.tight_layout() plt.show()