引言

马尔科夫链蒙特卡洛模拟（MCMC）

计量模型的贝叶斯估计

特定模型的似然函数或后验分布

贝叶斯预测

贝叶斯方法提供了与传统频数方法不同的视角，参数估计、模型选择等。

贝叶斯估计通过蒙特卡洛模拟获得参数的抽样分布，更适合于复杂的模型。

贝叶斯推断的优势：将多种不确定性，包括数据、模型和参数，融合到统一的框架内。

引言

马尔科夫链蒙特卡洛模拟（MCMC）

计量模型的贝叶斯估计

特定模型的似然函数或后验分布

贝叶斯预测

随机过程\(X_t\)取不同数值\(S=(1,2,\cdots, s)\)，转移概率为 \[ p_{ij} = P(X_{t+1}=j|X_t=i), i,j \in S. \]

比如， \[ P=\begin{bmatrix} p_{11} & p_{12} \\ p_{21} & p_{22} \end{bmatrix} = \begin{bmatrix} 1-\alpha & \alpha \\ \beta & 1-\beta \end{bmatrix}. \]

马尔可夫特征：随机变量的\(t\)期状态只取决于\(t-1\)期状态。即下一步去哪里只取决于变量当前所在的位置，而不取决于是怎么到达当前位置的。

马尔科夫链：满足马尔科夫特征的随机过程。

如果\(\pi_t=(\pi_1, \cdots, \pi_s)\)，那么\(\pi_{t+1} = \pi_t P\). \[ \begin{bmatrix} \pi_{1,t+1} & \pi_{2,t+1} \end{bmatrix} = \begin{bmatrix} \pi_{1t} & \pi_{2t} \end{bmatrix} \begin{bmatrix} p_{11} & p_{12} \\ p_{21} & p_{22} \end{bmatrix} \]

\(\pi_{t+n} = \pi_{t+n-1} P=\pi_{t+n-2} P^2=\pi_{t} P^n\).

平稳分布：如果分布\(\pi_t\)不再变化，即\(\pi = \pi P\)，\(\pi\)称为平稳分布。

\(\pi (I-P) = 0\)，所以\(\pi\)为P矩阵特征值1对应的特征向量。

如果转移矩阵为 \[ \begin{bmatrix} 0.75 & 0.25 \\ 0.125 & 0.875 \end{bmatrix} \] 可以求出其对应的平稳分布为(1/3, 2/3).

对于一个转移矩阵，平稳分布是不是唯一的？ \[ \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \]

无穷多平稳分布：转移矩阵为 \[ \begin{bmatrix} 1/3 & 2/3 & 0 \\ 1/3 & 2/3 & 0 \\ 0 & 0 & 1 \end{bmatrix} \]

令\(P^{(n)} = \pi_{t} P^n\)，如果\(P^{(n)}_{ij}>0\)，称\(j\)从\(i\)可达（accessible）。

如果\(P^{(n)}_{ij}>0\)，\(P^{(n)}_{ji}>0\)，称\(i\)与\(j\)互通(communicate)。一组互通的状态构成互通类(communicate class)。

只包含一个互通类的马尔科夫链称为不可约的(irreducible)。

定理：不可约的马尔科夫链具有唯一的平稳分布。

含义：对于目标分布，如果能找到其对应的马尔科夫链，那么总可以到达该分布。

如何找到马尔科夫链：MCMC

MCMC: Gibbs抽样(Geman and Geman, 1984), Metropolis-Hastings抽样(Metropolis et al., 1953; Hastings, 1970).

MH抽样；设目标分布为\(f(x)\)，

给定\(x_t\), 定义工具分布\(q\), 从\(q(y|x_t\))\(抽取随机数\)y_t$

定义接受概率(acceptance probability): \(\rho(x_t, y_t)=\min \left[ \frac{f(y_t)}{f(x_t)} \frac{q(x_t|y_t)}{q(y_t|x_t)}, 1 \right]\), \[ x_{t+1} = \begin{cases} y_{t} & \rho(x_t, y_t) \\ x_{t} & 1- \rho(x_t, y_t) \\ \end{cases} \]

对称情形：\(\rho(x_t, y_t)=\min \left[ \frac{f(y_t)}{f(x_t)}, 1 \right]\),

独立情形：\(\rho(x_t, y_t)=\min \left[ \frac{f(y_t)}{f(x_t)} \frac{q(x_t)}{q(y_t)}, 1 \right]\),

target: Beta(2.7, 6.3); proposal: uniform(0,1).

target: Cauchy; proposal: t(1).

将随机变量分为两组，\((x_t, y_t)\)，联合密度为\(f(x,y)\)。

抽取\(x_0\)

\(t=1,2,\cdots\)，抽取\(y_{t+1} \sim f(y|x_t)\), \(x_{t+1} \sim f(x|y_{t+1})\)

例：二元正态分布。

抽取\(x_0=0\)

\(t=1,2,\cdots\)，抽取 \[ \begin{eqnarray*} y_{t+1} | x_t & \sim & N \left( \frac{\sigma_1}{\sigma_2} \rho x_t,(1-\rho^2) \sigma_1^2 \right) \\ x_{t+1} | y_{t+1} & \sim & N \left( \frac{\sigma_2}{\sigma_1} \rho y_{t+1},(1-\rho^2) \sigma_2^2 \right) \\ \end{eqnarray*} \]

可以直接扩展到随机变量分为多组的情况。

对协方差矩阵的抽样一般采用Gibbs抽样，inverse Wishart为最广泛采用的先验分布。

设\(T_0\)为burn-in样本量，\(T\)为MCMC的样本量，整个MCMC迭代的次数为\(T_{total}=T_0+(T-1)\times thinning + 1\)。

MH抽样生成试验值的公式为: \(\theta_{*} = \theta_{t-1} + e, e \sim N(0, \rho_k^2 \Sigma_k)\)，所谓适应性MH抽样是指每隔一定的区间调整\(\rho_k\)和\(\Sigma_k\)，以实现最优的接受率(TAR)。Gelman, Gilks and Roberts (1997)证明，单一参数的TAR为0.44，多个参数的TAR为0.234。

设Alen为每个适应性区间的长度（adaptation(every())）。可以设置最小迭代次数（adaptation(miniter())）和最大迭代次数（adaptation(maxiter())）。

初始化：\(\theta_t = \theta_0^{(f)}\)，\(\rho_0=2.38/\sqrt{d}\)，其中d为参数个数。\(\Sigma_0 = I\)或者由covariance()选项自行设定。

生成试验值: \(\theta_{*} = \theta_{t-1} + e, e \sim N(0, \rho_k^2 \Sigma_k)\)。

接受概率为 \[ \alpha(\theta_{*}|\theta_{t-1}) = \text{min} \left{ \frac{p(\theta_{*}|y)}{p(\theta_{t-1}|y)}, 1 \right}. \]

其中，\(p(\theta_{*}|y) = f(y|\theta) p(\theta)\)为\(\theta\)的后验分布。

\(\rho_k\)和\(\Sigma_k\)的更新公式为： \[ \rho_k = \rho_{k-1} \exp( \beta_k [\Phi^{-1}(AR_k/2) - \Phi^{-1}(TAR/2) ] ) \]

其中， \[ AR_k = (1-\alpha) AR_{k-1} + \alpha AP \]

AP为接受概率\(\alpha(\theta_{*}|\theta_{t-1})\)的均值。\(\alpha\)的默认值为0.75，或者由adaptation(alpha())自行设定。

\[ \Sigma_k = (1-\beta_k) \Sigma_{k-1} + \beta_k \hat{\Sigma}_k, \beta_k = \beta_0/k^{\gamma}. \]

其中，\(\Sigma_0 = I\)或者由covariance()选项自行设定；\(\hat{\Sigma}_k\)为MCMC样本的协方差。\(\beta_0\)的默认值为0.8，或者由adaptation(beta())自行设定。\(\gamma\)的默认值为0,或者由adaptation(gamma())自行设定。

引言

马尔科夫链蒙特卡洛模拟（MCMC）

计量模型的贝叶斯估计与检验

特定模型的似然函数或后验分布

贝叶斯预测

\[ \pi(\theta|y) = \frac{f(y|\theta) \pi(\theta)}{f(y)} \propto f(y|\theta) \pi(\theta). \]

其中，\(f(y|\theta) \pi(\theta)\)叫做贝叶斯核(kernel)。

如果先验分布\(\pi(\theta)\)对后验分布的影响非常小， 称其是无信息的(noninformative, vague, diffuse, flat)。

如果\(\pi(\theta)\)的积分不等于1，称其是不恰当的(improper)。

Jeffreys先验：\(\pi(\theta) \propto |I(\theta)|^{1/2}\)，其中，\(I(\theta)\)为Fisher信息矩阵， \[ I(\theta) = -E\left[ \frac{\partial^2 \log(p(y|\theta))}{\partial \theta^2} \right]. \]

模型: \(y=x\beta + u\)

似然函数: \(y_i \sim Normal(x_i \beta, \sigma^2)\)

先验: \[ \begin{eqnarray*} \beta & \sim & \text{Mumltivariate-Normal}(\beta_0, V_0) \\ \sigma^2 & \sim & \text{Inverse-Gamma}(\alpha_0/2, \delta_0/2) \\ \end{eqnarray*} \]

一般的逆Gamma分布，InverseGamma[\(\alpha, \beta,\gamma,\mu\)] 定义域为(\(\mu, \infty\))，\(\mu\)为位置参数，\(\alpha, \gamma\)为形状参数，\(\beta\)为刻度参数。

InverseGamma[\(\alpha, \beta\)]表示InverseGamma[\(\alpha, \beta, 1, 0\)], pdf为 \[ pdf(x) = \frac{\exp(-\beta/x) (\beta/x)^{\alpha}}{x \Gamma(\alpha)} \]

\[ E(x) = \frac{\beta}{\alpha-1} (\text{ for } \alpha>1); Var(x) = \frac{\beta^2}{(\alpha-2)(\alpha-1)^2} (\text{ for } \alpha>2). \]

如果\(x\sim Gamma(\alpha, \beta)\)，那么\(1/x \sim Inv-Gamma(\alpha, 1/beta)\).

边际后验分布： \[ \beta | \sigma^2,y \sim \text{Mumltivariate-Normal} \left( \beta_1, V_1 \right) \]

其中， \[ V_1 = \left[ \sigma^{-2} X'X + V_0^{-1} \right]^{-1}, \beta_1 = V_1 \left[ \sigma^{-2} X'y + V_0^{-1} \beta_0 \right] \]

\[ \sigma^2 | \beta,y \sim \text{Inverse-Gamma}\left( (\alpha_0+n)/2, [\delta_0+(y-x\beta)'(y-x\beta)]/2 \right) . \]

\(\beta\) are sampled in one block (Greenberg, 2013).

设定初始值\(\sigma^{2(0)}\).

第\(t\)次迭代, 从下面的分布抽样 \[ \begin{eqnarray*} \beta^{(t)} | \sigma^2,y & \sim & \text{MNormal} \left( \beta_1, V_1 \right) \\ \sigma^2 | \beta,y & \sim & \text{IG}\left( (\alpha_0+n)/2, [\delta_0+(y-x\beta^{(t)})'(y-x\beta^{(t)})]/2 \right) . \end{eqnarray*} \]

bayes [, bayesopts] : estimation_command [, estopts]

其中，priorspec的形式为：{[eqname:]param[, matrix]}, priordist。比如，

引言

马尔科夫链蒙特卡洛模拟（MCMC）

计量模型的贝叶斯估计与检验

特定模型的似然函数或后验分布

贝叶斯预测

收敛性(是否达到平稳分布)

踪迹图（围绕常数均值波动）

核密度图（模拟样本分为前后两段，比较其分布，或者做等均值等方差检验）

自相关

Gelman-Rubin (shrink factor)统计量：多条马尔科夫链进行抽样，那么链之间的平均差异与链内的平均差异应该进行相等。二者的比例称作收缩因子（shrink factor），作为经验规则， 收缩因子应低于1.1。

Efficiency: 大于0.1，视作良好。低于0.01，严重关注。

bayesgraph type which

其中，type包括:diag, trace, histogram, kdensity, ac, cusum, matrix.

which包括： _all, {[*eqname*:]*param*}

所有参数：

所有参数：

单个参数：

bayesstats ic: 三个统计指标。

边际似然函数

dic (应选择dic最低的模型)

bf (贝叶斯因子=两个模型的边际似然函数的比（相同的样本）)

贝叶斯因子 (Jeffreys, 1961)： \[ BF_{12} = \frac{P(y|M_1)}{P(y|M_2)} = \frac{m_1(y)}{m_2(y)}. \]

posterior odds: \[ \frac{P(M_1|y)}{P(M_2|y)} = \frac{P(y|M_1)P(M_1)}{P(y|M_2)P(M_2)} = BF_{12} \frac{P(M_1)}{P(M_2)}. \]

Kass and Raftery (1995): 贝叶斯因子：

可以用于比较不同的先验、不同的后验分布、不同的回归函数、不同的解释变量等。

条件：不同模型使用完全相同的样本；具有恰当的后验分布；MC收敛。

对参数落在某个区间的检验：

Stata对不同模型的不同类型的参数设置了不同的先验分布。

回归系数的先验分布默认为独立的正态分布\(N(0, \sigma_{prior}^2)\)，其中\(\sigma_{prior}^2=10,000\)。

正值的参数（比如方差）默认先验分布为Inverse-Gamma(\(\alpha, \beta\))。默认值为\(\alpha=0.01\), \(\beta=0.01\)。

排序选择模型的切点：先验分布为Unifrom(0,1)。

多水平模型：随机效应协方差矩阵

independent, identity: Inverse-Gamma(\(\alpha, \beta\))

unstructured: Invser-Wishart(d+1, I(d))，其中d为协方差矩阵的维度，d+1为自由度。

引言

马尔科夫链蒙特卡洛模拟（MCMC）

计量模型的贝叶斯估计与检验

特定模型的似然函数或后验分布

贝叶斯预测

bayesmh可以用于估计单方程线性模型、单方程非线性模型、多方程线性模型、多方程非线性模型、多水平模型、概率分布等。

bayesmh depvar [indepvarspec] [if] [in] [weight], likelihood(modelspec) prior(priorspec) llevaluator(log-likelihood) evaluator(log-postdist) [options]

三种用法：

likelihood() + prior(): 内置似然函数 + 内置先验分布

llevaluator() + prior(): 自定义的对数似然函数 + 内置先验分布

evaluator(): 自定义的对数后验分布

似然函数即概率密度函数，不同分布涉及到的参数也不同。比如，正态分布有均值和方差两个参数；泊松分布只有一个参数，既为均值也是方差。

bayesmh中，分布的均值参数通过模型来设定。分布的形式和其它参数则通过likelihood()来设定。比如，

bayesmh y x1 x2, likelihood(normal({var}))，表示\(y\sim N(b_0+b_1 x_1 + b_2 x_2, var)\)

bayesmh y x1 x2, likelihood(poisson)，表示\(y\sim Poisson(exp(b_0+b_1 x_1 + b_2 x_2))\)

bayesmh y x1 x2, likelihood(probit)，表示\(y\sim Bernoulli(\Phi(b_0+b_1 x_1 + b_2 x_2))\)

Logit模型的对数似然函数。

似然函数程序文件: normalthresh.ado

引言

马尔科夫链蒙特卡洛模拟（MCMC）

计量模型的贝叶斯估计与检验

特定模型的似然函数或后验分布

贝叶斯预测

参数是随机的，预测时从参数的后验分布随机抽样，得到参数的R个随机数，进而得到对每个观测值的预测，即每个观测值有R个预测（模拟）数值。

bayespredict {_ysim#[numlist]}

得到模拟的因变量的名称为：_ysim1_1, _ysim1_2, …，表示对第1个方程中第1、2、…个预测值的预测。

numlist：设定对哪些观测值进行预测。

模拟的前5个观测值的统计指标。

预测值的其它统计函数：

bayesgraph histogram [label:]@*func(arg1* [, arg2)], saving(file)

bayesgraph histogram [label:]@*userprog, arg1* [arg2], saving(file)

其中，func为mata函数（该函数对列向量操作得到标量），比如mean, variance等。

arg1和arg2为{_ysim[#]}, {_resid[#]}（残差）, 或{_mu[#]}（因变量的期望）。 预测值的函数（比如均值、方差、偏度等）：

可以直接得到R次预测数值的均值（或中位数、标准差）等统计量。

bayespredict newvar, mean median std cri outcome(depvar)

贝叶斯预测：对每个观测值，随机生成对应的因变量的S个随机数，计算其均值。

直接预测用户自定义的其它指标：

bayespredict [label:]@*func(arg1* [, arg2)], saving(file)

bayespredict [label:]@*userprog, arg1* [arg2], saving(file)

计算统计量的概率值。统计量可以利用mata函数来计算。

这些统计量的概率值越接近0.5，表明模拟值与实际值越吻合。

可以通过mata函数或者ado程序（或者二者混合使用）定义自己的统计量。

例：偏度和峰度。

预测时所模拟生成的前3个变量为：

贝叶斯预测也可以用来评估模型的拟合优度。利用后验分布随机模拟因变量的取值，并与实际值进行比较。