数理工房

100本ノック / 確率統計 / 確率・統計Python100本ノック

需要変動を読み解き、生産計画へつなぐ：製造業の時系列分析10本ノック

需要変動を読み解き、生産計画へつなぐ：製造業の時系列分析10本ノック

精密部品工場の日次データを題材に、ACF、PACF、AR、MA、ARIMA、SARIMA、VAR、状態空間モデル、カルマンフィルタ、ProphetをPythonで実装します。予測精度だけを競うのではなく、在庫、生産能力、保全、欠品リスクをどう判断するかまで扱います。

本記事のコードとデータはすべて自己完結しており、外部データは使用しません。乱数seedを固定しているため、同じ環境で結果を再現できます。

[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。

はじめに：この記事で扱う製造業の実務課題

架空の精密部品工場では、翌月の需要に合わせて材料手配、人員配置、設備負荷を決めます。日次需要には増加傾向、曜日性、季節性、連休停止が重なり、単純な平均では欠品と過剰在庫の双方が起きます。また、受注・生産・設備温度も時間差を伴って影響します。

本記事では、時系列の構造を診断し、単変量・多変量・状態空間の予測へ段階的に進みます。各モデルを「当たるか」だけでなく、「どの仮定に基づき、どの意思決定に使えるか」で比較します。

現場でよくある状況

  • 月平均は計画内でも、曜日ごとの欠品と余剰が繰り返される
  • 需要実績、受注残、生産実績、設備データが別々に集計され、時間関係を見ていない
  • 高精度な予測値だけが共有され、予測区間や外れた場合の対応が決まっていない
  • 設備停止や価格改定などの構造変化後も、同じモデルを使い続ける
  • ランダム分割で評価し、未来情報が学習側へ混入している

なぜこの問題は判断が難しいのか

時系列では観測順序が情報です。今日の需要は昨日や先週と独立ではなく、トレンドと周期も含みます。そのため通常のランダム分割は使わず、過去で学習し未来で検証します。

自己相関があっても因果関係とは限りません。さらに、予測誤差が小さくても、欠品損失と在庫費用が非対称なら業務上の最適予測は変わります。モデル選択、予測区間、業務コスト、再学習条件を一体で設計する必要があります。

今回扱うノックの全体像

No.テーマ製造業での問い
081ACF需要は何日前まで似ているか
082PACF中間の影響を除くと重要な遅れはどれか
083ARモデル過去需要だけで短期需要を予測できるか
084MAモデル一時的な需要ショックの影響は何日残るか
085ARIMA非定常な需要水準を差分で扱えるか
086SARIMA曜日性を明示して予測できるか
087VAR受注・需要・生産の相互作用を使えるか
088状態空間モデル見えない基調と季節成分を分離できるか
089カルマンフィルタノイズの多い設備計測から真の状態を逐次推定できるか
090Prophetトレンド・複数季節性・休業日を運用しやすく表現できるか

Python 環境の準備

NumPyで架空データを生成し、pandasで整形、statsmodelsで時系列モデル、Prophetで加法モデル、matplotlibで可視化します。japanize_matplotlib は日本語表示にのみ使います。全ノックで同じ時系列順の訓練・検証分割を使い、RMSEとMAEを比較します。

import contextlib
import io
import logging
import platform
import warnings

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
import scipy
import statsmodels
from IPython.display import display
from sklearn.metrics import mean_absolute_error, root_mean_squared_error
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.ar_model import AutoReg
from statsmodels.tsa.api import VAR
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.statespace.structural import UnobservedComponents
with contextlib.redirect_stdout(io.StringIO()), contextlib.redirect_stderr(io.StringIO()):
    from prophet import Prophet

warnings.filterwarnings("ignore")
logging.getLogger("cmdstanpy").setLevel(logging.WARNING)
logging.getLogger("cmdstanpy").disabled = True
logging.getLogger("prophet").setLevel(logging.WARNING)
pd.set_option("display.max_columns", 20)
pd.set_option("display.float_format", "{:.2f}".format)
rng = np.random.default_rng(20260711)

print(f"Python: {platform.python_version()}")
print(f"NumPy: {np.__version__} / pandas: {pd.__version__} / SciPy: {scipy.__version__}")
print(f"statsmodels: {statsmodels.__version__}")
Python: 3.13.1
NumPy: 2.5.1 / pandas: 3.0.3 / SciPy: 1.18.0
statsmodels: 0.14.6

架空データの作成

2024年1月から2年間の日次データを生成します。需要は緩やかな増加、曜日性、年周期、自己回帰的な揺れ、年末年始の休業影響を持ちます。受注は需要に先行し、生産は需要に遅れて追随します。設備温度センサーには観測ノイズを加えます。

実務では欠損、締め時刻、単位、返品、特急注文、停止日を定義し、予測時点で実際に利用できる説明変数だけを使います。

dates = pd.date_range("2024-01-01", periods=730, freq="D")
n = len(dates)
t = np.arange(n)
weekday = dates.dayofweek.to_numpy()

trend = 118 + 0.035 * t
weekly = np.array([-4, 1, 4, 7, 5, -10, -14])[weekday]
annual = 9 * np.sin(2 * np.pi * t / 365.25) + 4 * np.cos(2 * np.pi * t / 365.25)
shutdown = np.asarray(((dates.month == 1) & (dates.day <= 4)) | ((dates.month == 12) & (dates.day >= 29)))

shock = np.zeros(n)
innovation = rng.normal(0, 5.0, n)
for i in range(2, n):
    shock[i] = 0.58 * shock[i - 1] - 0.18 * shock[i - 2] + innovation[i] + 0.25 * innovation[i - 1]

demand = trend + weekly + annual + shock - 38 * shutdown
orders = np.roll(demand, -2) + rng.normal(5, 6, n)
orders[-2:] = demand[-2:] + rng.normal(5, 6, 2)
production = 0.62 * np.roll(demand, 1) + 0.38 * demand + rng.normal(2, 4, n)
production[0] = demand[0] + rng.normal(2, 4)

true_temp = 61 + 0.012 * t + 1.8 * np.sin(2 * np.pi * t / 30)
sensor_temp = true_temp + rng.normal(0, 2.2, n)

df = pd.DataFrame({
    "需要個数": demand.round(1),
    "受注個数": orders.round(1),
    "生産個数": production.round(1),
    "真の設備温度": true_temp,
    "センサー温度": sensor_temp,
    "休業日": shutdown.astype(int),
}, index=dates)

test_days = 60
train = df.iloc[:-test_days].copy()
test = df.iloc[-test_days:].copy()

display(df.head())
print(f"期間: {df.index.min().date()}〜{df.index.max().date()} / 訓練: {len(train)}日 / 検証: {len(test)}日")

fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(df.index, df["需要個数"], lw=1, label="日次需要")
ax.axvline(test.index[0], color="tab:red", ls="--", label="検証開始")
ax.set_title("架空工場の日次需要（時系列順に訓練・検証を分割）")
ax.set_xlabel("日付")
ax.set_ylabel("需要個数 / 日")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
需要個数 受注個数 生産個数 真の設備温度 センサー温度 休業日
2024-01-01 80.00 91.60 83.30 61.00 61.93 1
2024-01-02 85.20 93.20 80.40 61.39 62.61 1
2024-01-03 90.00 128.50 91.60 61.76 60.63 1
2024-01-04 89.50 127.10 100.10 62.09 61.13 1
2024-01-05 118.30 121.10 99.70 62.39 65.17 0 
期間: 2024-01-01〜2025-12-30 / 訓練: 670日 / 検証: 60日

png

No.081：ACF — 需要に残る時間的な記憶を確認する

実務での意味

自己相関関数（ACF）は、需要と kk 日前の需要がどの程度似るかを測ります。ラグ7、14、21日で相関が残れば曜日性を、緩やかに減衰すればトレンドやAR構造を疑います。発注リードタイム内に相関が残るかは、短期予測を在庫補充へ使えるかの判断材料です。

分析・モデル化の考え方

ラグ kk の標本自己相関は概ね

rk=∑t=k+1T(yt−yˉ)(yt−k−yˉ)∑t=1T(yt−yˉ)2r_k=\frac{\sum_{t=k+1}^{T}(y_t-\bar y)(y_{t-k}-\bar y)}{\sum_{t=1}^{T}(y_t-\bar y)^2}

です。ただしトレンドや季節性がある系列のACFは見かけ上高くなります。原系列と前年差分・季節差分を比較し、定常化後の構造も確認します。

Pythonで確認する

acf_lags = [1, 2, 7, 14, 21, 28]
acf_table = pd.DataFrame({
    "ラグ（日）": acf_lags,
    "原系列ACF": [train["需要個数"].autocorr(lag=k) for k in acf_lags],
    "7日差分ACF": [train["需要個数"].diff(7).dropna().autocorr(lag=k) for k in acf_lags],
})
display(acf_table)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(train["需要個数"], lags=35, ax=axes[0], zero=False)
axes[0].set_title("日次需要のACF")
axes[0].set_xlabel("ラグ（日）")
axes[0].set_ylabel("自己相関")
axes[0].grid(alpha=0.3)
plot_acf(train["需要個数"].diff(7).dropna(), lags=35, ax=axes[1], zero=False)
axes[1].set_title("7日差分後のACF")
axes[1].set_xlabel("ラグ（日）")
axes[1].set_ylabel("自己相関")
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
ラグ（日） 原系列ACF 7日差分ACF
0 1 0.75 0.69
1 2 0.38 0.36
2 7 0.61 -0.49
3 14 0.61 -0.00
4 21 0.62 0.06
5 28 0.58 -0.05

png

結果の読み取り

原系列では短期の相関と7日周期が見え、昨日と同じ水準を置く単純予測にも一定の根拠があります。一方、7日差分後は多くの相関が弱まり、元の高い相関の一部が曜日性とトレンドに由来すると分かります。ACFだけで次数を確定せず、残差ACFと未来期間での誤差まで確認します。

No.082：PACF — 直接効いているラグを絞り込む

実務での意味

PACF（偏自己相関）は、間にある日の影響を統制したうえで、あるラグが需要へ直接持つ関係を見ます。ラグ1が強ければ直前日の情報、ラグ7が残れば同じ曜日の情報を生産計画へ入れる候補になります。

分析・モデル化の考え方

ラグ kk のPACFは、yty_t を yt−1,…,yt−ky_{t-1},\ldots,y_{t-k} で回帰したときの yt−ky_{t-k} の係数に対応します。AR(pp)では理論上PACFが pp より後で打ち切られるため、AR次数の候補選定に使えます。ただし季節性や非定常性が残ると単純な読み方はできません。

Pythonで確認する

stationary_demand = train["需要個数"].diff(7).dropna()
fig, ax = plt.subplots(figsize=(10, 4))
plot_pacf(stationary_demand, lags=28, method="ywm", ax=ax, zero=False)
ax.set_title("7日差分した需要のPACF")
ax.set_xlabel("ラグ（日）")
ax.set_ylabel("偏自己相関")
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()

from statsmodels.tsa.stattools import pacf
pacf_values = pacf(stationary_demand, nlags=14, method="ywm")
display(pd.DataFrame({"ラグ（日）": np.arange(1, 15), "PACF": pacf_values[1:]}).sort_values("PACF", key=abs, ascending=False).head(7))

png

ラグ（日） PACF
0 1 0.69
7 8 0.31
6 7 -0.26
5 6 -0.23
1 2 -0.23
4 5 -0.16
8 9 -0.15

結果の読み取り

大きな偏自己相関を持つラグはARモデルの候補ですが、「統計的に有意」と「生産計画上重要」は同じではありません。候補次数を少数に絞り、AICなどの情報量規準と時系列検証で比較します。ラグを増やしすぎると説明が難しくなり、構造変化への追随も遅くなります。

No.083：ARモデル — 過去需要から短期需要を予測する

実務での意味

自己回帰（AR）モデルは、直近の需要実績から次の需要を予測します。短い調達リードタイムや日々の要員微調整に向きます。一方、価格改定や顧客の新規採用など、過去にない変化は自力で予測できません。

分析・モデル化の考え方

AR(pp)は

yt=c+ϕ1yt−1+⋯+ϕpyt−p+εty_t=c+\phi_1y_{t-1}+\cdots+\phi_py_{t-p}+\varepsilon_t

と表します。ここでは曜日周期を含めるためラグ14を使います。検証期間を末尾60日に固定し、未来値を学習へ混ぜずにRMSE・MAEを算出します。

Pythonで確認する

ar_fit = AutoReg(train["需要個数"], lags=14, trend="ct", old_names=False).fit()
ar_pred = ar_fit.predict(start=len(train), end=len(df) - 1, dynamic=False)
ar_pred.index = test.index

ar_metrics = pd.DataFrame({
    "モデル": ["AR(14)"],
    "RMSE": [root_mean_squared_error(test["需要個数"], ar_pred)],
    "MAE": [mean_absolute_error(test["需要個数"], ar_pred)],
    "AIC": [ar_fit.aic],
})
display(ar_metrics)

fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(train.index[-90:], train["需要個数"].iloc[-90:], label="訓練実績", color="0.55")
ax.plot(test.index, test["需要個数"], label="検証実績", color="black")
ax.plot(test.index, ar_pred, label="AR(14)予測", color="tab:blue")
ax.set_title("ARモデルによる日次需要予測")
ax.set_xlabel("日付")
ax.set_ylabel("需要個数 / 日")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
モデル RMSE MAE AIC
0 AR(14) 12.14 9.96 4370.51

png

結果の読み取り

ARモデルは水準と短期的な波を追いますが、長い予測ほど同じ再帰構造へ収束しやすくなります。MAEは平均的な必要調整数、RMSEは大外しを強く罰する指標として読みます。能力超過や欠品に重い損失があるなら、対称な誤差指標だけで採用を決めません。

No.084：MAモデル — 一時的ショックの残り方を表現する

実務での意味

移動平均（MA）モデルは、過去の予測誤差が現在へ残る構造を表します。特急注文、短時間停止、納入前倒しのような一時ショックが数日で解消する工程に適します。ここでいうMAは、単純移動平均による平滑化とは別物です。

分析・モデル化の考え方

MA(qq)は

yt=μ+εt+θ1εt−1+⋯+θqεt−qy_t=\mu+\varepsilon_t+\theta_1\varepsilon_{t-1}+\cdots+\theta_q\varepsilon_{t-q}

です。今回はトレンド項を併用したMA(2)を当て、過去2日の未説明ショックを表します。誤差は直接観測できないため最尤法で推定します。

Pythonで確認する

ma_fit = ARIMA(train["需要個数"], order=(0, 0, 2), trend="ct").fit()
ma_forecast = ma_fit.get_forecast(steps=len(test))
ma_pred = pd.Series(ma_forecast.predicted_mean.to_numpy(), index=test.index)
ma_ci = ma_forecast.conf_int(alpha=0.05)

display(pd.DataFrame({
    "モデル": ["MA(2)+trend"],
    "RMSE": [root_mean_squared_error(test["需要個数"], ma_pred)],
    "MAE": [mean_absolute_error(test["需要個数"], ma_pred)],
    "AIC": [ma_fit.aic],
}))

fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(test.index, test["需要個数"], label="検証実績", color="black")
ax.plot(test.index, ma_pred, label="MA(2)予測", color="tab:orange")
ax.fill_between(test.index, ma_ci.iloc[:, 0], ma_ci.iloc[:, 1], color="tab:orange", alpha=0.2, label="95%予測区間")
ax.set_title("MAモデルによる需要予測と予測区間")
ax.set_xlabel("日付")
ax.set_ylabel("需要個数 / 日")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
モデル RMSE MAE AIC
0 MA(2)+trend 12.64 10.32 4712.80

png

結果の読み取り

MA(2)は短期ショックを表せますが、明示的な曜日性を持たないため週次の上下を捉え切れません。予測区間は平均需要の信頼区間ではなく、将来の個々の観測が入り得る範囲です。安全在庫では上側分位点が重要ですが、誤差分布とリードタイム中の累積需要を別途検討します。

No.085：ARIMA — 差分で変化量をモデル化する

実務での意味

需要水準が増加・減少していると、一定平均を仮定するモデルは古い水準へ引き戻されます。ARIMAは差分によって水準変化を扱い、立上げ期や緩やかな市場成長下の基準予測に使えます。

分析・モデル化の考え方

ARIMA(p,d,qp,d,q)は、dd 回差分した系列にARMA(p,qp,q)を当てます。1階差分は Δyt=yt−yt−1\Delta y_t=y_t-y_{t-1} です。ここではARIMA(2,1,2)を使います。差分次数を上げすぎるとノイズを増幅するため、グラフ、単位根検定、残差、未来精度を併用します。

Pythonで確認する

arima_fit = ARIMA(train["需要個数"], order=(2, 1, 2), trend="t").fit()
arima_fc = arima_fit.get_forecast(steps=len(test))
arima_pred = pd.Series(arima_fc.predicted_mean.to_numpy(), index=test.index)
arima_ci = arima_fc.conf_int()

display(pd.DataFrame({
    "モデル": ["ARIMA(2,1,2)"],
    "RMSE": [root_mean_squared_error(test["需要個数"], arima_pred)],
    "MAE": [mean_absolute_error(test["需要個数"], arima_pred)],
    "AIC": [arima_fit.aic],
}))

fig, axes = plt.subplots(1, 2, figsize=(13, 4))
axes[0].plot(train.index[-120:], train["需要個数"].iloc[-120:])
axes[0].set_title("需要水準（差分前）")
axes[0].set_xlabel("日付")
axes[0].set_ylabel("需要個数 / 日")
axes[0].grid(alpha=0.3)
axes[1].plot(test.index, test["需要個数"], label="実績", color="black")
axes[1].plot(test.index, arima_pred, label="ARIMA予測", color="tab:green")
axes[1].fill_between(test.index, arima_ci.iloc[:, 0], arima_ci.iloc[:, 1], color="tab:green", alpha=0.2)
axes[1].set_title("ARIMA予測")
axes[1].set_xlabel("日付")
axes[1].set_ylabel("需要個数 / 日")
axes[1].grid(alpha=0.3)
axes[1].legend()
plt.tight_layout()
plt.show()
モデル RMSE MAE AIC
0 ARIMA(2,1,2) 12.71 10.39 4658.08

png

結果の読み取り

ARIMAは局所的な水準変化へ追随しますが、週次パターンを明示していません。AICは同じデータに当てた候補間の相対比較に使い、異なる検証区間の業務損失とは分けて扱います。残差に7日周期が残るなら、次のSARIMAが自然な候補です。

No.086：SARIMA — 曜日性を含めて生産量を見通す

実務での意味

曜日ごとに受注・出荷パターンが違う工場では、季節性を無視すると毎週同じ日に欠品します。SARIMAは通常のARIMAに季節差分と季節AR・MAを加え、週次・月次など一定周期の繰返しを表現します。

分析・モデル化の考え方

SARIMAは (p,d,q)×(P,D,Q)s(p,d,q)\times(P,D,Q)_s と書きます。s=7s=7 は日次データの週周期です。今回は (1,1,1)×(1,0,1)7(1,1,1)\times(1,0,1)_7 とし、複雑さを抑えます。休日が移動する場合や設備停止が不規則な場合は、外生変数を追加します。

Pythonで確認する

sarima_fit = SARIMAX(
    train["需要個数"],
    order=(1, 1, 1),
    seasonal_order=(1, 0, 1, 7),
    trend="t",
    enforce_stationarity=False,
    enforce_invertibility=False,
).fit(disp=False)
sarima_fc = sarima_fit.get_forecast(steps=len(test))
sarima_pred = pd.Series(sarima_fc.predicted_mean.to_numpy(), index=test.index)
sarima_ci = sarima_fc.conf_int()

model_metrics = pd.DataFrame({
    "モデル": ["AR(14)", "MA(2)+trend", "ARIMA(2,1,2)", "SARIMA-weekly"],
    "RMSE": [
        root_mean_squared_error(test["需要個数"], ar_pred),
        root_mean_squared_error(test["需要個数"], ma_pred),
        root_mean_squared_error(test["需要個数"], arima_pred),
        root_mean_squared_error(test["需要個数"], sarima_pred),
    ],
    "MAE": [
        mean_absolute_error(test["需要個数"], ar_pred),
        mean_absolute_error(test["需要個数"], ma_pred),
        mean_absolute_error(test["需要個数"], arima_pred),
        mean_absolute_error(test["需要個数"], sarima_pred),
    ],
}).sort_values("RMSE")
display(model_metrics)

fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(test.index, test["需要個数"], label="検証実績", color="black", lw=1.5)
ax.plot(test.index, sarima_pred, label="SARIMA予測", color="tab:red")
ax.fill_between(test.index, sarima_ci.iloc[:, 0], sarima_ci.iloc[:, 1], color="tab:red", alpha=0.18, label="95%予測区間")
ax.set_title("週次季節性を含むSARIMA需要予測")
ax.set_xlabel("日付")
ax.set_ylabel("需要個数 / 日")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
モデル RMSE MAE
0 AR(14) 12.14 9.96
1 MA(2)+trend 12.64 10.32
2 ARIMA(2,1,2) 12.71 10.39
3 SARIMA-weekly 17.13 15.12

png

結果の読み取り

今回のSARIMAは、曜日性を明示したにもかかわらず単純なARよりRMSE・MAEが悪化しました。季節項を入れれば必ず改善するわけではなく、トレンド仕様や検証末尾の年末休業を説明できていないことが誤差要因です。1回の結果で棄却せず、次数と休日外生変数をローリング検証します。予測区間上限が生産能力を超える日には、前倒し生産や外注のルールを結び付けます。

No.087：VAR — 需要・受注・生産の相互作用をモデル化する

実務での意味

需要単独ではなく、先行する受注と追随する生産を同時に扱うと、需給の時間差を把握できます。VARは複数系列を互いの過去値で説明し、販売・生産管理間で共通のシナリオを作る基礎になります。

分析・モデル化の考え方

KK 変量のVAR(pp)は

yt=c+A1yt−1+⋯+Apyt−p+εt\mathbf y_t=\mathbf c+A_1\mathbf y_{t-1}+\cdots+A_p\mathbf y_{t-p}+\boldsymbol\varepsilon_t

です。非定常な水準のまま偽回帰を起こさないよう、ここでは1階差分を使います。系列数とラグ数を増やすほどパラメータが急増するため、データ量と業務解釈を優先します。

Pythonで確認する

var_cols = ["需要個数", "受注個数", "生産個数"]
var_train_diff = train[var_cols].diff().dropna()
var_fit = VAR(var_train_diff).fit(maxlags=7, ic="aic")
horizon = 14
diff_fc = var_fit.forecast(var_train_diff.to_numpy()[-var_fit.k_ar:], steps=horizon)
var_levels = train[var_cols].iloc[-1].to_numpy() + np.cumsum(diff_fc, axis=0)
var_forecast = pd.DataFrame(var_levels, index=test.index[:horizon], columns=var_cols)

print(f"AICで選ばれたラグ次数: {var_fit.k_ar}")
display(var_forecast.head())

fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)
for ax, col in zip(axes, var_cols):
    ax.plot(test.index[:horizon], test[col].iloc[:horizon], label="実績", color="black")
    ax.plot(var_forecast.index, var_forecast[col], label="VAR予測", color="tab:purple")
    ax.set_title(f"VARによる{col}の14日予測")
    ax.set_xlabel("日付")
    ax.set_ylabel("個数 / 日")
    ax.grid(alpha=0.3)
    ax.legend()
plt.tight_layout()
plt.show()
AICで選ばれたラグ次数: 7
需要個数 受注個数 生産個数
2025-11-01 124.63 145.22 133.53
2025-11-02 129.16 147.47 129.57
2025-11-03 137.48 145.05 136.03
2025-11-04 139.19 146.75 141.50
2025-11-05 139.17 143.29 142.63

png

結果の読み取り

3系列を同じ時間軸で予測すると、需要増に生産が追随できるかを確認できます。ただしVARの係数は因果効果ではありません。受注が将来需要を先取りして見えるのは生成上の時間関係であり、本番では入力確定時刻とデータ改訂を監査します。長期予測より、短期の需給ギャップ検知に向く構成です。

No.088：状態空間モデル — 観測の背後にある基調を推定する

実務での意味

日々の需要はノイズを含みますが、生産能力や購買契約は基調に合わせて決めたいものです。状態空間モデルは、観測できない需要レベル・傾き・季節成分を状態として分け、欠測や不規則な変化にも柔軟に対応します。

分析・モデル化の考え方

一般形は観測方程式と状態方程式

yt=Ztαt+εt,αt+1=Ttαt+ηty_t=Z_t\alpha_t+\varepsilon_t,\qquad \alpha_{t+1}=T_t\alpha_t+\eta_t

です。αt\alpha_t が潜在状態、εt\varepsilon_t が観測誤差、ηt\eta_t が状態変化です。ここでは局所線形トレンドと週次季節成分を推定します。

Pythonで確認する

ucm = UnobservedComponents(
    train["需要個数"],
    level="local linear trend",
    seasonal=7,
    autoregressive=2,
)
ucm_fit = ucm.fit(disp=False)
ucm_fc = ucm_fit.get_forecast(steps=len(test))
ucm_pred = pd.Series(ucm_fc.predicted_mean.to_numpy(), index=test.index)
ucm_ci = ucm_fc.conf_int()

display(pd.DataFrame({
    "モデル": ["状態空間（局所線形+週次+AR2）"],
    "RMSE": [root_mean_squared_error(test["需要個数"], ucm_pred)],
    "MAE": [mean_absolute_error(test["需要個数"], ucm_pred)],
    "AIC": [ucm_fit.aic],
}))

fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(test.index, test["需要個数"], label="検証実績", color="black")
ax.plot(test.index, ucm_pred, label="状態空間モデル予測", color="tab:brown")
ax.fill_between(test.index, ucm_ci.iloc[:, 0], ucm_ci.iloc[:, 1], color="tab:brown", alpha=0.2)
ax.set_title("状態空間モデルによる需要予測")
ax.set_xlabel("日付")
ax.set_ylabel("需要個数 / 日")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
モデル RMSE MAE AIC
0 状態空間（局所線形+週次+AR2） 11.26 8.90 4364.45

png

結果の読み取り

状態空間モデルは予測に加え、「基調がどの程度動くか」を状態分散として表現します。変化を許しすぎればノイズへ過適合し、許さなければ構造変化へ追随しません。モデル部品を業務知識と対応させ、状態分散、残差、予測区間の被覆率を監視することが導入の要点です。

No.089：カルマンフィルタ — 設備状態を逐次推定する

実務での意味

温度・振動・圧力センサーはノイズを含みます。閾値を観測値へ直接適用すると誤報が増えます。カルマンフィルタは新しい測定が届くたびに状態推定を更新し、予防保全や異常監視の安定した入力を作ります。

分析・モデル化の考え方

予測と更新を繰り返します。1次元局所レベルモデルでは、予測分散 Pt∣t−1=Pt−1∣t−1+QP_{t|t-1}=P_{t-1|t-1}+Q、カルマンゲイン

Kt=Pt∣t−1Pt∣t−1+RK_t=\frac{P_{t|t-1}}{P_{t|t-1}+R}

更新値 x^t∣t=x^t∣t−1+Kt(yt−x^t∣t−1)\hat x_{t|t}=\hat x_{t|t-1}+K_t(y_t-\hat x_{t|t-1}) です。QQ は真の状態変化、RR は観測ノイズの分散です。

Pythonで確認する

observations = df["センサー温度"].to_numpy()
q, r = 0.08, 2.2**2
x_est, p = observations[0], 10.0
filtered, gains = [], []

for y in observations:
    x_pred = x_est
    p_pred = p + q
    k = p_pred / (p_pred + r)
    x_est = x_pred + k * (y - x_pred)
    p = (1 - k) * p_pred
    filtered.append(x_est)
    gains.append(k)

df["カルマン推定温度"] = filtered
sensor_rmse = root_mean_squared_error(df["真の設備温度"], df["センサー温度"])
kalman_rmse = root_mean_squared_error(df["真の設備温度"], df["カルマン推定温度"])
display(pd.DataFrame({
    "系列": ["生センサー", "カルマン推定"],
    "真の状態に対するRMSE": [sensor_rmse, kalman_rmse],
    "直近カルマンゲイン": [np.nan, gains[-1]],
}))

fig, ax = plt.subplots(figsize=(12, 4))
window = df.iloc[-120:]
ax.plot(window.index, window["センサー温度"], color="0.7", lw=1, label="観測値")
ax.plot(window.index, window["真の設備温度"], color="black", lw=2, label="真の状態（検証用）")
ax.plot(window.index, window["カルマン推定温度"], color="tab:cyan", lw=2, label="カルマン推定")
ax.set_title("カルマンフィルタによる設備温度の逐次推定")
ax.set_xlabel("日付")
ax.set_ylabel("設備温度（℃）")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
系列 真の状態に対するRMSE 直近カルマンゲイン
0 生センサー 2.23 NaN
1 カルマン推定 1.17 0.12

png

結果の読み取り

生センサーより推定系列のRMSEが小さくなり、短期ノイズを抑えています。ただし急激な真の変化にも遅れるため、平滑化後の値だけで安全停止を判断してはいけません。QQ と RR は保全履歴、校正試験、擬似異常で調整し、生値のハード閾値と状態推定のソフト警報を役割分担します。

No.090：Prophet — トレンド・季節性・休業日を説明可能な部品に分ける

実務での意味

Prophetはトレンド、週次・年次季節性、休日効果を加法的に表し、カレンダー要因を説明しやすいモデルです。営業・生産管理と「何曜日に増えるか」「休業でどれだけ下がるか」を共有しやすく、定期的な再学習の基準モデルに向きます。

分析・モデル化の考え方

基本形は

y(t)=g(t)+s(t)+h(t)+εty(t)=g(t)+s(t)+h(t)+\varepsilon_t

で、g(t)g(t) は区分的トレンド、s(t)s(t) はFourier級数による季節性、h(t)h(t) は休日効果です。柔軟性を上げれば過適合しやすいため、changepoint・季節性・休日窓を業務変化と照合します。

Pythonで確認する

prophet_train = train.reset_index().rename(columns={"index": "ds", "需要個数": "y"})[["ds", "y"]]
holiday_dates = df.index[df["休業日"].eq(1)]
holidays = pd.DataFrame({"holiday": "工場休業", "ds": holiday_dates})

prophet_model = Prophet(
    holidays=holidays,
    weekly_seasonality=True,
    yearly_seasonality=True,
    daily_seasonality=False,
    interval_width=0.95,
    changepoint_prior_scale=0.05,
)
prophet_model.fit(prophet_train)
future = pd.DataFrame({"ds": test.index})
prophet_fc = prophet_model.predict(future).set_index("ds")

display(pd.DataFrame({
    "モデル": ["Prophet"],
    "RMSE": [root_mean_squared_error(test["需要個数"], prophet_fc["yhat"])],
    "MAE": [mean_absolute_error(test["需要個数"], prophet_fc["yhat"])],
    "区間内割合": [((test["需要個数"] >= prophet_fc["yhat_lower"]) & (test["需要個数"] <= prophet_fc["yhat_upper"])).mean()],
}))

fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(test.index, test["需要個数"], color="black", label="検証実績")
ax.plot(prophet_fc.index, prophet_fc["yhat"], color="tab:pink", label="Prophet予測")
ax.fill_between(prophet_fc.index, prophet_fc["yhat_lower"], prophet_fc["yhat_upper"], color="tab:pink", alpha=0.2, label="95%予測区間")
ax.set_title("Prophetによる需要予測（週次・年次・工場休業）")
ax.set_xlabel("日付")
ax.set_ylabel("需要個数 / 日")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()
モデル RMSE MAE 区間内割合
0 Prophet 7.55 6.05 0.88

png

結果の読み取り

Prophetは休業による需要低下と季節性を分けて表現し、予測区間の被覆率も確認できます。区間内割合が名目95%から大きく外れるなら、不確かさを過小・過大評価しています。便利な自動化があっても、将来の休業日登録、受注制度の変更、予測時点で未知の情報を混ぜない管理は利用側の責任です。

対象ノックを通して見える実務上の示唆

  1. 診断から始める：ACF・PACFで時間構造を見てから次数候補を作り、モデル名を先に決めない。
  2. 未来で評価する：時系列順のホールドアウトとローリング検証を使い、ランダム分割によるリーケージを防ぐ。
  3. 点予測を業務判断へ翻訳する：予測区間上限と能力、分位点と安全在庫、予測誤差と欠品損失を対応させる。
  4. 単変量と多変量を使い分ける：受注などの先行情報は有用だが、予測時点で確定しているかを必ず確認する。
  5. 状態推定と安全制御を分ける：カルマンフィルタはノイズ低減に有効だが、急変を遅れて捉えるリスクを持つ。
  6. 単一の勝者を固定しない：季節性、構造変化、予測期間に応じて、単純基準モデルを含む候補を継続比較する。

実務導入する場合に必要なこと

  • 意思決定の定義：予測対象、粒度、リードタイム、締め時刻、更新頻度、利用者を明文化する
  • データ契約：欠損、取消、返品、休業、特急注文、単位変更、マスタ改訂を記録する
  • 検証設計：ローリング検証で複数季節を評価し、RMSEだけでなく欠品・在庫・残業コストを測る
  • 不確かさの運用：予測区間ごとの前倒し生産、外注、追加発注、管理者承認ルールを決める
  • 監視：予測誤差、区間被覆率、残差自己相関、入力分布、欠測率、計算失敗を記録する
  • 変更管理：設備増強、価格改定、大口顧客、連休など構造変化をモデルへ反映し、再学習履歴を残す
  • 責任分界：モデルが提案し、現場が承認する範囲と、安全・品質上のハード制約を明確にする

まとめ

No.081〜No.090では、ACF・PACFによる診断から、AR・MA・ARIMA・SARIMA、VAR、状態空間モデル、カルマンフィルタ、Prophetまでを同じ架空工場データで実装しました。

時系列分析の価値は、最も複雑なモデルを選ぶことではありません。時間構造と不確かさを可視化し、在庫・能力・保全の判断ルールへ接続し、外れた後に改善できる仕組みを作ることです。まず単純な基準予測と時系列検証を整え、そのうえで業務上説明できる複雑さを追加してください。

法人向けのご相談

数理工房では、製造業の需要予測、生産・在庫計画、設備状態推定、時系列データ基盤、Python研修、PoCから運用定着までをご支援しています。「予測はあるが発注判断につながらない」「モデルの精度監視を設計したい」「現場データの時刻・欠損・マスタに課題がある」といった段階から、業務要件とデータの双方を整理します。

📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。

Related Posts

関連する記事