🍎 連享會主頁：lianxh.cn
🎦 2021 文本分析·爬蟲·機器學習📅 2021 年 3.17-4.11 (三個周末)🔑 主講：遊萬海 (福州大學)；司繼春 (上海對外經貿大學)🍓 課程主頁：https://gitee.com/lianxh/text
⏩ 2021 生存分析專題 (Survival Aanlysis)⌚ 2021 年 4.24-25 (周六、周日)⭐ 主講：王存同教授 (中央財經大學)⛪ 課程主頁：https://gitee.com/lianxh/ST
New！ lianxh 命令發布了： GIF 動圖介紹隨時搜索 Stata 推文、教程、手冊、論壇，安裝命令如下： . ssc install lianxh
作者：賓健博 (中山大學)郵箱：binjb@mail2.sysu.edu.cn
目錄
1. 引言
2. 貝葉斯估計
3. 伯努利實驗
4. Stata 中貝葉斯估計
5. MCMC 方法和線性回歸的貝葉斯估計
6. 參考資料
7. 相關推文
近年來，貝葉斯估計 (Bayesian estimation) 在宏觀經濟分析和複雜的計量經濟模型中應用越來越廣泛。比如，以貝葉斯估計為基礎的 VAR (向量自回歸) 模型和 DSGE (動態隨機一般均衡) 模型的前沿研究，就經常出現在主流經濟學期刊上。因此，掌握貝葉斯方法逐漸成為了宏觀計量研究者的必修課。
遺憾的是，一般對 VAR 模型和 DSGE 模型的貝葉斯估計都是在 MATLAB 中完成。當然，我們也可以利用 Stata 內置的貝葉斯估計命令處理一些較簡單的模型。
基於上述考慮，本文將簡要介紹貝葉斯估計和 MCMC 方法的基本概念、如何在 Stata 中使用 bayes 命令前綴、標準估計命令 bayesmh 及其自帶的交互菜單。
在統計學和計量經濟學中，參數估計都扮演著非常重要的角色。頻率學派認為概率是在長期重複的隨機實驗中產生的性質。但在貝葉斯學派看來，概率是不確定性的一種表達。這使得我們可以對無法重複的事件的概率進行考察，但相應的問題是將主觀性引入了概率之中。反映到參數估計上，頻率學派認為模型參數是一個定值，應當使用觀測數據去接近它；而貝葉斯學派認為模型參數是一個隨機變量，有自己的分布，應當使用先驗信息並結合觀測值去修正它的分布。
因此，貝葉斯估計的核心任務就是結合先驗信息和觀測值信息，獲得模型參數的概率分布 (後驗分布，posterior)。我們將所有觀測值記為 ，將模型參數記為 。在估計之前，需要確定先驗分布 (prior) 和似然函數：
在確定了先驗分布和似然函數後，可以使用貝葉斯定理 (Bayes' Law) 求解出後驗分布：
其中：
在實際的貝葉斯分析中，由於 是常數，一般會忽略這一項，只考慮後驗分布的分子，即：
在本小節，我們將介紹最簡單的貝葉斯估計實例——伯努利實驗。假設現有一枚硬幣，將其拋起，然後估計正面朝上的概率 。
根據前面的討論，我們需要明確先驗分布和似然函數。由於似然函數不具有主觀性，故先予以分析。假設拋了 次硬幣，其中正面朝上的次數為 ，則似然函數為：
再來看先驗分布。在實際操作中，我們往往選擇似然函數的共軛先驗 (conjugate prior)。由於共軛先驗具有與似然函數相乘後分布族保持不變的特徵，我們可以通過簡單的數學推導求解出 。在本例中，似然函數是由伯努利分布 (Bernoulli distribution) 推導而來，故其對應的共軛先驗為 Beta 分布。
在選定了先驗分布的分布族之後，還需要調整分布的參數以符合先驗信念。Beta 分布有兩個參數 和 ，通過改變兩者的大小，其形狀會發生很大的變化。在這裡，我們主要關注兩種情況：
，即均勻分布 (uniform distribution)。均勻分布意味著參數在每一個值的概率密度相同，也稱為無信息先驗 (non-informative prior) 或者 扁平先驗 (flat prior)。如果選用這個分布，就代表完全沒有參數的先驗信息；
。由下圖，可以發現當 時， 分布會出現 「中間高，兩頭低」 的情況，且概率密度的最大值出現在 0.5 處。如果選用這樣的分布，意味著我們認為模型參數 的可能性最高，這和我們一般的常識相符 (拋硬幣正面朝上的概率應該在 50% 附近)。同時還可以發現 和 越大，概率密度越向中間集中，因此 和 可以作為衡量我們先驗信息確信程度的參數。這裡選擇的 20 並不是唯一的，代表了我們對先驗信息是比較確信的。
在 Stata 中，我們可以使用上述先驗分布對原問題進行估計，具體代碼如下：
. clear //清空工作區. set seed 20 //設定隨機種子. set obs 10000 //設定觀測值數量，在這裡 n=10000. gen y = runiform(0,1) > 0.5 //生成服從伯努利分布的樣本，記為 y. bayesmh y, likelihood(dbernoulli({theta})) prior({theta}, beta(1,1)) //估計命令Burn-in ...Simulation ...Model summary--Likelihood: y ~ bernoulli({theta})Prior: {theta} ~ beta(1,1)--Bayesian Bernoulli model MCMC iterations = 12,500Random-walk Metropolis-Hastings sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 10,000 Acceptance rate = .442Log marginal likelihood = -6935.8231 Efficiency = .1992 -- | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval]--+---- theta | .4991434 .0051053 .000114 .499113 .4892574 .5092233--可以看出，當設定先驗分布為 時，平均值的點估計結果為 0.4991434，和真實的數據生成過程中的參數 0.5 很接近。如果選擇具有較強先驗信息的 ，同樣的設定下點估計值達到了 0.4992197，這證明了先驗信息對後驗分布的影響。 4. Stata 中貝葉斯估計上例中，核心代碼如下：. bayesmh y, likelihood(dbernoulli({theta})) prior({theta}, beta(1,1))其中，likelihood() 指定了似然函數的形式，在本例中似然函數由伯努利分布導出，同時我們將參數稱為 theta (在本例中只有一個參數，因此不需要進行區分，只是單純的命名，在更複雜的例子中需要通過名字區分參數)。prior() 則指定了先驗分布，本例中將 theta 的先驗分布定為 。上述內容的編程實現起來並不困難，但是在更複雜的問題中，根據似然函數和先驗分布的不同，很多時候需要指定額外的參數或選項。此時，貝葉斯估計命令 bayesmh 並不適合，需要通過 Stata 自帶的交互菜單生成。具體來看：我們可以通過點擊 Stata 上方的 Statistics 菜單，選擇底部的 Bayesian analysis 選項，再根據具體問題選擇相應選項，打開貝葉斯估計的交互界面。由於上文的例子不涉及回歸分析，選擇 General estimation and regression。由於模型是一個簡單的獨立實驗，每一個觀測值都是一次獨立實驗的結果，因此在 Syntax 處選擇 Univariate distributions，Dependent variable 處選擇生成的變量 y。Distribution 用於確定似然函數的分布，在本例中選擇伯努利分布，同時指定參數名為 theta。最後需要指定模型的先驗分布，首先點擊 Create 按鈕，然後在彈出的菜單中選擇相應的分布族和參數，便可以生成一個想要的先驗分布。最終結果如下：界面截圖點擊 OK 鍵，系統便會自動開始估計，可以得到前文中的結果。為了說明交互界面的便利性，我們考慮一個稍微複雜一些的例子。假設樣本不再是由伯努利分布生成，而是由二項分布 生成，每一個觀測值代表進行了 10 次成功率為 0.3 的伯努利實驗，觀測值序列記為 。我們依舊可以導出其似然函數：在這個例子中，原分布和似然函數中有兩個參數 和 ，其中， 代表二項分布的第一個參數， 代表二項分布的第二個參數。並且，參數 是預先給定的，需要在使用交互界面時需要輸入這個參數。沿用之前的先驗分布 (二項分布的共軛先驗也是 Beta 分布)，此時的代碼、結果和菜單操作如下：. clear //清空工作區. set seed 20 //設定隨機種子. set obs 10000 //設定觀測值數量，在這裡 n=10000. gen z = rbinomial(10,0.3) //生成服從伯努利分布的樣本，記為 z. bayesmh z, likelihood(dbinomial({theta},10)) prior({theta}, normal(0,10000))Burn-in ...Simulation ...Model summary--Likelihood: z ~ binomial({theta},10)Prior: {theta} ~ beta(1,1)--Bayesian binomial model MCMC iterations = 12,500Random-walk Metropolis-Hastings sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 10,000 Acceptance rate = .4526Log marginal likelihood = -17700.763 Efficiency = .2339 -- | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval]--+---- theta | .2987419 .0014841 .000031 .2987601 .2959134 .3016467--可以看到，當似然函數為二項分布之後，在菜單界面中多了一個 Bernoulli trials 的輸入框，也就是上文指定的 。在更複雜的問題和模型中，可能會要求輸入更多的參數和其他內容，此時菜單界面可以很好地代替複雜的語句指令，免去了繁雜的語法記憶。不僅如此，我們還可以通過鍵入命令 db postest 打開後驗分析的交互菜單進行分析。常用的功能如 Bayesian analysis 下的 Graphical summaryies and convergence diagnostics 菜單可以繪製很多與後驗分布相關的圖，Hypothesis testing using model posterior probabilities 和 Interval hypothesis testing 菜單可以利用生成的後驗分布進行假設檢驗。 5. MCMC 方法和線性回歸的貝葉斯估計在前面的例子中，可能會有讀者感到疑惑，當先驗分布選為共軛先驗時，後驗分布是有解析表達式的，這意味著其概率密度函數應該是確定的。但在結果顯示的界面中，明顯有模擬 (simulation)、迭代 (iteration) 和抽樣 (sampling) 這樣的常見於數值方法的字眼。這其實是因為 Stata 並非使用解析方法對後驗分布進行求解。為了進一步深入了解現代貝葉斯方法的操作流程，我們有必要對 MCMC (Markov chain Monte Carlo，馬爾科夫鏈蒙特卡洛) 方法進行介紹。對於 MCMC 方法，維基百科定義為：馬爾可夫鏈蒙特卡洛方法包括一類用於從概率分布中採樣的算法。通過構建具有所需分布作為其平衡分布的馬爾可夫鏈，可以通過記錄鏈中的狀態來獲得所需分布的樣本。這一段話有些抽象，實際上從應用層面上來說，我們只需要了解兩個最常用的 MCMC 方法：Metropolis–Hastings 算法 (簡稱 M-H 算法) 和吉布斯採樣 (Gibbs sampling)。前者是 Stata 內置的對後驗分布進行抽樣的算法，其主要作用是通過後驗分布的概率密度之比構造了一個具有接受-拒絕機制的馬爾科夫鏈，並通過它抽取出所需的後驗分布。這一算法也廣泛應用於 DSGE 模型的貝葉斯估計中，但它不是本文的重點。在宏觀計量經濟學中最常用的 (主要是因為其應用於 VAR 模型的貝葉斯估計) MCMC 方法是吉布斯採樣。為了讓大家對吉布斯採樣的應用背景有一個簡單的了解，下面將以線性回歸為例進行介紹。假設有如下的線性回歸模型：其中， 是一個 的矩陣， 是一個 的矩陣， 是一個 的矩陣。其似然函數的表達式如下：這裡有兩個待估參數，即係數矩陣 和擾動項方差 。分別考慮以下三種情況：在第一種情況中，由於擾動項方差已知，此時的共軛先驗是正態分布，因此假定係數矩陣 的先驗分布為 。先驗分布：似然函數：後驗分布：可以發現後驗分布仍然是正態分布：在第二種情況中，由於似然函數中的給定變量改變，其共軛先驗也隨之發生變化。查表知此時擾動項方差 應服從倒伽馬分布 (Inverse Gamma distribution)，具體運算過程和第一種情況類似，可以獲得相同分布族的後驗分布 。在實際中，更加常見的是第三種情況，即兩個待估參數都未知。此時，沒有辦法獲得單一的共軛先驗，一種解決辦法是設定聯合先驗分布：其中， 的分布仍是無條件的，和第二種情況中的設定一樣，但 則是有條件的，我們將其設定為：在這樣的設定下，聯合先驗分布還是一個正態分布，和似然函數的分布族一致，這種情況稱之為自然共軛先驗 (natural conjugate prior)。需要注意的是，最後解出的後驗分布是一個聯合後驗分布 ，但我們想考察的後驗分布則是邊緣後驗分布 和 。在自然共軛先驗的設定下，可以通過積分的方法求得邊緣後驗分布的解析形式。但這也對先驗分布的選擇施加了很強的限制。在較複雜的模型中，我們往往不希望限制先驗分布的選擇，但非自然共軛先驗難以通過積分求解出邊緣後驗分布，因此需要應用吉布斯採樣。吉布斯採樣的思想其實很簡單，當有了兩個能條件分布的抽取器時 (即第一種和第二種情況)，可以通過迭代的方法獲得逼近邊緣分布的抽樣。在本例中，具體操作為：設定係數矩陣 的無條件先驗分布和擾動項方差 的無條件先驗分布，並從 的無條件先驗分布中抽取一個樣本；根據第一種情況的做法，以第 1 步中抽取的 為條件，從獲得的後驗分布中抽取一個 的樣本；根據第二種情況的做法，以第 2 步中抽取的 為條件，從獲得的後驗分布中抽取一個 的樣本；重複第 2 步和第 3 步，獲得一個經驗樣本，可用於估計邊緣後驗分布。這一操作步驟看似複雜，但在 Stata 中的實現其實很容易，只需要在線性回歸的命令前加上 bayes: 即可，但需要注意的是，此時回歸命令不能縮寫，並且如果要用吉布斯採樣，需要使用 bayes,gibbs:。. clear //清空工作區. set seed 20 //設定隨機種子. sysuse auto //使用系統自帶的 auto 數據集. bayes, gibbs: regress price mpg weight //執行吉布斯採樣的貝葉斯估計Burn-in ...Simulation ...Model summary----Likelihood: price ~ normal(xb_price,{sigma2})Priors: {price:mpg weight _cons} ~ normal(0,10000) (1) {sigma2} ~ igamma(.01,.01)----(1) Parameters are elements of the linear form xb_price.Bayesian linear regression MCMC iterations = 12,500Gibbs sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 74 Acceptance rate = 1 Efficiency: min = .9457 avg = .9864Log marginal likelihood = -696.85658 max = 1 ---- | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval]----+----price | mpg | -5.282435 27.88942 .278894 -5.410969 -59.53671 49.12283 weight | 2.074039 .198759 .001988 2.076242 1.682404 2.464258 _cons | 1.928181 99.76441 .997644 2.99603 -202.2644 194.9716----+---- sigma2 | 6448502 1099995 11311.5 6331241 4641478 8928296----Note: Default priors are used for model parameters. 6. 參考資料溫馨提示： 文中連結在微信中無法生效。請點擊底部「閱讀原文」。M-H 算法和共軛先驗的維基百科 -Link1- -Link2-Blake A P, Mumtaz H. Applied Bayesian econometrics for central bankers[J]. Technical Books, 2012. -Link-Koop G, Korobilis D. Bayesian multivariate time series methods for empirical macroeconomics[M]. Now Publishers Inc, 2010. -Link-Bill Rising. Bayesian Analysis using Stata. 2015. -Link- 7. 相關推文Note：產生如下推文列表的命令為： lianxh VAR DSGE 分布, m安裝最新版 lianxh 命令： ssc install lianxh, replace溫馨提示： 文中連結在微信中無法生效。請點擊底部「閱讀原文」。Stata程序：Monte-Carlo-模擬之產生符合特定分布的隨機數VaR 風險價值: Stata 及 Python 實現 🎦 2021 文本分析·爬蟲·機器學習🍓 課程主頁：https://gitee.com/lianxh/text ⏩ 2021 生存分析專題⛪ 課程主頁：https://gitee.com/lianxh/ST 關於我們🍎 連享會 ( 主頁：lianxh.cn ) 由中山大學連玉君老師團隊創辦，定期分享實證分析經驗。👉 直達連享會：【百度一下：連享會】即可直達連享會主頁。亦可進一步添加 主頁，知乎，面板數據，研究設計 等關鍵詞細化搜索。連享會主頁 lianxh.cn New！ lianxh 命令發布了： GIF 動圖介紹隨時搜索連享會推文、Stata 資源，安裝命令如下： . ssc install lianxh使用詳情參見幫助文件 (有驚喜)： . help lianxh
可以看出，當設定先驗分布為 時，平均值的點估計結果為 0.4991434，和真實的數據生成過程中的參數 0.5 很接近。如果選擇具有較強先驗信息的 ，同樣的設定下點估計值達到了 0.4992197，這證明了先驗信息對後驗分布的影響。
上例中，核心代碼如下：
. bayesmh y, likelihood(dbernoulli({theta})) prior({theta}, beta(1,1))其中，likelihood() 指定了似然函數的形式，在本例中似然函數由伯努利分布導出，同時我們將參數稱為 theta (在本例中只有一個參數，因此不需要進行區分，只是單純的命名，在更複雜的例子中需要通過名字區分參數)。prior() 則指定了先驗分布，本例中將 theta 的先驗分布定為 。上述內容的編程實現起來並不困難，但是在更複雜的問題中，根據似然函數和先驗分布的不同，很多時候需要指定額外的參數或選項。此時，貝葉斯估計命令 bayesmh 並不適合，需要通過 Stata 自帶的交互菜單生成。具體來看：我們可以通過點擊 Stata 上方的 Statistics 菜單，選擇底部的 Bayesian analysis 選項，再根據具體問題選擇相應選項，打開貝葉斯估計的交互界面。由於上文的例子不涉及回歸分析，選擇 General estimation and regression。由於模型是一個簡單的獨立實驗，每一個觀測值都是一次獨立實驗的結果，因此在 Syntax 處選擇 Univariate distributions，Dependent variable 處選擇生成的變量 y。Distribution 用於確定似然函數的分布，在本例中選擇伯努利分布，同時指定參數名為 theta。最後需要指定模型的先驗分布，首先點擊 Create 按鈕，然後在彈出的菜單中選擇相應的分布族和參數，便可以生成一個想要的先驗分布。最終結果如下：界面截圖點擊 OK 鍵，系統便會自動開始估計，可以得到前文中的結果。為了說明交互界面的便利性，我們考慮一個稍微複雜一些的例子。假設樣本不再是由伯努利分布生成，而是由二項分布 生成，每一個觀測值代表進行了 10 次成功率為 0.3 的伯努利實驗，觀測值序列記為 。我們依舊可以導出其似然函數：在這個例子中，原分布和似然函數中有兩個參數 和 ，其中， 代表二項分布的第一個參數， 代表二項分布的第二個參數。並且，參數 是預先給定的，需要在使用交互界面時需要輸入這個參數。沿用之前的先驗分布 (二項分布的共軛先驗也是 Beta 分布)，此時的代碼、結果和菜單操作如下：. clear //清空工作區. set seed 20 //設定隨機種子. set obs 10000 //設定觀測值數量，在這裡 n=10000. gen z = rbinomial(10,0.3) //生成服從伯努利分布的樣本，記為 z. bayesmh z, likelihood(dbinomial({theta},10)) prior({theta}, normal(0,10000))Burn-in ...Simulation ...Model summary--Likelihood: z ~ binomial({theta},10)Prior: {theta} ~ beta(1,1)--Bayesian binomial model MCMC iterations = 12,500Random-walk Metropolis-Hastings sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 10,000 Acceptance rate = .4526Log marginal likelihood = -17700.763 Efficiency = .2339 -- | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval]--+---- theta | .2987419 .0014841 .000031 .2987601 .2959134 .3016467--可以看到，當似然函數為二項分布之後，在菜單界面中多了一個 Bernoulli trials 的輸入框，也就是上文指定的 。在更複雜的問題和模型中，可能會要求輸入更多的參數和其他內容，此時菜單界面可以很好地代替複雜的語句指令，免去了繁雜的語法記憶。不僅如此，我們還可以通過鍵入命令 db postest 打開後驗分析的交互菜單進行分析。常用的功能如 Bayesian analysis 下的 Graphical summaryies and convergence diagnostics 菜單可以繪製很多與後驗分布相關的圖，Hypothesis testing using model posterior probabilities 和 Interval hypothesis testing 菜單可以利用生成的後驗分布進行假設檢驗。 5. MCMC 方法和線性回歸的貝葉斯估計在前面的例子中，可能會有讀者感到疑惑，當先驗分布選為共軛先驗時，後驗分布是有解析表達式的，這意味著其概率密度函數應該是確定的。但在結果顯示的界面中，明顯有模擬 (simulation)、迭代 (iteration) 和抽樣 (sampling) 這樣的常見於數值方法的字眼。這其實是因為 Stata 並非使用解析方法對後驗分布進行求解。為了進一步深入了解現代貝葉斯方法的操作流程，我們有必要對 MCMC (Markov chain Monte Carlo，馬爾科夫鏈蒙特卡洛) 方法進行介紹。對於 MCMC 方法，維基百科定義為：馬爾可夫鏈蒙特卡洛方法包括一類用於從概率分布中採樣的算法。通過構建具有所需分布作為其平衡分布的馬爾可夫鏈，可以通過記錄鏈中的狀態來獲得所需分布的樣本。這一段話有些抽象，實際上從應用層面上來說，我們只需要了解兩個最常用的 MCMC 方法：Metropolis–Hastings 算法 (簡稱 M-H 算法) 和吉布斯採樣 (Gibbs sampling)。前者是 Stata 內置的對後驗分布進行抽樣的算法，其主要作用是通過後驗分布的概率密度之比構造了一個具有接受-拒絕機制的馬爾科夫鏈，並通過它抽取出所需的後驗分布。這一算法也廣泛應用於 DSGE 模型的貝葉斯估計中，但它不是本文的重點。在宏觀計量經濟學中最常用的 (主要是因為其應用於 VAR 模型的貝葉斯估計) MCMC 方法是吉布斯採樣。為了讓大家對吉布斯採樣的應用背景有一個簡單的了解，下面將以線性回歸為例進行介紹。假設有如下的線性回歸模型：其中， 是一個 的矩陣， 是一個 的矩陣， 是一個 的矩陣。其似然函數的表達式如下：這裡有兩個待估參數，即係數矩陣 和擾動項方差 。分別考慮以下三種情況：在第一種情況中，由於擾動項方差已知，此時的共軛先驗是正態分布，因此假定係數矩陣 的先驗分布為 。先驗分布：似然函數：後驗分布：可以發現後驗分布仍然是正態分布：在第二種情況中，由於似然函數中的給定變量改變，其共軛先驗也隨之發生變化。查表知此時擾動項方差 應服從倒伽馬分布 (Inverse Gamma distribution)，具體運算過程和第一種情況類似，可以獲得相同分布族的後驗分布 。在實際中，更加常見的是第三種情況，即兩個待估參數都未知。此時，沒有辦法獲得單一的共軛先驗，一種解決辦法是設定聯合先驗分布：其中， 的分布仍是無條件的，和第二種情況中的設定一樣，但 則是有條件的，我們將其設定為：在這樣的設定下，聯合先驗分布還是一個正態分布，和似然函數的分布族一致，這種情況稱之為自然共軛先驗 (natural conjugate prior)。需要注意的是，最後解出的後驗分布是一個聯合後驗分布 ，但我們想考察的後驗分布則是邊緣後驗分布 和 。在自然共軛先驗的設定下，可以通過積分的方法求得邊緣後驗分布的解析形式。但這也對先驗分布的選擇施加了很強的限制。在較複雜的模型中，我們往往不希望限制先驗分布的選擇，但非自然共軛先驗難以通過積分求解出邊緣後驗分布，因此需要應用吉布斯採樣。吉布斯採樣的思想其實很簡單，當有了兩個能條件分布的抽取器時 (即第一種和第二種情況)，可以通過迭代的方法獲得逼近邊緣分布的抽樣。在本例中，具體操作為：設定係數矩陣 的無條件先驗分布和擾動項方差 的無條件先驗分布，並從 的無條件先驗分布中抽取一個樣本；根據第一種情況的做法，以第 1 步中抽取的 為條件，從獲得的後驗分布中抽取一個 的樣本；根據第二種情況的做法，以第 2 步中抽取的 為條件，從獲得的後驗分布中抽取一個 的樣本；重複第 2 步和第 3 步，獲得一個經驗樣本，可用於估計邊緣後驗分布。這一操作步驟看似複雜，但在 Stata 中的實現其實很容易，只需要在線性回歸的命令前加上 bayes: 即可，但需要注意的是，此時回歸命令不能縮寫，並且如果要用吉布斯採樣，需要使用 bayes,gibbs:。. clear //清空工作區. set seed 20 //設定隨機種子. sysuse auto //使用系統自帶的 auto 數據集. bayes, gibbs: regress price mpg weight //執行吉布斯採樣的貝葉斯估計Burn-in ...Simulation ...Model summary----Likelihood: price ~ normal(xb_price,{sigma2})Priors: {price:mpg weight _cons} ~ normal(0,10000) (1) {sigma2} ~ igamma(.01,.01)----(1) Parameters are elements of the linear form xb_price.Bayesian linear regression MCMC iterations = 12,500Gibbs sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 74 Acceptance rate = 1 Efficiency: min = .9457 avg = .9864Log marginal likelihood = -696.85658 max = 1 ---- | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval]----+----price | mpg | -5.282435 27.88942 .278894 -5.410969 -59.53671 49.12283 weight | 2.074039 .198759 .001988 2.076242 1.682404 2.464258 _cons | 1.928181 99.76441 .997644 2.99603 -202.2644 194.9716----+---- sigma2 | 6448502 1099995 11311.5 6331241 4641478 8928296----Note: Default priors are used for model parameters. 6. 參考資料溫馨提示： 文中連結在微信中無法生效。請點擊底部「閱讀原文」。M-H 算法和共軛先驗的維基百科 -Link1- -Link2-Blake A P, Mumtaz H. Applied Bayesian econometrics for central bankers[J]. Technical Books, 2012. -Link-Koop G, Korobilis D. Bayesian multivariate time series methods for empirical macroeconomics[M]. Now Publishers Inc, 2010. -Link-Bill Rising. Bayesian Analysis using Stata. 2015. -Link- 7. 相關推文Note：產生如下推文列表的命令為： lianxh VAR DSGE 分布, m安裝最新版 lianxh 命令： ssc install lianxh, replace溫馨提示： 文中連結在微信中無法生效。請點擊底部「閱讀原文」。Stata程序：Monte-Carlo-模擬之產生符合特定分布的隨機數VaR 風險價值: Stata 及 Python 實現 🎦 2021 文本分析·爬蟲·機器學習🍓 課程主頁：https://gitee.com/lianxh/text ⏩ 2021 生存分析專題⛪ 課程主頁：https://gitee.com/lianxh/ST 關於我們🍎 連享會 ( 主頁：lianxh.cn ) 由中山大學連玉君老師團隊創辦，定期分享實證分析經驗。👉 直達連享會：【百度一下：連享會】即可直達連享會主頁。亦可進一步添加 主頁，知乎，面板數據，研究設計 等關鍵詞細化搜索。連享會主頁 lianxh.cn New！ lianxh 命令發布了： GIF 動圖介紹隨時搜索連享會推文、Stata 資源，安裝命令如下： . ssc install lianxh使用詳情參見幫助文件 (有驚喜)： . help lianxh
其中，likelihood() 指定了似然函數的形式，在本例中似然函數由伯努利分布導出，同時我們將參數稱為 theta (在本例中只有一個參數，因此不需要進行區分，只是單純的命名，在更複雜的例子中需要通過名字區分參數)。prior() 則指定了先驗分布，本例中將 theta 的先驗分布定為 。
上述內容的編程實現起來並不困難，但是在更複雜的問題中，根據似然函數和先驗分布的不同，很多時候需要指定額外的參數或選項。此時，貝葉斯估計命令 bayesmh 並不適合，需要通過 Stata 自帶的交互菜單生成。具體來看：
我們可以通過點擊 Stata 上方的 Statistics 菜單，選擇底部的 Bayesian analysis 選項，再根據具體問題選擇相應選項，打開貝葉斯估計的交互界面。由於上文的例子不涉及回歸分析，選擇 General estimation and regression。
由於模型是一個簡單的獨立實驗，每一個觀測值都是一次獨立實驗的結果，因此在 Syntax 處選擇 Univariate distributions，Dependent variable 處選擇生成的變量 y。Distribution 用於確定似然函數的分布，在本例中選擇伯努利分布，同時指定參數名為 theta。最後需要指定模型的先驗分布，首先點擊 Create 按鈕，然後在彈出的菜單中選擇相應的分布族和參數，便可以生成一個想要的先驗分布。最終結果如下：
點擊 OK 鍵，系統便會自動開始估計，可以得到前文中的結果。
為了說明交互界面的便利性，我們考慮一個稍微複雜一些的例子。假設樣本不再是由伯努利分布生成，而是由二項分布 生成，每一個觀測值代表進行了 10 次成功率為 0.3 的伯努利實驗，觀測值序列記為 。我們依舊可以導出其似然函數：
在這個例子中，原分布和似然函數中有兩個參數 和 ，其中， 代表二項分布的第一個參數， 代表二項分布的第二個參數。並且，參數 是預先給定的，需要在使用交互界面時需要輸入這個參數。沿用之前的先驗分布 (二項分布的共軛先驗也是 Beta 分布)，此時的代碼、結果和菜單操作如下：
. clear //清空工作區. set seed 20 //設定隨機種子. set obs 10000 //設定觀測值數量，在這裡 n=10000. gen z = rbinomial(10,0.3) //生成服從伯努利分布的樣本，記為 z. bayesmh z, likelihood(dbinomial({theta},10)) prior({theta}, normal(0,10000))Burn-in ...Simulation ...Model summary--Likelihood: z ~ binomial({theta},10)Prior: {theta} ~ beta(1,1)--Bayesian binomial model MCMC iterations = 12,500Random-walk Metropolis-Hastings sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 10,000 Acceptance rate = .4526Log marginal likelihood = -17700.763 Efficiency = .2339 -- | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval]--+---- theta | .2987419 .0014841 .000031 .2987601 .2959134 .3016467--可以看到，當似然函數為二項分布之後，在菜單界面中多了一個 Bernoulli trials 的輸入框，也就是上文指定的 。在更複雜的問題和模型中，可能會要求輸入更多的參數和其他內容，此時菜單界面可以很好地代替複雜的語句指令，免去了繁雜的語法記憶。不僅如此，我們還可以通過鍵入命令 db postest 打開後驗分析的交互菜單進行分析。常用的功能如 Bayesian analysis 下的 Graphical summaryies and convergence diagnostics 菜單可以繪製很多與後驗分布相關的圖，Hypothesis testing using model posterior probabilities 和 Interval hypothesis testing 菜單可以利用生成的後驗分布進行假設檢驗。 5. MCMC 方法和線性回歸的貝葉斯估計在前面的例子中，可能會有讀者感到疑惑，當先驗分布選為共軛先驗時，後驗分布是有解析表達式的，這意味著其概率密度函數應該是確定的。但在結果顯示的界面中，明顯有模擬 (simulation)、迭代 (iteration) 和抽樣 (sampling) 這樣的常見於數值方法的字眼。這其實是因為 Stata 並非使用解析方法對後驗分布進行求解。為了進一步深入了解現代貝葉斯方法的操作流程，我們有必要對 MCMC (Markov chain Monte Carlo，馬爾科夫鏈蒙特卡洛) 方法進行介紹。對於 MCMC 方法，維基百科定義為：馬爾可夫鏈蒙特卡洛方法包括一類用於從概率分布中採樣的算法。通過構建具有所需分布作為其平衡分布的馬爾可夫鏈，可以通過記錄鏈中的狀態來獲得所需分布的樣本。這一段話有些抽象，實際上從應用層面上來說，我們只需要了解兩個最常用的 MCMC 方法：Metropolis–Hastings 算法 (簡稱 M-H 算法) 和吉布斯採樣 (Gibbs sampling)。前者是 Stata 內置的對後驗分布進行抽樣的算法，其主要作用是通過後驗分布的概率密度之比構造了一個具有接受-拒絕機制的馬爾科夫鏈，並通過它抽取出所需的後驗分布。這一算法也廣泛應用於 DSGE 模型的貝葉斯估計中，但它不是本文的重點。在宏觀計量經濟學中最常用的 (主要是因為其應用於 VAR 模型的貝葉斯估計) MCMC 方法是吉布斯採樣。為了讓大家對吉布斯採樣的應用背景有一個簡單的了解，下面將以線性回歸為例進行介紹。假設有如下的線性回歸模型：其中， 是一個 的矩陣， 是一個 的矩陣， 是一個 的矩陣。其似然函數的表達式如下：這裡有兩個待估參數，即係數矩陣 和擾動項方差 。分別考慮以下三種情況：在第一種情況中，由於擾動項方差已知，此時的共軛先驗是正態分布，因此假定係數矩陣 的先驗分布為 。先驗分布：似然函數：後驗分布：可以發現後驗分布仍然是正態分布：在第二種情況中，由於似然函數中的給定變量改變，其共軛先驗也隨之發生變化。查表知此時擾動項方差 應服從倒伽馬分布 (Inverse Gamma distribution)，具體運算過程和第一種情況類似，可以獲得相同分布族的後驗分布 。在實際中，更加常見的是第三種情況，即兩個待估參數都未知。此時，沒有辦法獲得單一的共軛先驗，一種解決辦法是設定聯合先驗分布：其中， 的分布仍是無條件的，和第二種情況中的設定一樣，但 則是有條件的，我們將其設定為：在這樣的設定下，聯合先驗分布還是一個正態分布，和似然函數的分布族一致，這種情況稱之為自然共軛先驗 (natural conjugate prior)。需要注意的是，最後解出的後驗分布是一個聯合後驗分布 ，但我們想考察的後驗分布則是邊緣後驗分布 和 。在自然共軛先驗的設定下，可以通過積分的方法求得邊緣後驗分布的解析形式。但這也對先驗分布的選擇施加了很強的限制。在較複雜的模型中，我們往往不希望限制先驗分布的選擇，但非自然共軛先驗難以通過積分求解出邊緣後驗分布，因此需要應用吉布斯採樣。吉布斯採樣的思想其實很簡單，當有了兩個能條件分布的抽取器時 (即第一種和第二種情況)，可以通過迭代的方法獲得逼近邊緣分布的抽樣。在本例中，具體操作為：設定係數矩陣 的無條件先驗分布和擾動項方差 的無條件先驗分布，並從 的無條件先驗分布中抽取一個樣本；根據第一種情況的做法，以第 1 步中抽取的 為條件，從獲得的後驗分布中抽取一個 的樣本；根據第二種情況的做法，以第 2 步中抽取的 為條件，從獲得的後驗分布中抽取一個 的樣本；重複第 2 步和第 3 步，獲得一個經驗樣本，可用於估計邊緣後驗分布。這一操作步驟看似複雜，但在 Stata 中的實現其實很容易，只需要在線性回歸的命令前加上 bayes: 即可，但需要注意的是，此時回歸命令不能縮寫，並且如果要用吉布斯採樣，需要使用 bayes,gibbs:。. clear //清空工作區. set seed 20 //設定隨機種子. sysuse auto //使用系統自帶的 auto 數據集. bayes, gibbs: regress price mpg weight //執行吉布斯採樣的貝葉斯估計Burn-in ...Simulation ...Model summary----Likelihood: price ~ normal(xb_price,{sigma2})Priors: {price:mpg weight _cons} ~ normal(0,10000) (1) {sigma2} ~ igamma(.01,.01)----(1) Parameters are elements of the linear form xb_price.Bayesian linear regression MCMC iterations = 12,500Gibbs sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 74 Acceptance rate = 1 Efficiency: min = .9457 avg = .9864Log marginal likelihood = -696.85658 max = 1 ---- | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval]----+----price | mpg | -5.282435 27.88942 .278894 -5.410969 -59.53671 49.12283 weight | 2.074039 .198759 .001988 2.076242 1.682404 2.464258 _cons | 1.928181 99.76441 .997644 2.99603 -202.2644 194.9716----+---- sigma2 | 6448502 1099995 11311.5 6331241 4641478 8928296----Note: Default priors are used for model parameters. 6. 參考資料溫馨提示： 文中連結在微信中無法生效。請點擊底部「閱讀原文」。M-H 算法和共軛先驗的維基百科 -Link1- -Link2-Blake A P, Mumtaz H. Applied Bayesian econometrics for central bankers[J]. Technical Books, 2012. -Link-Koop G, Korobilis D. Bayesian multivariate time series methods for empirical macroeconomics[M]. Now Publishers Inc, 2010. -Link-Bill Rising. Bayesian Analysis using Stata. 2015. -Link- 7. 相關推文Note：產生如下推文列表的命令為： lianxh VAR DSGE 分布, m安裝最新版 lianxh 命令： ssc install lianxh, replace溫馨提示： 文中連結在微信中無法生效。請點擊底部「閱讀原文」。Stata程序：Monte-Carlo-模擬之產生符合特定分布的隨機數VaR 風險價值: Stata 及 Python 實現 🎦 2021 文本分析·爬蟲·機器學習🍓 課程主頁：https://gitee.com/lianxh/text ⏩ 2021 生存分析專題⛪ 課程主頁：https://gitee.com/lianxh/ST 關於我們🍎 連享會 ( 主頁：lianxh.cn ) 由中山大學連玉君老師團隊創辦，定期分享實證分析經驗。👉 直達連享會：【百度一下：連享會】即可直達連享會主頁。亦可進一步添加 主頁，知乎，面板數據，研究設計 等關鍵詞細化搜索。連享會主頁 lianxh.cn New！ lianxh 命令發布了： GIF 動圖介紹隨時搜索連享會推文、Stata 資源，安裝命令如下： . ssc install lianxh使用詳情參見幫助文件 (有驚喜)： . help lianxh
可以看到，當似然函數為二項分布之後，在菜單界面中多了一個 Bernoulli trials 的輸入框，也就是上文指定的 。在更複雜的問題和模型中，可能會要求輸入更多的參數和其他內容，此時菜單界面可以很好地代替複雜的語句指令，免去了繁雜的語法記憶。
不僅如此，我們還可以通過鍵入命令 db postest 打開後驗分析的交互菜單進行分析。常用的功能如 Bayesian analysis 下的 Graphical summaryies and convergence diagnostics 菜單可以繪製很多與後驗分布相關的圖，Hypothesis testing using model posterior probabilities 和 Interval hypothesis testing 菜單可以利用生成的後驗分布進行假設檢驗。
在前面的例子中，可能會有讀者感到疑惑，當先驗分布選為共軛先驗時，後驗分布是有解析表達式的，這意味著其概率密度函數應該是確定的。但在結果顯示的界面中，明顯有模擬 (simulation)、迭代 (iteration) 和抽樣 (sampling) 這樣的常見於數值方法的字眼。這其實是因為 Stata 並非使用解析方法對後驗分布進行求解。為了進一步深入了解現代貝葉斯方法的操作流程，我們有必要對 MCMC (Markov chain Monte Carlo，馬爾科夫鏈蒙特卡洛) 方法進行介紹。
對於 MCMC 方法，維基百科定義為：馬爾可夫鏈蒙特卡洛方法包括一類用於從概率分布中採樣的算法。通過構建具有所需分布作為其平衡分布的馬爾可夫鏈，可以通過記錄鏈中的狀態來獲得所需分布的樣本。
這一段話有些抽象，實際上從應用層面上來說，我們只需要了解兩個最常用的 MCMC 方法：Metropolis–Hastings 算法 (簡稱 M-H 算法) 和吉布斯採樣 (Gibbs sampling)。前者是 Stata 內置的對後驗分布進行抽樣的算法，其主要作用是通過後驗分布的概率密度之比構造了一個具有接受-拒絕機制的馬爾科夫鏈，並通過它抽取出所需的後驗分布。這一算法也廣泛應用於 DSGE 模型的貝葉斯估計中，但它不是本文的重點。
在宏觀計量經濟學中最常用的 (主要是因為其應用於 VAR 模型的貝葉斯估計) MCMC 方法是吉布斯採樣。為了讓大家對吉布斯採樣的應用背景有一個簡單的了解，下面將以線性回歸為例進行介紹。
假設有如下的線性回歸模型：
其中， 是一個 的矩陣， 是一個 的矩陣， 是一個 的矩陣。其似然函數的表達式如下：
這裡有兩個待估參數，即係數矩陣 和擾動項方差 。分別考慮以下三種情況：
在第一種情況中，由於擾動項方差已知，此時的共軛先驗是正態分布，因此假定係數矩陣 的先驗分布為 。
先驗分布：
似然函數：
後驗分布：
可以發現後驗分布仍然是正態分布：
在第二種情況中，由於似然函數中的給定變量改變，其共軛先驗也隨之發生變化。查表知此時擾動項方差 應服從倒伽馬分布 (Inverse Gamma distribution)，具體運算過程和第一種情況類似，可以獲得相同分布族的後驗分布 。
在實際中，更加常見的是第三種情況，即兩個待估參數都未知。此時，沒有辦法獲得單一的共軛先驗，一種解決辦法是設定聯合先驗分布：
其中， 的分布仍是無條件的，和第二種情況中的設定一樣，但 則是有條件的，我們將其設定為：
在這樣的設定下，聯合先驗分布還是一個正態分布，和似然函數的分布族一致，這種情況稱之為自然共軛先驗 (natural conjugate prior)。
需要注意的是，最後解出的後驗分布是一個聯合後驗分布 ，但我們想考察的後驗分布則是邊緣後驗分布 和 。在自然共軛先驗的設定下，可以通過積分的方法求得邊緣後驗分布的解析形式。但這也對先驗分布的選擇施加了很強的限制。在較複雜的模型中，我們往往不希望限制先驗分布的選擇，但非自然共軛先驗難以通過積分求解出邊緣後驗分布，因此需要應用吉布斯採樣。
吉布斯採樣的思想其實很簡單，當有了兩個能條件分布的抽取器時 (即第一種和第二種情況)，可以通過迭代的方法獲得逼近邊緣分布的抽樣。在本例中，具體操作為：
這一操作步驟看似複雜，但在 Stata 中的實現其實很容易，只需要在線性回歸的命令前加上 bayes: 即可，但需要注意的是，此時回歸命令不能縮寫，並且如果要用吉布斯採樣，需要使用 bayes,gibbs:。
. clear //清空工作區. set seed 20 //設定隨機種子. sysuse auto //使用系統自帶的 auto 數據集. bayes, gibbs: regress price mpg weight //執行吉布斯採樣的貝葉斯估計Burn-in ...Simulation ...Model summary----Likelihood: price ~ normal(xb_price,{sigma2})Priors: {price:mpg weight _cons} ~ normal(0,10000) (1) {sigma2} ~ igamma(.01,.01)----(1) Parameters are elements of the linear form xb_price.Bayesian linear regression MCMC iterations = 12,500Gibbs sampling Burn-in = 2,500 MCMC sample size = 10,000 Number of obs = 74 Acceptance rate = 1 Efficiency: min = .9457 avg = .9864Log marginal likelihood = -696.85658 max = 1 ---- | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval]----+----price | mpg | -5.282435 27.88942 .278894 -5.410969 -59.53671 49.12283 weight | 2.074039 .198759 .001988 2.076242 1.682404 2.464258 _cons | 1.928181 99.76441 .997644 2.99603 -202.2644 194.9716----+---- sigma2 | 6448502 1099995 11311.5 6331241 4641478 8928296----Note: Default priors are used for model parameters. 6. 參考資料溫馨提示： 文中連結在微信中無法生效。請點擊底部「閱讀原文」。M-H 算法和共軛先驗的維基百科 -Link1- -Link2-Blake A P, Mumtaz H. Applied Bayesian econometrics for central bankers[J]. Technical Books, 2012. -Link-Koop G, Korobilis D. Bayesian multivariate time series methods for empirical macroeconomics[M]. Now Publishers Inc, 2010. -Link-Bill Rising. Bayesian Analysis using Stata. 2015. -Link- 7. 相關推文Note：產生如下推文列表的命令為： lianxh VAR DSGE 分布, m安裝最新版 lianxh 命令： ssc install lianxh, replace溫馨提示： 文中連結在微信中無法生效。請點擊底部「閱讀原文」。Stata程序：Monte-Carlo-模擬之產生符合特定分布的隨機數VaR 風險價值: Stata 及 Python 實現 🎦 2021 文本分析·爬蟲·機器學習🍓 課程主頁：https://gitee.com/lianxh/text ⏩ 2021 生存分析專題⛪ 課程主頁：https://gitee.com/lianxh/ST 關於我們🍎 連享會 ( 主頁：lianxh.cn ) 由中山大學連玉君老師團隊創辦，定期分享實證分析經驗。👉 直達連享會：【百度一下：連享會】即可直達連享會主頁。亦可進一步添加 主頁，知乎，面板數據，研究設計 等關鍵詞細化搜索。連享會主頁 lianxh.cn New！ lianxh 命令發布了： GIF 動圖介紹隨時搜索連享會推文、Stata 資源，安裝命令如下： . ssc install lianxh使用詳情參見幫助文件 (有驚喜)： . help lianxh
溫馨提示： 文中連結在微信中無法生效。請點擊底部「閱讀原文」。
Note：產生如下推文列表的命令為： lianxh VAR DSGE 分布, m安裝最新版 lianxh 命令： ssc install lianxh, replace
🎦 2021 文本分析·爬蟲·機器學習🍓 課程主頁：https://gitee.com/lianxh/text
⏩ 2021 生存分析專題⛪ 課程主頁：https://gitee.com/lianxh/ST
New！ lianxh 命令發布了： GIF 動圖介紹隨時搜索連享會推文、Stata 資源，安裝命令如下： . ssc install lianxh使用詳情參見幫助文件 (有驚喜)： . help lianxh