因果推論與貝氏結構時間序列模型
INFERRING CAUSAL IMPACT USING BAYESIAN STRUCTURAL TIME-SERIES MODELS (Google Inc. — 2015)
前言
**不想看廢話的可以直接往下拉到正文XD
因果推斷是現代分析中幾乎最重要的一環，隨著更多複雜模型的興起，人類已經難以辨識在建模中所經手的資訊是否有足夠的因果要素，這使得最後得出的結論很常因為不恰當的資訊使用與處理方法，而得出完全不一樣的結果。歷史上有許多經典的橋段演繹了這一問題，比如，著名的倖存者偏差問題，錯誤地下了結論認為機翼防護由於平均彈孔數量最多而最需要防護，忽略了實際上回收樣本已是倖存樣本這一偏誤問題。
另一種經典而廣為人知的情況則是相關性與因果性並不互相蘊含，比如在流行病學的研究下發現，接受激素療法的女性冠心病發作機率比一般人低，忽略了當時可接種的婦女通常社會地位較高，發作機率本就較低。(實際上這也是一種倖存者偏差的體現)
上面的結論聽起來頗為荒謬可笑，然而實際的案例與分析方法往往複雜地令人無可避免地掉入這類陷阱裡頭。舉例而言，大型的零售量販店想要測試新貨架擺放位置對於該類商品銷售額的增長效益，很顯然地我們無法單純用擺放改變前後的差異進行判斷，否則會忽略了由於時間變化所帶來的外部因素變化；又或是在電商的場景下你進行了一個基於 Logistic Regression 的消費者終生價值預測任務，得到結論說明年齡層高的消費者傾向有更多的人均消費，而這可能也是不正確的，因為這些人的瀏覽行為本質上可能就不同，或是有特定離群值的存在造成該關聯性沒有辦法在廣義範疇上都通用 (當然你可以做更多更嚴謹的檢驗，但隨著你想應用更複雜的模型得到精確結果時，檢驗的難度也隨之上升)。
當然，現時最可行的方法似乎是 A/B Testing，許多基於其概念衍生出來的方法也能更加快速適應這個大環境，比如 Netflix 提到的快速線上評估方法 :
然而雖然說有長足的進步，相較於蓬勃發展的建模技術還是顯得因果推論這一塊有些單薄，因此今天我想來介紹一篇 Google 在 2015 年發表的好文章，連結與名稱如下 :
INFERRING CAUSAL IMPACT USING BAYESIAN STRUCTURAL TIME-SERIES MODELS (2015)
https://projecteuclid.org/download/pdfview_1/euclid.aoas/1430226092Predicting the Present with Bayesian Structural Time Series (2013)
https://storage.googleapis.com/pub-tools-public-publication-data/pdf/41335.pdf
這篇文章提出了一種叫作 Bayesian Structural Time-Series Model (BSTSM) 的方法，用以推斷變因造成的時間性演變，並通過貝氏方法整合先驗知識，童時能兼容多種因素來源如趨勢、季節性等。但在開始介紹以前，我們先來對何謂因果推論進行一些回顧。
回顧一 : 因果推論的關鍵要素
因果推論當然是統計學裡面舉足輕重的一環，實際上包含了許多的先備知識與研究，但這邊我想要很快速的總結前人成果，並嘗試用簡單、好理解的方式說明因果推論須包含哪些關鍵因素，參考書目主要為 :
Element of Casual Inferencing (2017)
https://library.oapen.org/handle/20.500.12657/26040Casual Inference in Statistics : An Overview (2009)
https://projecteuclid.org/euclid.ssu/1255440554http://projecteuclid.org/euclid.ssu/1255440554
根據基因學家 Sewall Wright 的觀點，我們可以透過一些方程或是圖的方式來描述因果關係，假設我們先以單變量、線性方程的情況來描述因果關係，應該寫成 :
前面 y = Bx 非常好理解，描述 y 這個變量，無論是連續或是離散的標籤值，某種程度上可以被寫成是 x 變量乘上某個權重的結果，而後面的 uY 則代表著是一切的外在因素，或也可以理解成自然情況下 y 的發生結果，當 X 不變時該值仍在影響著 Y (注意這裡的概念並不完全等同於常數項)。而既然 y 變量有一個 uY 在影響其結果，相對來說 X 也應該有個 uX 代表所有其他因素影響著他的結果，這些變數我們統一稱作外生變量 (Exogeneous Variables)。
但是只透過方程來描述顯然不夠正確，如果我們調整等式位置改成 :
難不成這是在說 y 同時也影響著 x ? 顯然不對，我們缺乏了一些路徑的方向性，如果我們視覺化該方向性，實際上應該會像是 :
其中，箭頭所代表的係數 Beta，一般我們會寫成共變數，也就是 Cov(X, Y) = B，我們要確保導致 X 與 Y 的外生變量沒有互相影響時，也等同於確保 Cov(Ux, Uy) = 0，此時便可更加確保 X 對 Y 因果影響的獨立性。這類的寫法不證自明，且這種方程搭配圖的表達方式非常易於拓展，比方說你可以寫成以下多因素的模式 :
但這樣的方式顯然不夠，問題在於變數增加時的拓展仍舊十分繁瑣，且線性的表達並不足以完美反應真實情況，因此基於機率學的表達方法因應而生，旨在抽離原先用係數的方式來反應因果影響的效果，進而重新定義該效果為一種將改變發送至不同變量的能力，聽起來很拗口，我們來看實際例子，以上圖 A 的三種情況來說，我們應該寫成 :
此時，Z 在 Fy(.) 中不存在的這一情況，反應了 Z 的變動並不使 Y 變動，只要 X 與 Uy 也一併保持不動，當我們用這樣的方式來描述變量之間的交互關係時，我們就說這些方程有結構性，每個方程並不因其他方程的形式(線性、非線性、高維低維) 而有所轉移，就可以很好地透過某些變化是否對變量有影響，而描繪出其因果關係，並便於我們後續用機率的方式進行推演。
機率的推斷等不等同於因果的推斷呢 ? 即機率的推斷可能通常是經過某些統計學習而得來的，比如對一個分配的期望值的求取，而因果推斷僅單純是針對觀察到的情況、中間干擾因素及最終結果，做出的一個決定性判斷，這兩者之間是等價的嗎 ?
當我們能夠將機率推斷所應用的數值，嚴格來說為不同變量的 Joint Distribution，歸入 (subsume) 因果推斷所需的觀察值與結果標籤時，我們同步也能把機率模型歸入因果推斷模型之中，進而加強我們對於兩者等價的信心，當然，凡事無法盡善盡美，但仍舊有一些原則我們可以遵循。
很顯然聯繫程度(Association, Correlation)是最基本我們會想到的，當 P(Y | X) > P(Y)，我們會開始對 X 與 Y 的因果關係抱持一些期待，而這樣的理解也是上面提及許多偏誤的產生，畢竟發生 X 的情況下 Y 的機率變高，實際上也可能是一種倒果為因，比方說我們可能發現智商高的人有較高的比例是名校學生，那有沒有可能是進了名校才讓他們智商變高? You never know。因此，以下還有兩個原則值得我們多注意 :
Intervention 干預
很多時候我們並不只希望看見 P(Y | X)，更希望知道的是如果人為干預 X 對 Y 的影響情況如何，也就是 P(Y|do(X)) 的情況，如果人工干預對另一方有效，我們就不必在糾結於兩者之間的共同關係，而是可以抽離出來單獨看成 兩個不同的機制 : 當我透過某種辦法調高一群人的智商時，這些人進入名校的機率變高了，而當我想辦法給予一群人實際的名校學位時，他們的智商並不見得同步變高，這時我們才可以推論說，起始智商與進入名校是有因果關係存在的。
而在統計的角度，如果講得抽象一些，干預結果會導致因果效應程度的分佈改變，而干預該因果效應本身並不使結果的分佈改變，換言之如果你干預中名校的機率(假設以名校學生/總學生來說)，那對於智商是否影響成為名校學生這件事的程度分佈可能有所改變，比如會變成 :
很好理解，比方因為當你改變了名校的定義，那這時候當然因為大家都容易進名校了，此時 IQ 與名校的因果關係分佈就不如以前。但反過來說，干預 IQ 的分佈並不使得進名校的機率分佈有所改動，這是因為不論該影響的分佈如何變化，仍舊是只有具備一定IQ的人更有機會進名校。
CounterFactuals 反事實條件 (虛擬條件)
設想一個情況 (這邊我就直接按照原文翻譯，有興趣的讀者請參照 Elements in Casual Inference 的 3.3)，有一個神奇的眼疾療法存在，對於 99% 的病人來說，該療法有效且病人被治癒了，而如果不做該治療，病人會在一天內失明；反之，對於剩下的 1% 來說，該治療會使他們在一天內失明，而如果不做治療，他們的視力會自然恢復。當然，醫生在開始前可不知道病人是屬於 99% 或是那 1%。
今天有一個眼疾病人走進來了，在醫生給予其療程後一天內失去了視力，很明顯我們並沒有辦法說明任何因果關係，畢竟病人在不做任何處理的情況下，也會自然失去視力，我們如何下定論說是因為醫生錯誤的處置導致，抑或是醫生處置得當，只是療程沒有效益? 此時，我們可以問這樣的一個問題 : 「如果醫生沒有給予療程，該病人會如何?」 很顯然地，該病人並不會因此失去視力，因此，我們方可以得到結論說醫生處置不當。 (當然，現實情況我們並不能這樣說，畢竟先決條件已知被治癒的機率要高得多)
這道理也應用在前面的例子上，當我們用某種方法提高某人智商，最終他進了名校，我們並不可以馬上下結論說智商一定與進名校有因果聯繫，畢竟不做任何事，他還是有可能進名校，但如果我們發現不提高其智商，他進不了名校的時候，某種程度上我們可以推斷這樣的因果關係，在統計的角度可能則是當有一定樣本的時候，出現顯著的機率下降。
回顧二 : 時間序列的因果推論方法
前面雖然講了這麼多假象案例，實際的情況卻更多的是，當我們做了某種行動，可能是行銷策略、新功能的上線、使用者介面呈現的調整等等，我們想知道是否該行動跟我們所要監測的指標 (比如用戶增長、停留時長、人均消費金額等) 有因果聯繫。回顧剛才講過的兩大原則，Intervention & Counterfactual，即做了會怎樣 vs. 不做會怎樣，自然會將我們的思維引導到 A/B Testing 上去，對測試組別做該活動並且對控制組別不做該活動，此外須確保測試組別與控制組別一開始是可以比較的。
在一個固定、離散的情況下，這種實驗比較好處理，舉例說我有 N 個使用者給予其原本的 DM 跟新設計的 DM，詢問其偏好，這時我們可以隨機地將其分成一半一半去做實驗，頂多是後續再人為確保某些特定的指標上兩個組別分佈一致，比方說不希望有一組年齡層平均特別高等等。然而，在更實際的場景來說，我們沒有辦法在如此離散的條件下區分實驗組與對照組，像是我們要測定的是以店面為單位的活動設計，或是雖然有使用者，我們只能得到其消費紀錄，無從判斷起其離散指標，這時候，我們的問題就會變成在連續的時間序列指標上進行 A/B Testing。
如上圖所示，而當中最關鍵的一環是，如何在連續值的維度下，還能準確抓出可比較的 Test 與 Control ? 隨機抓取當然是一個最直觀的方式，但這衍生出兩個問題 : 一、我如何確定兩者是可比較的；二、如果兩者因為各種因素而不可比較，如何進行修正。
一、時間序列的可比較性檢驗
我們有一系列的方式可以進行統計檢定，比方說進行一些相似度的檢定或是相關係數的監測，比較常見的則是通過 ANOVA 的方法去衡量兩者序列的平均值是否相等，但這有一個缺點是並不能反映序列的時間特質，比方說如果兩組序列正好在平均值上互相顛倒，那麼他們的序列平均數一致且標準差也一致，穩穩地通過 ANOVA 測試；或是在一個序列之中出現了 Outlier，可能該 Outlier 不會對整體結果有太大影響，但在實際情況中我們特別在意這些 Outlier 以確保資料蒐集正確且沒有外部因素影響。
有一些方法可以改進這類檢驗，其中一種是首先我們透過某種方式將時間序列拆解成 Trends, Seasonal, Residual 序列，比如透過基於 LOESS 的 STL 方法等等，然後只對 Trends, Seasonal 進行 ANOVA，並確保 Residual 是隨機即可，這通常代表你不是這麼在意序列中的一些微幅波動，而是關注在整體趨勢與週期性的方面兩者是否完美 match :
比較複雜的方法則像是動態時間規整 (Dynamic Time Warping)，對於時序資料的變形如平移、上下震盪或隨機震盪等十分敏感 :
DTW 方法的核心思想十分簡單，也就是計算兩時序資料每個點的距離，比較不同的是其通過一個嶄新概念叫做規整路徑距離 (Warping Path Distance) 來評估兩序列的相似度，而計算的方法則是採用動態規劃 (Dynamic Programming)。對於兩個長度不同的時間序列 X = {x1, x2, …, xn}; Y = {y1, y2, …, ym} 來說，我們可以計算其規整路徑 W = {w1, w2, …, wk}，其中 :
max(|X|, |Y|) <= K <= |X| + |Y|規整路徑是甚麼概念呢? 想像一下我們順著任一個時間序列上的不同資料點往下走，比如說 x1, x2, …，我們希望能夠對每一個點找到另一個時間序列 Y 上面對應的最近點，點距離的計算就是歐幾里得距離，這樣計算的好處就是可以確保序列的時間性，以免在點對點的 matching 中失去其順序因素。那麼，如何找到這條最佳的規整路徑 ? 這就是我們剛才提到的動態規畫概念了，可以直接看如下的圖片 :
簡單來說，當我順著一條路徑往下走，以 i 代表序列 X 順序且 j 代表序列 Y 的順序，因為 i, j 一定是單調遞增的，因此對於一個序列資料 (i, j)，其最佳路徑一定是來自於三個點 : (i-1, j-1), (i-1, j) 或 (i, j-1)，這時就可以在 O(NM) 的時間複雜度下完成計算。
最終，得到這樣的路徑後，便可以使用 ANOVA 檢驗的方法檢測其與完美擬和路徑 (理論最短距離為 0) 的差異，此時由於不會因為變形而造成檢驗失真，評估的效果往往較好，當然，可比較性本身沒有公定的標準可以判斷，大多數情況下仍就需要不少主觀判斷來幫助決策，導致不論再怎麼擬和，總是可以找到某種檢定讓其失敗，因此建議僅需針對當前的時間成本、可用工具等考量進行盡可能擬真即可，無須追求完美。
二、不可比較的時間序列修正
實際上這個問題應該被轉化為，我如何確保我用以模擬測試組表現的控制組在干預前的效果足夠好? 如果控制組與測試組在干預前差異過大，我們都應該優先思考重新抽樣、模擬的可能性，因為此時一切的修正都會遭遇到人工與現實差異的風險，而往往我們都會先執行隨機抽樣的方式去得出兩個組別，展開其時序資料藉以評估抽樣效果，然而問題就在於，往往由於各式各樣的因素，隨機抽樣的難度比起離散樣本要大得多，這時我們有哪些方式可以改善這樣的情況呢?
首先，永遠都應該檢測隨機抽樣的機制是否真如理想情況，許多人會使用用戶 id 的奇偶數來進行分割，然而這樣的方式往往在執行了幾次活動後效果越來越差 : 奇數群組永遠被選定為實驗組而出現行為偏差，或是在抽樣的過程中由於每個時間段的人數不同，造成最終時間序列出現不同期間的長相不同等情況，可能就要考慮基於時間段的電腦隨機抽樣等等。
不過，真實情況總是沒有這麼令人順心，舉例來說消費者本身的差異偏大，即使通過隨機抽樣也難以盡如人意，或是像量販店希望衡量特定活動效果，而由於店面本身差距過大而難以找到控制組等等問題，以下就來介紹兩種方法改善抽樣方法或人工合成抽樣樣本。
【Clustering Based Cohort Analysis】
許多時候我們更希望得出的測試組與對照組在特定指標上的分佈一致，這時我們就可以考慮基於 Clustering 的方式進行抽樣，這改善了我們如果僅針對單一參數分佈通過劃分區間抽樣時，隨著參數數量增多，可能造成後續劃分受到前面劃分的影響而無法達到理想效果等問題。
Clustering 的方法非常多，如 K-Means, Spectral Clustering, T-SNE 等等，而一般而言對於 Financial Patterns 來說，只要經過統計檢驗，一般的 K-Means 已經十分足夠，你往往可以從最簡單將其轉化為離散特徵開始，比如用其平均消費金額、成為用戶時長等進行初步 Cluster，如果需要一次用到整個序列進行 Cluster，則可以考慮特定 embedding 方法用少數維度表達長序列情況等等。而當我們做完 Clustering 以後，由於在分析上大多數還是需要將多種 Clusters 整合成測試與對照兩組，因此後續仍需要手動整合，整合的方式也不困難，比如加權平均就是一種簡單的方法，也就是說在拆分 Cluster 並兩兩進行對比後，根據該 Cluster 的人數將增長合理的還原出來 :
【Synthetic Control 合成控制組】
有時候問題則出在根本沒有現成的控制組可以使用，比如剛才提到的在某特定店面進行行銷活動，或是政府在新制訂計畫時選定特定縣市做為實驗對象，很顯然你根本難以找到一組可以比較的控制對象進行 Counterfactual Inferrence。這時你就可以參考在計量經濟領域常用的合成控制組方法進行分析，這邊分享 David 的文章，我認為寫得十分簡潔易懂，就不再多說了 :
** 備註 : 在行銷上用 Synthetic Control 也是一種過度樂觀的估計，SC 的思想本質上就是人為去挑選那些跟實驗組長得像的控制組出來當對照，這本身就有透過限縮條件來滿足分析需要的意涵在裡面。
【合成方法與因果推論】
其實有些讀者應該看出來了，實際上 Synthetic Control 就是透過某種方法去模擬出一組真實樣本來，可以想像成我有一些規則可以預測整段時間序列，我們將這個預測出來的時間序列資料拿來做為控制對象，那麼，我們可不可以用一些更具代表性的預測方法來進行實作呢? 比如將測試組在干預前的時間序列資料作為輸入值，通過如 ARIMA, GARCH, LSTM, TCN 等經典模型進行預測，就可以得到更準確的控制組別了。
Well，是，但也不是。你當然可以通過其他方法去得出控制組，但不論是剛才的 Synthetic Control 與其他經典模型，本質上都是有一些假說在背後的，比如 ARIMA, GARCH 的許多前提都在於如序列平穩、觀測週期夠長，這些在經濟或金融數據上也許可行，放到離散的行銷場景上就行不通的前提，使得透過這些方法做出預測的結果不大可信。進一步說，哪怕透過這些方法得出了可靠的控制組序列，後續進行分析的方式往往使用 DiD (Difference-in-difference) 方法，而 DiD 的限制在於 : 一、幾組序列真的是要互相獨立的 (which 在行銷上有太多可能性打破該假設)，對於序列相關的資料，容易誇飾其影響造成我們對結果過度樂觀；二、大多數應用場景下使用者成為控制組或對照組並不是參照一個固定時間的，而是隨著產品迭代而不斷加入新資訊，這時候就不適合用 DiD 分析；最後一點則是，上面說的 SCM, ARIMA 實際上都可以當作是一種 convex combination 的方式，比如在原始 SCM 的論文就提到用權重加總為 1 的限制條件去挑選控制對象，先天上排除了非 convex 的可能性 (which 又一個理想但不務實的假設)，而真的要做 LSTM, TCN，更大的關卡在於從頭到尾就沒有控制對象的 input 輸入，你很難界定 overfit 的程度是否影響以及多大程度影響分析結果等等。
因此，本文要介紹另外一種基於貝氏的時間序列方法，基於貝氏機率的計算可以有效彙總後驗知識，並且能利用序列中離散或是 emission process (像是趨勢序列、季節性序列) 等資訊，此外，由於是貝氏方法，因此能夠給出預測的不確定性，同步避免了 overfitting。
貝氏結構時間序列 (Bayesian Structural Time Series)
好的，又是一個寫得有夠長的前言，which 鐵定是流失了許多讀者XD 要是我是自媒體經營者肯定餓死但還好不是XD
該方法的提出是在 2013 年的一篇論文中 :
Predicting the Present with Bayesian Structural Time Series (2013)
https://storage.googleapis.com/pub-tools-public-publication-data/pdf/41335.pdf在論文的簡介就很明確地提到，其實很多經濟數據的蒐集都是非常不乾淨的，比如下游就是quarterly的回報，回報期間不一樣長等等，此外這些數據即使被 report 了，也常被修改，導致許多現代的計量經濟學很多就是希望在這些數據被 report 之前，基於我們已經觀測到的資訊，包含過去的序列(序列本身的演化)以及一些離散訊號(外在因子的影響)，做出 Nowcasting (預測當下)。
BSTS 一共包含兩個部分 : 第一是 time-series based 的方法去捕捉序列中的趨勢以及季節性訊號；第二則是通過其他影響因子對序列的影響能力進行迴歸分析，由於這兩個部分是可加的，我們可以很簡單透過 Bayesian 方法得到其 joint model。
1. Structural Time Series
這個部份很好理解，由兩個部分組成 :
首先，y 當然是我們觀測到的值，我們知道一定有某種隱含狀態 α 在影響他，因此我們將其寫成線性組合的形式表達，也就是式 (1)。而式 (2) 則說明隱含狀態本身也會隨時間改進，受到上一個隱含狀態的影響以及外在因素影響 (到這邊有點像 markov process)。而上面的 Z, T, R 則包含了已知值與未知參數 (通常是寫成 0 與 1)，用這樣的表達形式來描述一個模型，我們稱其有 state space form (狀態空間形式)。
為甚麼要寫成這樣的形式 ? 如果你對時間序列分解有一些概念，你一定知道一組序列可以用如下方式表達 (Basic Structural Model) :
由上到下分別是 :
1. 原始觀測序列 : y = μ + τ + βx + ϵ
2. 趨勢水平(Level of trend): μ = μ_last + slope_last(δ) + constant
3. 趨勢斜率(Slope of trend) : δ = δ_last + constant
4. 季節項(Seasonality) : τ = zero expectation of the seasonality terms in last period + constant至於 βx 在下個 section 會說明，就是我們上述提到的離散資訊對序列影響的回歸估計。一般而言，求解這個問題會使用 Kalman Filter (卡爾曼濾波)，基礎思想如我剛才提到是基於 Markov Chain，這邊就不再多談，只說明跟該算法有關的部分 :
該 Filter 每次透過序列觀測資訊對於隱含狀態推斷的先驗機率，得出下一個時間段隱含狀態的分布 :
注意由於這邊的 Formulation 都是透過 Gaussian Distribution，因此上述分布其實就是個多變量的常態分佈而已。
2. Spike-and-slab Regression
對於剛才提到的 βx ，代表其外部因子對序列的影響，當然這個因子數目可以很大，事實上在論文裡頭就提到大多數時候可能的外部因子數量是遠超過序列觀測數目的，比如在 Google Search Trend 的預測任務中，透過用戶 query 的文字來構造外部因子，可想而知會面臨高度稀疏的問題。
要透過貝氏方法來表達這個問題，一種作法為通過 spike and slab prior，這東西是幹嘛的 ? 想像在一個回歸問題裡面，參數一定會帶著某種 prior distribution，這跟一般的回歸問題我們透過梯度求解最小損失函數的方式來找最佳參數組合不同，在 Bayesian Regression 裏頭，對於參數組合的估計是通過從先驗分布與觀測值去構造後驗分布，再不斷往下迭代的結果 (將參數分布的預測區間不斷優化)，而對於參數過多難以估計的問題，在梯度求解裡面我們使用的是如 Lasso 之類帶正則化的方式來挑選變量，在 Bayesian Regression 裡頭就變成是我們將問題轉化為 :
多出來的參數 γ 也會自帶先驗分布，通常我們是假設 Bernoulli 分布，也就是非 0 即 1，這時候我們針對其作後驗分布推斷的時候，如果某參數 x 不重要，那麼 γ 為 0，自動幫我們扮演了篩選參數的角色，而以上的這些過程就是一個 Spike and Slab Regression。因此，我們的問題表示為 :
由於我們的假設是獨立分布，因此寫成連乘的形式，而 γ 由於是 Bernoulli，因此寫成 :
πk 可以視需求，一般工程便利上，直接假設為一樣的數值即可，該數值的決定可以視工程人員覺得模型大小應該限制在怎樣的範圍而決定，越高的值代表越高的機率 γ = 1，模型的參數也會隨之增多。對以上的過程有興趣請參考 :
The Bayesian Lasso (2008)
http://www.math.chalmers.se/Stat/Grundutb/GU/MSA220/S16/bayeslasso.pdf當然，以上還有一些分布的細節計算我這邊就不一一推導了，主要部分就是不斷的求 given 前一參數下的 marginal posterier，這邊直接列公式 :
其中 :
這些過程也許現在看起來很複雜 (對實際上也是挺複雜的)，但可以用圖的形式來表達 :
由上圖可以看到各個狀態與參數視怎麼被逐步估計，以及基於 Bayesian 的出的 Counterfactual 推斷與從中得出的 Prediction (Inference) 部分。說明至此，我們就可以用如下的形式去迭代求解 :
1. 通過 p(α|y, θ, β, σ^2) 模擬隱含狀態 α，採用 Kalman Filter 或論文提及為了特殊需求改寫的 simulation smoother from Durbin and Koopman2. 擬合參數 θ ∼ p(θ|y, α, β, σ^2)3. 通過 Markov Chain 的方式模擬 (β, σ) ~ p(β, σ^2 |y, α, θ)
當然，今天的重點也不是該模型怎麼求解或有哪些變體，而是在於用其來做因果推斷，因此，這邊就簡化這些流程了，具體的使用可以參照論文作者自己開源的 R Package (Python 的我目前還沒找到，後面有找到補上或哪天看能不能簡單寫一下 XD) :
BSTS 因果推斷
實際上這篇的重點在於因果推斷，只是不知道為何先備知識一寫就沒完沒了，只能說 Bayesian 的東西雖然強大但現今非主流，因此許多東西真的得要從頭看過一遍。
那麼，怎麼用 BSTS 做因果推斷呢 ? 實際應用起來也非常簡單，直接將其 Prediction 的部分當作控制組就行了，如下圖 :
此時我們關注的點在於序列的 post-period 與 pre-period 之間的增量效應 (incremental effect)，也就是 :
由於 Bayesian Spike-and-slab 的 Formulation 方法，此時我們可以很簡單地透過無干預因子，也就是其中 βx 的 β = 0 方式來篩選 (當然也不一定是 0 啦，看你怎麼 Formulate 這一切 XD。當然看上去簡單，實際上我們還是需要思考一些問題，也就是同樣都是預測序列作為對照組，用 BSTS 的方法就是比較好，他能夠解決哪些問題?
一、分離序列因素與外在影響因素
如本文前面所說，因果推論的重點在於僅考量干預作用造成的影響，過往用 ARIMA 等方法合成的控制組由於缺乏將序列本身的發展因素如趨勢、季節性與外在變量整合在一起的方案 (比如 ARIMA 考慮了前者但壓根沒有將 Latent State 考量在裡頭)，因此有可能 overestimate 或 underestimate 干預因子表現，而 Synthetic Control 就更不用說，從頭到尾只是對於控制組的一個權重篩選，挑出 Post-period 比較像的控制對象而已。
二、對於先驗序列的要求
由於改成使用 Bayesian Inference 的方式實作，推斷的結果不需要滿足一些條件如序列平穩等等，可以直接在參數分布上去優化區間，在許多實際的問題上效果可能優於 ARIMA, TCN 等方法 (當然反之就有可能更差)。
不過 Bayesian 為人詬病的地方在於需要花費更多的時間去研究先驗分布的設計，這可能意味著你需要花更多力氣在數據清洗、探索性數據分析、數學建模 (將數據合理地轉成分布形式且要易於推導)，這部分在現今追求減小特徵工程負擔的機器學習界來說可以說是背道而馳。
三、預測隨機過程
如果你採用一些非時間序列特化的算法，比如 multiple linear regression (which synthetic control 由於有真實 control group 反而不會有這問題)，是沒有辦法在足夠長的預測時間內還能給出 Stochastic Processes 的，這多半是因為 Serial Autocorrelation 沒有被衡量在模型裏頭。
四、放鬆 Convex 條件
相對於 Synthetic Control 來說，Bayesian 是無須考量 Convex 問題的 (梯度求解與後驗分布求解的差異)，這種最大的好處就在於我無需考量輸入參數的 scale 問題，如同決策樹一般由於下次的分裂的衡量基準與 scale 無關，不需要做一些 Normalization 或是 EDA 的操作。
五、區間估計
由於 Bayesian Prediction 最後是給出一個區間的預測範圍，因此在分析上也可以同步得出 Significant Lift 與 Insignificant Lift，這樣的估計有時候也能避免我們過度樂觀地就下決策，在保守的情境下我們也可以等到出現 Significant Lift 後再進行大規模佈署。
以上這些因素都是 BSTS 的優勢，至於在實際應用中具體需要用哪種方式，其實還是見仁見智，各個方法其實缺乏一個公定的比較基準在於何者才能最好地反應真實 control group 應該有的情況，因此考量的點都變成是工程的實現性考量與假設的合理性考量二者。
結論
其實 Bayesian 的東西對我而言一直都是相對困難的，困難的點倒不是在於用怎麼樣的數學運算去得出結論，而是在於邏輯與數學形式之間的強連結性，從本篇的例子其實可以看得出來 Bayesian 以及其與 Markov Chain, State-Space 等想法之所以結合，都是因為在邏輯推斷上這樣的撰寫合理，我們再用合理的數學形式去表達，相對於傳統機器學習或是深度學習在已經給定一個欲優化目標的前提下去推導要來得複雜許多。
當然，這也是 Bayesian 之所以被認為是因果推論一大神器的關鍵，甚至有人認為該領域是下一次機器學習的突破口，也就是讓 AI 具備強因果推論能力，而該方法可能就是由 Bayesian 的方式透過充分利用先驗去推斷後驗的過程，來不斷加強機器對於各式新事件的適應能力。
關於本篇文章由於篇幅實在過長，因此中間省略了不少推導過程，也請各位讀者如果文章內容有誤不吝賜教，雖說 medium 文章的中心思想一直都是以自學為主，但也很歡迎各方前輩後進們的共同交流！
