RのCausalImpactをPythonで完全再現する

Google CausalImpact（R）と数値同等なPythonライブラリ。
GibbsサンプラーをRustで実装し、Rの30倍高速。TensorFlow不要。

既存Python実装（pycausalimpact, tfcausalimpact）はアルゴリズムがRと異なり、
推定結果が乖離する根本的な問題がありました。
本ライブラリはR bstsと同一のGibbsサンプリング・Spike-and-Slab変数選択を
Rustで忠実に再実装し、5シナリオ×6指標のCI自動テストでR互換性を保証します。

pip install bsts-causalimpact
https://github.com/YuminosukeSato/bsts-causalimpact

Transcript

1. bsts-causalimpact Bayesian Structural Time Series for Causal Inference in Python v1.5.1 | MIT License | Python 3.10+
2. 施策効果定量化のための反事実予測 Google CausalImpact(R)のPython移植。Rustで10-30x高速化 特徴 内容 R bsts完全互換 同一Gibbsアルゴリズム。CI自動検証で数値同等性を保証 Rust高速推論 PyO3経由のRust Gibbs sampler。Rの10-30倍の速度 Spike-and-Slab 多数の共変量から重要変数を自動選択 BSTS=Bayesian Structural Time Series, CI=Continuous Integration
3. CausalImpactの分析フロー 介入前後の時系列データから因果効果を推定する3ステップ Step 処理 説明 1 Pre期間でモデル構築 BSTSをGibbs samplerで推定 2 Post期間で反事実予測 介入なしの時系列を事後分布から生成 3 差分 = 因果効果 実測値と反事実予測の差分を効果量として定量 化 Brodersen et al. (2015). MCMC=Markov Chain Monte Carlo
4. 手法の適用条件と限界 以下の3条件を全て満たす場合にのみ有効 前提条件 内容 違反時の代替手法 コントロール非汚染 統制群が介入の影響を受けない DiD / Synthetic Control 関係安定性 被説明-統制間の関係が不変 構造変化検定 十分なPre期間 Post期間の3倍以上(目安) Pre期間の延長
5. 既存Python実装の根本的制約 既存Python移植は全てアルゴリズムか依存関係に問題を抱える ライブラリ アルゴリズム 根本的制約 pycausalimpact MLE (statsmodels) MCMCでなくMLEのため不確実性評価 が不正確 tfcausalimpact VI / HMC (TF) TF 3GB+必須。VIはGibbsと異なる近似 tfp-causalimpact TFP-based R数値検証未公開。Spike-and-slab未 実装 MLE=Maximum Likelihood Estimation, VI=Variational Inference, GB=Gigabytes
6. 主要ライブラリの機能・品質比較 R CI This tfp-ci tfCI pyCI Algorithm Gibbs/C++ Gibbs/Rust TFP VI/HMC MLE Depend. R, bsts np, pd, mpl TF 3GB+ TF 3GB+ statsm Spike-slab Yes Yes ? No No Seasonal Yes Yes ? Yes No Dyn. reg. Yes Yes ? No No Speed 2.1s 0.07s sec min sub-s R test Ref. CI auto None Visual None Speed: T=1000. TFP=TensorFlow Probability, HMC=Hamiltonian Monte Carlo
7. 本ライブラリの 3つの設計原則 原則 実現手段 R bsts完全互換 同一Gibbsアルゴリズム(Kalman + Sim. smoother)をRust実 装 TF不要の軽量依存 numpy/pandas/matplotlib のみ。計算コアをRustに集約 CI自動検証 Fixture-based + Live R比較 の2層CIで全コミット検証
8. Python-Rustハイブリッドアーキテクチャ Python API層 -> PyO3 -> Rust Gibbs sampler LTO=Link Time Optimization, API=Application Programming Interface
9. R対比10-30倍の推論速度 Rust LTO最適化 + Cholesky precision sampler T k niter Rust R (bsts) 倍率 100 0 1000 0.008s 0.213s 26x 500 0 1000 0.033s 0.997s 30x 1,000 0 1000 0.069s 2.108s 31x 1,000 5 1000 0.197s 2.171s 11x 5,000 0 1000 0.330s 10.264s 31x Median of 3 runs. python benchmarks/benchmark.py で再現可能
10. R bsts機能の互換状況一覧 Matching=CI検証済み(+-3%以内) / Supported=実装済み R bsts機能 Status 補足 Local level model Matching Kalman + Simulation smoother SdPrior (s.size=32) Matching InvGamma(16, 16*sigma^2) Spike-and-slab Matching Coordinate-wise + StudentSpikeSlabPrior Seasonal component Matching nseasons + season_duration Data standardization Matching Pre期間 mean/sd で標準化 prior_level_sd Matching デフォルト 0.01 (R同一) expected.model.size Matching デフォルト 2 expected.r2 / prior.df Matching 0.8 / 50 Dynamic regression Supported FFBS + random-walk beta Local linear trend Supported state_model設定で切替 FFBS=Forward Filtering Backward Sampling
11. 5シナリオ x 6指標のR数値同等性テスト R CausalImpact 1.4.1 (bsts 0.9.10, R 4.5) と比較 Scenario point effect cum effect ci lower ci upper rel effect p value basic +-3% +-3% +-1% +-1% +-3% a=0.05 covariates +-3% +-3% +-1% +-1% +-3% a=0.05 strong +-3% +-3% +-1% +-1% +-3% a=0.05 no_effect <2.0 <2.0 <2.0 <2.0 <0.5 a=0.05 seasonal +-1% +-1% +-1% +-1% +-1% a=0.05 seasonal は最も厳しい+-1%基準。no_effect は絶対値許容
12. 2層CI戦略によるリグレッション防止 全コミットでR数値同等性を自動検証 Laye r CI名 トリガー 内容 1 Fixture-based (ci.yml) 全PR/push R参照データとPython出力を比較 2 Live R comparison 週次自動 R 4.5でfixture再生成して比較 Layer 1: 高速(R不要) / Layer 2: R環境での完全再現を保証 PR=Pull Request ごとにLayer 1でブロッキング検証
13. 5行で完結する分析コード pip install bsts-causalimpact 後すぐ実行可能 import pandas as pd from causal_impact import CausalImpact data = pd.read_csv('data.csv', index_col='date', parse_dates=True) pre = ['2020-01-01', '2020-03-14'] post = ['2020-03-15', '2020-04-14'] ci = CausalImpact(data, pre, post) print(ci.summary()) ci.plot().savefig('result.png') data: 第1列=被説明変数、残りの列=共変量
14. summary()とreport()の出力例 Posterior inference表で効果量・CI・p値を一覧表示 Average Cumulative Actual 136.32 3953.19 Prediction (s.d.) 125.42 (0.66) 3637.07 (19.08) 95% CI [124.18, 126.71] [3601.33, 3674.59] Abs. effect (s.d.) 10.90 (0.66) 316.13 (19.08) 95% CI [9.61, 12.13] [278.60, 351.86] Rel. effect (s.d.) 8.69% (0.57%) 8.69% (0.57%) p = 0.001 / Posterior prob. of causal effect: 99.90%
15. 3パネル出力による効果の可視化 Original / Pointwise / Cumulative の3パネル構成 ci.plot() で生成。青=実測、破線=反事実予測、陰影=95%CI
16. 共変量なしでの単変量分析 Local levelモデルのみで介入効果を推定可能 共変量なしでもLocal levelモデルで反事実を推定
17. Spike-and-Slab変数選択の自動化 多数の共変量から重要変数を自動選択 expected_model_size で事前スパース性制御。posterior_inclusion_probs で確認
18. 季節性モデルへの対応 nseasons / season_duration でR AddSeasonal()互換 State-space seasonal: sum-to-zero制約の回転行列で表現
19. Dynamic Regressionによる時変係数推定 dynamic_regression=True で係数の時間変化を捕捉 Random-walk FFBS (Durbin-Koopman simulation smoother) で実装
20. R Vignetteデータでの完全再現 R公式チュートリアルと同一データ・同一結果 R CausalImpact vignette の ArrowAds データを使用
21. 導入方法と今後のロードマップ 項目 内容 Install pip install bsts-causalimpact Python 3.10+ License MIT GitHub github.com/YuminosukeSato/bsts-causalimpact Roadmap Rayon並列 / MCMC diagnostics / Prior sensitivity Contributions welcome! Issues / PRs on GitHub MIT=Massachusetts Institute of Technology