用于后验估计的调谐切片采样器

打开脚本

这个例子说明了如何改进切片采样器来进行贝叶斯线性回归模型的后验估计和推断。

在估计由数据似然和半共轭或自定义先验模型组成的后验模型时，估计使用MCMC采样器。如果样本的跟踪图显示短暂行为或非常高的序列相关性，或者您希望从后部存储少量样本，则可以指定老化样本或细化。

考虑预测美国实际国民生产总值的多元线性回归模型(国民生产总值)使用工业生产指数的线性组合(IPI公司)，总就业人数(E类)，以及实际工资(总重).

$$\texttt{GNPR}\t=\beta\u0+\beta\u1\texttt{IPI}\beta\u2\texttt{E}\beta\u3\texttt{WR}\varepsilon\t$$

为所有人美元$时间点，$\varepsilon_t美元$是一系列平均值为0且方差为0的独立高斯扰动$\sigma^2码$. 假设先前的分配是：

  • $\beta_j\vert\sigma^2版$是四维的t型各分量50自由度分布和相关矩阵的恒等式矩阵。而且，分布集中在${\left[{\begin{array}{{{20}{c}{-25}&}38；4&}38；0&}38；3\end{array}右]^\prime}$每个分量都由向量的相应元素缩放${\left[{\begin{array}{{{20}{c}{1}&；1&；1&；1\end{array}}\right]^\prime}$.

  • $\sigma^2\sim IG（3,1个）$.

贝耶斯姆将这些假设和数据可能性视为对应的后验是分析上难以处理的。

声明一个MATLAB®函数：

  • 接受的价值观$\beta版$和$\sigma^2码$在一个列向量中，以及超参数的值。

  • 返回联合先验分布的值，$\pi\左（\beta、\sigma ^2\右）$，给定$\beta版$和$\sigma^2码$.

功能logPDF=priorMVTIG（参数，ct，st，dof，C，a，b）%多元t次反γ的priorMVTIG测井密度%priorMVTIG将参数（1:end-1）传递给多元t密度%每个组件的自由度为自由度且为正的函数%定相关矩阵C.priorMVTIG返回%两个评估密度。%%params：计算密度的参数值%m乘1数值向量。%%ct：多元t分布成分中心，an（m-1）-by-1%数值向量。元素对应于第一个m-1元素%参数的。%%st：多元t分布分量表，an（m-1）-by-1％数值（M-1）×1数值向量。元素对应%params的第一个m-1元素。%%dof：多元t分布的自由度，a%数值标量或（m-1）-x-1数值向量。priorMVTIG扩张%标量，使dof=dof*one（m-1,1）。自由度要素%对应于params的元素（1:end-1）。%%C：多元t分布的相关矩阵%（m-1）-by-（m-1）对称正定矩阵。行和%列对应于params的元素（1:end-1）。%%a：反伽马形状参数，正数标量。%%反伽马标度参数，正标量。%β=params（1:（end-1））；sigma2=params（end）；tVal=（β-ct）。/st；mvtDensity=mvtpdf（tVal，C，dof）；igDensity=sigma2^（-a-1）*exp（-1/（sigma2*b））/（gamma（a）*b^a）；logPDF=log（mvtDensity*igDensity）；结束

创建一个匿名函数，其操作方式如下普里奥姆夫蒂格，但只接受参数值，并将超参数值固定在其值上。

dof=50；C=eye（4）；ct=zeros（4,1）；st=one（4,1）；a=3；b=1；logPDF=@（params）priorMVTIG（params，ct，st，dof，C，a，b）；

为线性回归参数创建自定义关节先验模型。指定预测数，第页. 另外，指定普里奥姆夫蒂格以及变量名。

p=3；PriorMdl=bayeslm（p，'模型类型','自定义','日志PDF'，日志pdf，...'变量名',[“IPI”“E”“总重”]);

普里奥姆德尔是一个客户贝叶斯线性回归模型对象表示回归系数和扰动方差的先验分布。

加载纳尔逊 - 普洛瑟数据集。创建响应和预测序列变量。

负载数据处理器X=数据表{：，PriorMdl.VarNames公司（2:结束）}；y=数据表{：，'GNPR'};

估计$\beta版$和$\sigma^2码$. 若要确定适当的老化期，请指定“无老化期”。指定切片采样器的宽度，该宽度接近假定漫反射先前模型的参数的后验标准差。绘制15000个样品。

宽度=[20,0.5,0.01,1,20]；numDraws=15e3；rng（1）%再现性PosteriorMdl1=估计值（PriorMdl，X，y，“伯宁”，0，'宽度'，宽度，...'裸体'，努德拉）；
警告：MCMC的高自相关得出：IPI，E，WR，西格玛-2。方法：MCMC 15000抽样得出的观测总数62号码预测的：4 |均值标准CI95正分布-----------------------------------------------------------------------拦截|-0.6669 1.5325 [-3.202，1.621] 0.309实证IPI |4.6196 0.1021 [4.437，4.839] 1.000实证E |0.0004 0.0002 [0.000，0.001] 0.995实证WR |2.6135 0.3025 [1.982，3.152] 1.000实证西格玛-2 |48.9990 9.1926 [33.147，67.993] 1.000实证

后周期DL1是一个经验主义模型对象。这个贝塔德劳斯和信号图属性存储切片采样器中的样本。

绘制马尔可夫链的迹图和自相关图。

数字；对于j=1:4子块（2,3，j）图（后周期dl1.BetaDraws（j，：））；轴紧的xlabel公司('示例索引')标题（sprintf(['跟踪图'字符（8212）'%s'],PriorMdl.VarNames公司{j} ）；结束子块（2,3,5）绘图（后周期性1.Sigma2Draws）；轴紧的xlabel公司('示例索引')标题(['跟踪图'字符（8212）'信号2'])；数字；对于j=1:4子块（2,3，j）自相关（后周期dl1.BetaDraws（j，：））；轴紧的标题（sprintf([“相关图”字符（8212）'%s'],PriorMdl.VarNames公司{j} ）；结束子块（2,3,5）自相关（后周期性1.Sigma2Draws）；轴紧的标题([“相关图”字符（8212）'信号2']);

除截距外，所有参数均表现出显著的瞬态效应和高自相关。估计检查高相关性，如果存在，则发出警告。而且，这些参数似乎没有收敛到它们的平稳分布。

再次估计后验分布。为了减少短暂的影响，使用默认的磨合期(5000个，这是默认设置）。为了减少串行相关性的影响，将其减薄20倍并指定15e4/薄有效抽签次数。此外，数值稳定估计$\sigma^2码$通过指定重新参数化为$\log（\sigma^2）$. 绘图跟踪和自相关绘图。

thin=20；numDraws=15e4/thin；PosteriorMdl2=估计值（PriorMdl，X，y，“薄”，薄，'宽度'，宽度，...'重新参数化'，对，'裸体'，numDraws）；数字；对于j=1:4子块（2,3，j）图（后周期dl2.BetaDraws（j，：））；轴紧的xlabel公司('示例索引')标题（sprintf(['跟踪图'字符（8212）'%s'],PriorMdl.VarNames公司{j} ）；结束子块（2,3,5）图（后周期性2.Sigma2Draws）；轴紧的xlabel公司('示例索引')标题(['跟踪图'字符（8212）'信号2'])；数字；对于j=1:4子块（2,3，j）自动更正（PosteriorMdl2.BetaDraws（j，：））；标题（sprintf([“相关图”字符（8212）'%s'],PriorMdl.VarNames公司{j} ）；结束子批次（2,3,5）自动相关（后周期性2.Sigma2Draws）；标题([“相关图”字符（8212）'信号2']);
方法：MCMC抽样7500份，抽取观察数据62份，预测因子：4 |平均Std-CI95正态分布，即截距|-0.4777 1.2225[-2.911，1.955]0.338经验IPI | 4.6506 0.1128[4.434，4.884]1.000经验E | 0.0005 0.0002[0.000，0.001]0.991经验WR | 2.5243 0.3468[1.815，3.208]1.000经验Sigma2 | 51.3444 9.3170[36.319，72.469]1.000经验

跟踪图表明，马尔可夫链：

  • 已达到其平稳分布，因为样本的分布具有相对恒定的均值和方差。具体来说，点的分布似乎并没有改变它。

  • 混合相当快，因为对于某些参数来说，序列相关性是中等到非常低的。

  • 不要表现出任何短暂的行为。

这些结果表明，你可以进行敏感性分析，或者进行后验推断。

另见

功能

物体

相关话题

计量经济学工具箱文档
金宝app
用MATLAB建立财务风险模型的实用指南

用MATLAB建立财务风险模型的实用指南

下载电子书