【ベンチマーク】間欠需要データ

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
n = 72  # 6年分の月次
dates = pd.date_range("2018-01-01", periods=n, freq="MS")

# 約70%の確率でゼロ、30%の確率で需要発生
demand_occurs = rng.random(n) < 0.30
demand_size = rng.lognormal(mean=2.5, sigma=0.8, size=n)
values = np.where(demand_occurs, demand_size, 0.0)
series = pd.Series(values, index=dates)

h = 12
train_raw = series.iloc[:-h]
test_raw = series.iloc[-h:]

fig, axes = plt.subplots(1, 2, figsize=(13, 4))

ax = axes[0]
ax.bar(series.index, series.values, color="#2563eb", width=20)
ax.axvline(series.index[-h], color="#ef4444", linestyle="--", label="テスト開始")
ax.set_title("間欠需要データ（赤破線より右がテスト期間）")
ax.set_ylabel("需要量")
ax.legend()
ax.grid(alpha=0.2)

ax = axes[1]
zero_pct = (series == 0).sum() / len(series) * 100
nonzero = series[series > 0]
ax.hist(nonzero, bins=15, color="#10b981", edgecolor="white", alpha=0.8)
ax.set_title(f"非ゼロ時の需要分布（ゼロ率 {zero_pct:.0f}%）")
ax.set_xlabel("需要量")
ax.set_ylabel("頻度")
ax.grid(alpha=0.2)

fig.tight_layout()
plt.show()

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108

import warnings
from statsmodels.tsa.holtwinters import ExponentialSmoothing, SimpleExpSmoothing
from statsmodels.tsa.arima.model import ARIMA

warnings.filterwarnings("ignore")

def compute_mase(actual, predicted, train, m=12):
    e = np.abs(actual - predicted)
    scale = np.abs(train.values[m:] - train.values[:-m]).mean()
    return e.mean() / scale if scale > 0 else np.nan

def croston_forecast(train, h, alpha=0.1):
    """Croston法: 需要間隔と需要サイズを分離してSES。"""
    vals = train.values
    # 需要発生時点を抽出
    demand_idx = np.where(vals > 0)[0]
    if len(demand_idx) < 2:
        mean_demand = vals[vals > 0].mean() if (vals > 0).any() else 0
        return np.full(h, mean_demand)

    # 需要サイズ
    sizes = vals[demand_idx]
    # 需要間隔（発生間の期間数）
    intervals = np.diff(demand_idx).astype(float)

    # SESで平滑化
    z = sizes[1:]  # 最初のサイズは間隔と対応しないので除く
    p = intervals

    z_hat = z[0]
    p_hat = p[0]
    for j in range(1, len(z)):
        z_hat = alpha * z[j] + (1 - alpha) * z_hat
        p_hat = alpha * p[j] + (1 - alpha) * p_hat

    # 予測: 需要サイズ / 需要間隔 = 期待需要/期
    forecast_value = z_hat / p_hat if p_hat > 0 else 0
    return np.full(h, forecast_value)

results = {}

# 季節ナイーブ
last12 = train_raw.values[-12:]
results["季節ナイーブ"] = np.tile(last12, int(np.ceil(h / 12)))[:h]

# SES
try:
    results["SES"] = SimpleExpSmoothing(
        train_raw, initialization_method="estimated"
    ).fit().forecast(h).values
except Exception:
    results["SES"] = np.full(h, np.nan)

# HW加法
try:
    results["HW加法"] = ExponentialSmoothing(
        train_raw, trend="add", seasonal=None,
        initialization_method="estimated",
    ).fit().forecast(h).values
except Exception:
    results["HW加法"] = np.full(h, np.nan)

# ARIMA(1,0,0)
try:
    results["ARIMA"] = ARIMA(train_raw, order=(1, 0, 0)).fit().forecast(steps=h).values
except Exception:
    results["ARIMA"] = np.full(h, np.nan)

# Croston法
results["Croston"] = croston_forecast(train_raw, h, alpha=0.1)

# MASE + RMSE
metrics = {}
for name, pred in results.items():
    mase = compute_mase(test_raw.values, pred, train_raw, m=12)
    rmse = np.sqrt(np.mean((test_raw.values - pred) ** 2))
    metrics[name] = {"MASE": mase, "RMSE": rmse}

# 棒グラフ
names = list(metrics.keys())
mase_vals = [metrics[n]["MASE"] for n in names]
rmse_vals = [metrics[n]["RMSE"] for n in names]

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

colors = ["#10b981" if v < 1 else "#ef4444" for v in mase_vals]
bars = axes[0].barh(names, mase_vals, color=colors, edgecolor="white")
axes[0].axvline(1.0, color="#1e40af", linestyle="--", linewidth=1)
axes[0].set_xlabel("MASE")
axes[0].set_title("MASE 比較（低いほど良い）")
for bar, v in zip(bars, mase_vals):
    if not np.isnan(v):
        axes[0].text(bar.get_width() + 0.03, bar.get_y() + bar.get_height() / 2,
                     f"{v:.2f}", va="center", fontsize=9)
axes[0].grid(alpha=0.2, axis="x")

axes[1].barh(names, rmse_vals, color="#2563eb", edgecolor="white")
axes[1].set_xlabel("RMSE")
axes[1].set_title("RMSE 比較（低いほど良い）")
for bar, v in zip(axes[1].patches, rmse_vals):
    if not np.isnan(v):
        axes[1].text(bar.get_width() + 0.1, bar.get_y() + bar.get_height() / 2,
                     f"{v:.1f}", va="center", fontsize=9)
axes[1].grid(alpha=0.2, axis="x")

fig.suptitle("間欠需要データ: モデル別精度比較", fontsize=13)
fig.tight_layout()
plt.show()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

fig, axes = plt.subplots(2, 2, figsize=(13, 8))
test_idx = test_raw.index

# 実測 vs 予測（2モデル）
for i, (name, color) in enumerate([("季節ナイーブ", "#2563eb"), ("SES", "#ef4444")]):
    ax = axes[0, i]
    width = 8
    ax.bar(test_idx - pd.Timedelta(days=width), test_raw.values,
           width=width, color="#94a3b8", label="実測")
    ax.bar(test_idx + pd.Timedelta(days=width), results[name],
           width=width, color=color, label=name)
    ax.set_title(f"{name} vs 実測")
    ax.legend(fontsize=8)
    ax.grid(alpha=0.2)

# Croston vs ARIMA
for i, (name, color) in enumerate([("Croston", "#10b981"), ("ARIMA", "#f97316")]):
    ax = axes[1, i]
    width = 8
    ax.bar(test_idx - pd.Timedelta(days=width), test_raw.values,
           width=width, color="#94a3b8", label="実測")
    ax.bar(test_idx + pd.Timedelta(days=width), results[name],
           width=width, color=color, label=name)
    ax.set_title(f"{name} vs 実測")
    ax.legend(fontsize=8)
    ax.grid(alpha=0.2)

fig.suptitle("間欠需要: 各モデルの予測パターン", fontsize=13)
fig.tight_layout()
plt.show()