Prophet（分解＋ベイズ）

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：ベイズ構造時系列・STL分解　|　関連：季節ARIMA（SARIMA）（自己相関で予測）・予測の評価指標と時系列CV（ホールドアウト）

⚙️ 要最新確認：本ノートは prophet 1.3.0（2026-01-27 リリース、cmdstanpy バックエンド・各プラットフォーム向けプリコンパイル済 Stan モデル同梱）で検証。Prophet は API・既定値が版で変わり得るので、実行時に import prophet; print(prophet.__version__) で確認してください。既定の不確実区間は interval_width=0.8（80%）、既定は MAP 推定（mcmc_samples=0）です。

要点（BLUF）

Prophet＝区分線形トレンド（changepoints）＋フーリエ級数の季節＋休日効果を加法分解し、ベイズ（Stan）で当てる実務向け予測ツール。成分が解釈しやすく、欠測・外れ値に強いのが売り。
ベイズ構造時系列と同じ「分解＋ベイズ」の発想を、自動化・頑健化して使いやすくしたもの。make_future_dataframe → predict で **yhat と yhat_lower/yhat_upper（不確実区間）**つき予測が一行で出ます。
トレンド屈折（changepoint）を捉えるのが強み。一方で、純粋な自己相関構造は季節ARIMA（SARIMA）の方が得意——分解志向の Prophet と自己回帰の ARIMA は使い分けます。擬似系列で成分復元・屈折検出・予測区間・SARIMA 比較を確かめます。

1. Prophet のモデル：3成分の加法分解

Prophet は時系列を、時刻 $t$ の回帰として組み立てます（自己回帰でなく曲線当て）：

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

$g(t)$ … トレンド（区分線形 or ロジスティック）
$s(t)$ … 季節（フーリエ級数、週・年など複数可）
$h(t)$ … 休日・イベント効果
$\varepsilon_t$ … 観測ノイズ

**トレンド $g(t)$ （区分線形）**が Prophet の特徴です。基準勾配 $k$ に対し、あらかじめ置いた候補点（changepoints） $s_j$ で勾配が $\delta_j$ だけ変わってよいとします：

g(t) = \Big(k + \textstyle\sum_{j:\,s_j<t}\delta_j\Big)\,t + \big(m + \textstyle\sum_{j:\,s_j<t}\gamma_j\big)

ここで勾配調整 $\delta_j$ にラプラス事前 $\delta_j\sim\text{Laplace}(0,\tau)$ を置きます（ $\tau=$ changepoint_prior_scale）。ラプラス事前はスパースを促す（正則化との関係の L1 と同型）ので、本当に必要な少数の点でだけ勾配が折れ、残りは 0 に潰れる——これが「自動 changepoint 検出」の正体です。 $\tau$ を上げるとトレンドは柔軟（折れやすい）、下げると滑らか。

季節 $s(t)$ は周期 $P$ のフーリエ級数（ベイズ構造時系列と同じ）：

s(t) = \sum_{k=1}^{K}\Big[a_k\cos\frac{2\pi k t}{P} + b_k\sin\frac{2\pi k t}{P}\Big]

係数 $a_k,b_k$ に正規事前を置いてベイズ回帰で当てます。年次・週次を別々に足せるのが実務で便利。

flowchart LR
    Y["観測 y(t)（欠測・外れ値があってよい）"] --> FIT["Stan で MAP/MCMC 当てはめ"]
    CP["区分線形トレンド g(t)：changepoint にラプラス事前"] --> FIT
    SE["フーリエ季節 s(t)：年次・週次"] --> FIT
    HO["休日・イベント h(t)"] --> FIT
    FIT --> PR["predict：yhat + yhat_lower/upper"]
    FIT --> CMP["成分分解：trend / seasonal / holidays"]

推定は既定で MAP（mcmc_samples=0、cmdstanpy で最適化）。mcmc_samples>0 にするとフルベイズ（NUTS で季節の不確実性まで事後に反映）になりますが遅くなります。

2. 擬似系列で当てる：成分復元と予測区間

コード①：トレンド＋季節に Prophet を当て、成分と予測区間を出す

線形トレンド（傾き $0.15$ ）＋年次季節を仕込んだ月次系列に Prophet を当て、make_future_dataframe → predict で予測します。forecast から trend と yearly 成分を取り出して真値と比べ、yhat_lower/yhat_upper（95%区間）でホールドアウトを評価します。予測区間の MC 標本は np.random.seed(0) で固定し再現可能に。

import warnings, logging
warnings.simplefilter("ignore")
logging.getLogger("cmdstanpy").setLevel(logging.ERROR)
logging.getLogger("prophet").setLevel(logging.ERROR)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
import prophet
from prophet import Prophet
print(f"prophet version = {prophet.__version__}")

# 真の構造：線形トレンド + 年次季節(周期12) + 観測ノイズ
np.random.seed(3)
n_total, h = 132, 12
n = n_total - h
dates = pd.date_range("2014-01-01", periods=n_total, freq="MS")
t = np.arange(n_total)
trend_true = 10 + 0.15*t
seasonal_true = 6*np.sin(2*np.pi*t/12) + 3*np.cos(2*np.pi*2*t/12)
y_all = trend_true + seasonal_true + np.random.normal(0, 1.2, n_total)
df = pd.DataFrame({"ds": dates[:n], "y": y_all[:n]})
test = y_all[n:n_total]

# Prophet：加法分解（区分線形トレンド + フーリエ年次季節）をベイズで当てる
np.random.seed(0)                      # 予測区間のMC標本を再現可能に固定
m = Prophet(interval_width=0.95, yearly_seasonality=True,
            weekly_seasonality=False, daily_seasonality=False,
            seasonality_mode="additive")
m.fit(df)
future = m.make_future_dataframe(periods=h, freq="MS")
fc = m.predict(future)                 # yhat, yhat_lower, yhat_upper, trend, yearly ...

fut = fc.iloc[n:n_total]
yhat, lo, hi = fut["yhat"].values, fut["yhat_lower"].values, fut["yhat_upper"].values
print(f"予測RMSE = {np.sqrt(np.mean((yhat-test)**2)):.3f}  "
      f"95%予測区間カバレッジ = {int(np.sum((test>=lo)&(test<=hi)))}/{h}")

# 成分の復元（trend と yearly を真値と対比）
trend_hat = fc["trend"].values[:n]
seas_hat = fc["yearly"].values[:n]
slope_hat = (trend_hat[-1]-trend_hat[0])/(n-1)
print(f"トレンド成分RMSE vs 真値 = {np.sqrt(np.mean((trend_hat-trend_true[:n])**2)):.3f}  "
      f"推定傾き = {slope_hat:.4f}（真値 0.15）")
print(f"季節成分RMSE vs 真値     = {np.sqrt(np.mean((seas_hat-seasonal_true[:n])**2)):.3f}")

# 図：予測区間つき予測 と 成分分解（trend / yearly）
fig, ax = plt.subplots(3, 1, figsize=(9, 8))
ax[0].plot(dates[:n], df["y"], color="0.6", lw=1, label="訓練データ")
ax[0].plot(dates[n:n_total], test, color="k", lw=1.2, marker="x", ms=4, label="実測（検証）")
ax[0].plot(dates[n:n_total], yhat, color="C1", lw=2, label="Prophet 点予測 yhat")
ax[0].fill_between(dates[n:n_total], lo, hi, color="C1", alpha=0.25, label="95%不確実区間")
ax[0].axvline(dates[n-1], ls=":", color="k")
ax[0].set_ylabel("値"); ax[0].legend(fontsize=8); ax[0].set_title("Prophet：予測区間つき予測")
ax[1].plot(dates[:n], trend_true[:n], color="k", ls="--", lw=1.4, label="真のトレンド")
ax[1].plot(dates[:n], trend_hat, color="C0", lw=1.8, label="Prophet トレンド成分")
ax[1].set_ylabel("トレンド"); ax[1].legend(fontsize=8); ax[1].set_title("成分①トレンド")
ax[2].plot(dates[:n], seasonal_true[:n], color="k", ls="--", lw=1.4, label="真の季節")
ax[2].plot(dates[:n], seas_hat, color="C2", lw=1.6, label="Prophet 年次季節成分")
ax[2].set_xlabel("年月"); ax[2].set_ylabel("季節"); ax[2].legend(fontsize=8); ax[2].set_title("成分②年次季節")
plt.tight_layout(); plt.show()

出力：

prophet version = 1.3.0
予測RMSE = 1.444  95%予測区間カバレッジ = 9/12
トレンド成分RMSE vs 真値 = 0.149  推定傾き = 0.1541（真値 0.15）
季節成分RMSE vs 真値     = 0.515

出力の意味：Prophet は成分をきれいに復元します。トレンド成分の傾きは $0.1541$ （真値 $0.15$ ）、トレンド RMSE $0.149$ 、季節成分 RMSE $0.515$ ——forecast の trend・yearly 列を取り出すだけで、解釈できる分解が手に入ります（plot_components も同じ中身）。予測は yhat と yhat_lower/yhat_upper で一行。ただしカバレッジは 9/12（75%）と名目 95% より低めです。これは Prophet の既定（MAP）の不確実区間が、トレンド変化と観測ノイズは織り込むものの季節の推定誤差を含まないため。フルに不確実性を出すには mcmc_samples>0（フルベイズ）にします。この「既定区間はやや狭い」性質は実務で要注意の落とし穴です。

コード②：トレンド屈折（changepoint）を捉え、SARIMA と比較

途中で傾きが折れる区分線形トレンド（前半 $+0.5$ → 後半 $-0.2$ ）＋季節の系列で、Prophet が屈折を捉えるかを見ます。推定トレンドの前半・後半の傾きを測り、ホールドアウト RMSE を SARIMA と比較します（時間順分割、予測の評価指標と時系列CV）。

import warnings, logging
warnings.simplefilter("ignore")
logging.getLogger("cmdstanpy").setLevel(logging.ERROR)
logging.getLogger("prophet").setLevel(logging.ERROR)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
from prophet import Prophet
from statsmodels.tsa.statespace.sarimax import SARIMAX

# 真の構造：途中で傾きが屈折する区分線形トレンド + 季節 + ノイズ
np.random.seed(4)
n_total, h = 144, 24
n = n_total - h                       # 訓練 120
cp = 70                               # トレンド屈折点（傾きが切り替わる）
t = np.arange(n_total)
slope1, slope2 = 0.5, -0.2
trend_true = np.where(t < cp, 10 + slope1*t, 10 + slope1*cp + slope2*(t-cp))
seasonal_true = 5*np.sin(2*np.pi*t/12)
y_all = trend_true + seasonal_true + np.random.normal(0, 1.5, n_total)
dates = pd.date_range("2012-01-01", periods=n_total, freq="MS")
df = pd.DataFrame({"ds": dates[:n], "y": y_all[:n]})
test = y_all[n:n_total]

# Prophet：区分線形トレンド＋自動changepointが屈折を捉える
np.random.seed(0)                     # 予測区間のMC標本を固定
m = Prophet(interval_width=0.95, yearly_seasonality=True,
            weekly_seasonality=False, daily_seasonality=False,
            changepoint_prior_scale=0.1)
m.fit(df)
future = m.make_future_dataframe(periods=h, freq="MS")
fc = m.predict(future)
yhat_p = fc["yhat"].values[n:n_total]
lo, hi = fc["yhat_lower"].values[n:n_total], fc["yhat_upper"].values[n:n_total]
rmse_p = np.sqrt(np.mean((yhat_p - test)**2))

# 推定トレンドの前半・後半の傾き（屈折を捉えたか）
trend_hat = fc["trend"].values[:n]
slope_early = (trend_hat[cp-5]-trend_hat[5])/(cp-10)
slope_late = (trend_hat[n-1]-trend_hat[cp+5])/(n-1-cp-5)
print(f"推定トレンド傾き  前半 = {slope_early:.3f}（真 {slope1}）  "
      f"後半 = {slope_late:.3f}（真 {slope2}）")

# SARIMA と比較
sar = SARIMAX(df["y"].values, order=(1,1,1), seasonal_order=(0,1,1,12),
              trend="n").fit(disp=False)
yhat_s = sar.get_forecast(steps=h).predicted_mean
rmse_s = np.sqrt(np.mean((yhat_s - test)**2))

print(f"ホールドアウトRMSE  Prophet = {rmse_p:.3f}   SARIMA = {rmse_s:.3f}")
print(f"Prophet 95%予測区間カバレッジ = {int(np.sum((test>=lo)&(test<=hi)))}/{h}")

# 図：屈折を捉えたトレンドと両モデルの予測
fig, ax = plt.subplots(figsize=(9, 4.8))
ax.plot(dates[:n], df["y"], color="0.65", lw=1, label="訓練データ")
ax.plot(dates[n:n_total], test, color="k", lw=1.2, marker="x", ms=4, label="実測（検証）")
ax.plot(dates[:n], trend_hat, color="C0", lw=1.6, ls="--", label="Prophet 推定トレンド（屈折）")
ax.plot(dates[n:n_total], yhat_p, color="C1", lw=2, label="Prophet 予測")
ax.fill_between(dates[n:n_total], lo, hi, color="C1", alpha=0.22, label="Prophet 95%区間")
ax.plot(dates[n:n_total], yhat_s, color="C3", lw=2, ls=":", label="SARIMA 予測")
ax.axvline(dates[cp], ls=":", color="g", label="真の屈折点")
ax.axvline(dates[n-1], ls=":", color="k")
ax.set_xlabel("年月"); ax.set_ylabel("値"); ax.legend(fontsize=8, ncol=2)
ax.set_title("changepoint（トレンド屈折）を Prophet が捉える")
plt.tight_layout(); plt.show()

出力：

推定トレンド傾き  前半 = 0.512（真 0.5）  後半 = -0.188（真 -0.2）
ホールドアウトRMSE  Prophet = 1.710   SARIMA = 4.155
Prophet 95%予測区間カバレッジ = 21/24

出力の意味：Prophet の自動 changepoint が屈折を捉えました——推定トレンドの前半傾き $0.512$ （真値 $0.5$ ）、後半 $-0.188$ （真値 $-0.2$ ）。ラプラス事前が「屈折点でだけ勾配を折る」ので、向きが反転したトレンドを素直に表現できます。後半の傾きを正しく外挿するため、ホールドアウト RMSE は Prophet $1.710$ vs SARIMA $4.155$ と Prophet が勝ちます。SARIMA は差分（ $I(1)$ ）で「直近の局所勾配」を引き継ぐので、屈折後の系列でも前のモメンタムを引きずって誤差が膨らむ。ただしこれは「トレンド屈折が支配する系列だから」で、Prophet の区分線形という構造が真の構造に合致したのが勝因です。区間カバレッジは $21/24$ （ $87.5\%$ ）。逆に、トレンドが素直で残差に強い自己相関がある系列なら、それを直接モデル化する SARIMA が勝ちます。「データの主構造に合うモデルを選ぶ」——比較の教訓はいつもこれです。

3. 数式の直観

Prophet＝時間特徴への（ベイズ）回帰：自己回帰（過去の $y$ を使う）でなく、時刻 $t$ から作った特徴（トレンド基底・フーリエ・休日ダミー）への回帰。だから欠測があっても、その時刻の特徴さえあれば当てられる——欠測・不等間隔に強い理由です。
ラプラス事前＝スパースな屈折： $\delta_j\sim\text{Laplace}(0,\tau)$ は L1 罰則と同型で、多くの $\delta_j$ を 0 に潰す（正則化との関係）。だから候補 changepoint を多めに置いても、効くものだけ残る。 $\tau$ がトレンドの柔らかさのつまみ。
不確実性の3源：Prophet の予測区間は (1) 将来のトレンド屈折（過去の屈折頻度から将来の屈折を生成）、(2) 季節推定誤差（フルベイズ時のみ）、(3) 観測ノイズ、の3つから来ます。遠未来ほど (1) が効いて区間が広がる設計。既定 MAP では (2) が抜けるので狭めになりがち。
分解志向 vs 自己相関志向：Prophet/STL（STL分解）は「トレンド＋季節＋…」へ分解して各成分を当てる。ARIMA（季節ARIMA（SARIMA））は差分後の自己相関を AR/MA で当てる。前者は解釈・自動化・欠測耐性、後者は自己相関の効率的なモデル化が強み。

⚠️ よくある誤解・落とし穴

「Prophet の既定区間は正しいカバレッジ」ではない：既定は MAP（mcmc_samples=0）で、季節の不確実性が区間に入りません（コード①でカバレッジ 9/12）。フルに出すなら mcmc_samples=300 等でフルベイズに。いずれにせよホールドアウトでカバレッジ点検を（予測の評価指標と時系列CV）。
「Prophet は常に ARIMA より強い」ではない：コード②で勝ったのはトレンド屈折が主構造だったから。残差に自己相関が強い系列では SARIMA が上。汎用の万能薬ではありません。Prophet 自身、残差の自己相関は基本モデル化しません。
「changepoint_prior_scale は触らなくていい」ではない：小さすぎるとトレンドが鈍く屈折を見逃し、大きすぎると過剰に折れて将来トレンドが暴れます。屈折のあり方に応じて調整・交差検証を。
「点予測の遠未来トレンドは信頼できる」ではない：Prophet は最後の区間の勾配を直線外挿します。遠未来でトレンドが構造変化すれば外れる。区間は将来の屈折を生成して広げますが、過信は禁物（ベイズ構造時系列と同じ外挿の注意）。
「ロジスティックトレンドは自動」ではない：飽和する成長は growth="logistic" と cap（上限）を自分で指定します。既定の線形トレンドのまま頭打ちデータに当てると、青天井に外挿してしまいます。
「バージョン差を無視してよい」ではない：Prophet はバックエンド（pystan→cmdstanpy）や既定値が版で変わってきました（本ノートは prophet 1.3.0）。コードを使い回すときは prophet.__version__ と API を要最新確認。