可以从许多统计软件包中运行Stan。到目前为止，我一直在从R运行Stan。

由Kaizong Ye，Weilong Zhang撰写

第一步是为Stan模型编写文件。

这包含一个文件linreg.stan：

可下载资源

 data {
  int N;
  vector[N] x;
  vector[N] y;
}

model {
  y ~ normal(alpha + beta * x, sigma);
}

该文件的第一部分称为数据，它声明了将作为输入传递给Stan的标量，向量和矩阵。

利用贝叶斯定理求出：

线性回归模型中的权重 $\boldsymbol{w}$ 的后验概率 $p(\boldsymbol{w}|Data)$ 的概率分布。
新数据点的预测值 $f(\boldsymbol{x}^*)|Data, \boldsymbol{x}^*$ 的概率分布。

线性回归模型

$\begin{aligned} & f(\boldsymbol{x})=\boldsymbol{x}^T \boldsymbol{w} \\ & y = f(\boldsymbol{x}) + \epsilon, \epsilon \sim N(0, \sigma^2) \end{aligned} \\$

其中 $\epsilon$ 为噪声

理解线性回归

对于线性回归这个问题，可以分别从频率派和贝叶斯派的观点来理解它。

在频率派的观点中，权值 $\boldsymbol{w}$ 是一个未知的常数，因此将问题转化为最优化问题，并对权值进行点估计。做点估计的方法又分为两种：

最大似然估计(Maximum Likelihood Estimation, MLE):

$\boldsymbol{w}_{MLE}=\mathop{\arg\min}_\boldsymbol{w} p(Data|\boldsymbol{w}) \\$

该问题即等价于最小二乘估计(Least Squares Estimator)

最大后验估计(Maximum A Posterior Estimation, MAP):

$\boldsymbol{w}_{MAP}=\mathop{\arg\min}_\boldsymbol{w} p(\boldsymbol{w}|Data) \\$

该问题即等价于正则化(Regularized)的最小二乘估计

而在贝叶斯派的观点中，权值 $\boldsymbol{w}$ 是一个随机变量，因此求的是该随机变量的条件分布： $p(\boldsymbol{w}|Data)$ 。本文的主题——贝叶斯线性回归模型就是沿用贝叶斯派的思路进行分析的。

条件分布 $p(\boldsymbol{w}|Data)$

根据贝叶斯定理 $p(\theta|x) = \dfrac{p(\theta)p(x|\theta)}{p(x)}$ ，后验概率可以由先验概率和似然概率共同求出。

首先，为方便起见，设先验概率(prior)为 $\boldsymbol{w} \sim N(0, \Sigma_p)$ 。

$\begin{aligned} & p(\boldsymbol{w}|Data) = p(\boldsymbol{w}|X,Y) = \dfrac{p(\boldsymbol{w},Y|X)}{p(Y|X)} = \dfrac{p(Y|\boldsymbol{w},X) p(\boldsymbol{w}|X)}{\int p(Y|\boldsymbol{w},X) p(\boldsymbol{w}|X)\mathrm{d}\boldsymbol{w}} \\ \end{aligned} \\$

由于先验概率 $\boldsymbol{w}$ 与训练集 $Data$ 独立，上式可写为

$p(\boldsymbol{w}|Data) = \dfrac{p(Y|\boldsymbol{w},X) p(\boldsymbol{w})}{\int p(Y|\boldsymbol{w},X) p(\boldsymbol{w})\mathrm{d}\boldsymbol{w}} \\$

由于 $N$ 个数据相互独立，可得

$\begin{aligned} p(Y|\boldsymbol{w},X) & = \Pi_i p(y_i|\boldsymbol{w}, \boldsymbol{x}_i) = \Pi_i N(y_i|\boldsymbol{w}^T\boldsymbol{x}_i, \sigma^2) \\ & = C \cdot exp(\dfrac{1}{-2\sigma^2}\sum_i (y_i-\boldsymbol{w}^T \boldsymbol{x}_i)^2) = C \cdot exp(\dfrac{1}{-2\sigma^2}(Y-XW)^T(Y-XW)) \\ \end{aligned} \\$

根据高斯分布的共轭性质， $\boldsymbol{w}|Data$ 也服从高斯分布，记为 $\boldsymbol{w}|Data \sim N(\mu_\boldsymbol{w}, \Sigma^2_\boldsymbol{w})$ 。将 $p(Y|\boldsymbol{w},X)$ 和 $p(\boldsymbol{w})$ 代入贝叶斯公式，可得

$\begin{aligned} p(\boldsymbol{w}|Data) & \propto exp(\dfrac{1}{-2\sigma^2}(Y-XW)^T(Y-XW) - \dfrac{1}{2} W^T \Sigma_p W) \\ & \propto exp(\dfrac{1}{-2\sigma^2}(Y^TY-\underbrace{2Y^TXW}_{w的一次项}+\underbrace{W^TX^TXW}_{w的二次项}) - \underbrace{\dfrac{1}{2}W^T\Sigma_pW}_{w的二次项}) \end{aligned} \\$

辅助结论

对于任意一个服从高斯分布的随机变量 $x\sim N(\mu, \Sigma)$ ，

$\begin{aligned} p(x) & \propto exp(-\dfrac{1}{2} (X-\mu)^T\Sigma^{-1}(X-\mu)) \\ & \propto exp(-\dfrac{1}{2} (\underbrace{X^T\Sigma^{-1}X}_{x的二次项} - \underbrace{2\mu^T\Sigma^{-1}X}_{x的一次项} + C)) \end{aligned} \\$

可见：二次项中含有方差的形式；一次项中既包含均值，又包含方差。由此，经过对比相应的项，即可从一次项和二次项中求出均值和方差。

将以上结论应用于 $p(\boldsymbol{w}|Data)$ ，可以求出 $\mu_\boldsymbol{w}$ 和 $\Sigma^2_\boldsymbol{w}$ ：

$\left\{ \begin{aligned} & \mu_\boldsymbol{w} = \sigma^{-2}A^{-1}X^TY \\ & \Sigma_\boldsymbol{w} = A^{-1} \\ & A = \sigma^{-2}X^TX + \Sigma_p^{-1} \end{aligned} \right. \\$

预测

已知 $\boldsymbol{w}|Data \sim N(\mu_\boldsymbol{w}, \Sigma^2_\boldsymbol{w})$ 以及 $f(\boldsymbol{x}^*) = (\boldsymbol{x}^*)^T\boldsymbol{w}$ ，

根据线性回归的模型以及高斯分布的性质，可得

$f(\boldsymbol{x}^*)|Data, \boldsymbol{x}^* \sim N((\boldsymbol{x}^*)^T \mu_\boldsymbol{w}, (\boldsymbol{x}^*)^T \Sigma_{\boldsymbol{w}}\boldsymbol{x}^*) \\$

最后加上噪声的影响，可得

$y^*|Data, \boldsymbol{x}^* \sim N((\boldsymbol{x}^*)^T \mu_\boldsymbol{w}, (\boldsymbol{x}^*)^T \Sigma_{\boldsymbol{w}}\boldsymbol{x}^*+\sigma^2) \\$

接下来，我们可以通过运行以下R代码来模拟数据集，并使用Stan和我们的文件linreg.stan来拟合模型：

stan(file = 'linreg. ', data = mydata, iter = 1000,   = 4)

第一次安装Stan模型时，模型编译成C ++时会有几秒钟的延迟。然而，一旦编译了模型，就可以将其应用于新的数据集而无需重复编译过程（执行模拟研究具有很大的优势）。

视频

R语言中RStan贝叶斯层次模型分析示例

在上面的代码中，我们要求Stan运行4个独立的链，每个链有1000次迭代。运行后，我们可以通过以下方式汇总输出：

Inference for Stan model: linreg.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
alpha  -0.10    0.00 0.10  -0.29  -0.16  -0.10  -0.04   0.09  1346    1
beta    0.95    0.00 0.11   0.75   0.88   0.95   1.02   1.17  1467    1
sigma   0.98    0.00 0.07   0.85   0.93   0.98   1.03   1.12  1265    1
lp__  -47.54    0.06 1.24 -50.77 -48.02 -47.24 -46.68 -46.17   503    1

Samples were drawn using NUTS(diag_e) at Mon Jun 08 18:35:58 2015.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

对于回归斜率β，我们的后验均值为0.95（接近用于模拟数据的真实值1）。为了形成95％的后验置信区间，我们简单地采用取样后验的2.5％和97.5％的百分位数，这里是0.75到1.17。

您可以从拟合的模型中获取各种其他数量。一种是绘制其中一个模型参数的后验分布。要获得回归斜率，我们可以执行以下操作：

 hist(result$beta)

β后验分布直方图

最受欢迎的见解

1.R语言基于ARMA-GARCH-VaR模型拟合和预测实证研究

2.R语言时变参数VAR随机模型

3.R语言估计时变VAR模型时间序列的实证研究

4.R语言基于ARMA-GARCH过程的VAR拟合和预测

5.GARCH（1,1），MA以及历史模拟法的VaR比较

6.R语言用向量自回归（VAR）进行经济数据脉冲响应

7.R语言实现向量自动回归VAR模型

8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型

9.R语言VAR模型的不同类型的脉冲响应分析

现在让我们使用标准普通最小二乘拟合线性模型：

Residuals:
    Min      1Q  Median      3Q     Max 
-1.9073 -0.6835 -0.0875  0.5806  3.2904 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.10280    0.09755  -1.054    0.295    
x            0.94753    0.10688   8.865  3.5e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9707 on 98 degrees of freedom
Multiple R-squared:  0.4451,    Adjusted R-squared:  0.4394 
F-statistic:  78.6 on 1 and 98 DF,  p-value: 3.497e-14

这给出了我们对斜率0.95的估计，与Stan的后验平均值相差2位小数，标准误差为0.11，这与Stan的后验SD相同。