基于哈密顿蒙特卡罗的贝叶斯线性回归

打开脚本

这个例子展示了如何使用哈密顿蒙特卡罗(HMC)采样器对线性回归模型进行贝叶斯推断。

在贝叶斯参数推理中，目标是通过纳入模型参数的先验知识来分析统计模型。自由参数的后部分布 $\ theta$ 结合似然函数 $P (y | \θ)美元$ 先验分布 $P(\θ)美元$ ，使用贝叶斯定理：

$p（\ theta | y）= \ frac {p（y | \ theta）p（\ theta）} {p（y）}。$

通常，总结后验分布的最好方法是用蒙特卡罗方法从该分布中获得样本。使用这些样本，您可以估计边缘后验分布和派生的统计数据，如后验均值、中位数和标准偏差。HMC是一种基于梯度的马尔可夫链蒙特卡罗采样器，它比标准采样器更高效，特别是对于中维和高维问题。

线性回归模型

用截距分析线性回归模型 $\α美元$ ，线性系数 $β\美元$ (一列向量)，和噪声方差 $\σ^ 2美元$ 将数据分布作为自由参数。假设每个数据点都有独立的高斯分布:

$$ $ y_i | \θ\ sim \ mathcal {N} (\ mu_iθ(\)\σ^ 2)。＄＄$

模型的平均 $\ mu_i美元$ 高斯分布的函数 x_i美元和模型参数为

$$ mu_i = + x_i^T， $$

在贝叶斯分析中，您还必须为所有自由参数分配先验分布。在截距和线性系数上分配独立的高斯先验:

$\ alpha \ sim \ mathcal {n}（\ alpha_0，\ sigma_ \ alpha ^ 2）$ ，

$\ beta_i \ sim \ mathcal {n}（\ beta__ {0}，\ sigma _ {\ beta} ^ 2）$ ．

为了使用HMC，所有的抽样变量必须是无约束的，这意味着后验密度及其梯度必须对所有的实参数值有良好的定义。如果你有一个参数被限制在一个区间内，那么你必须把这个参数转换成一个无界的。为了保留概率，你必须将先验分布乘以相应的雅可比因子。同时，在计算后验梯度时也要考虑这个因素。

噪声方差是（平方）刻度参数，只能是正的。然后，它可以更容易，更自然地将其对数视为无界限的自由参数。在噪声方差对数之前分配正常情况：

$$ \ log \σ^ 2 \ sim \ mathcal {N} (\ kappa \ω^ 2)美元$ ．

写出自由参数的后验密度的对数 $\ theta =（\ alpha; \ beta; \ log \ sigma ^ 2）$ 作为

$logp (\theta | y) = mathrm{const。} + log P(y | \theta) + log P(\theta).$

忽略常量术语并调用最后两种术语的总和 $g（\ theta）$ ．要使用HMC，请创建一个计算函数句柄 $g（\ theta）$ 及其梯度 ${\部分g} /{\部分\θ}$ 对于任何价值 $\ theta$ ．用于计算的函数 $g（\ theta）$ 位于脚本的末尾。

创建数据集

定义截距的真实参数值，即线性系数β以及噪声标准偏差。知道真正的参数值使得可以与HMC采样器的输出进行比较。只有第一个预测器影响响应。

NumPredictors = 2;trueIntercept = 2;trueBeta = (3; 0);trueNoiseSigma = 1;

使用这些参数值在两个预测器的随机值上创建一个正态分布的样本数据集。

NumData = 100;rng ('默认'）重复性的％x = rand（numdata，numpredictors）;mu = x * truebeta + trueintercept;y = normrnd（mu，truenoiseigma）;

定义后验概率密度

选择高斯先验的均值和标准差。

拦截priormean = 0;拦截产品= 10;betapriormean = 0;Betapriorigma = 10;lognoisevariAcemean = 0;lognoisevariaresigma = 2;

保存一个函数logPosterior在MATLAB®路径上，返回产品的乘积对数和可能性的对数，以及该对数的梯度。的logPosterior函数在此示例的末尾定义。然后，调用参数的函数来定义logpdf.的输入参数HMCSampler.函数。

logpdf = @(参数)logPosterior(参数,X, y,．..InterceptPriorMean InterceptPriorSigma,．..Betapriormean，Betapriorigma，．..LogNoiseVarianceMean LogNoiseVarianceSigma);

创建HMC采样器

定义开始抽样的初始点，然后调用HMCSampler.函数来创建哈密顿采样器HamiltonianSampler对象。显示采样器属性。

拦截= randn;β= randn (NumPredictors, 1);LogNoiseVariance = randn;曾经繁荣=[拦截,β;LogNoiseVariance];smp = hmcSampler (logpdf,曾经繁荣,“NumSteps”，50）;

SMP.

smp = HamiltonianSampler with properties: StepSize: 0.1000 NumSteps: 50 MassVector: [4x1 double] JitterMethod: 'jitter-both' StepSizeTuningMethod: 'dual-averaging' MassVectorTuningMethod: ' iter- sampling' LogPDF: [function_handle] VariableNames: {4x1 cell} StartPoint: [4x1 double]

估计映射点

估计后验密度的MAP(最大后验概率)点。您可以从任何点开始采样，但通常更有效的方法是估计MAP点，然后将其作为调优采样器和绘制样本的起点。估计和显示MAP点。属性可以在优化过程中显示更多信息“VerbosityLevel”值为1。

[mappars，fitinfo] = viettatemap（SMP，“VerbosityLevel”，0）;MapIntercept = Mappars（1）Mapbeta = Mappars（2：end-1）Maplognoisevariance = mappars（END）

MAPIntercept = 2.3857 MAPBeta = 2.5495 -0.4508 mapplognoise variance = -0.1007

为了检查优化是否收敛到局部最优，绘制fitInfo。客观的字段。这个字段包含函数优化的每次迭代的负对数密度的值。最终的值都是相似的，所以优化已经收敛。

情节(fitInfo.Iteration fitInfo.Objective,'ro-'）;包含(“迭代”）;ylabel (“负对数密度”）;

调整取样器

为获得有效的采样，选择合适的采样参数值是非常重要的。找到好的值的最好方法是自动调优MassVector.，一步的大小，和numsteps.参数使用tuneSampler方法。使用方法：

1.调优MassVector.取样器。

2.调优一步的大小和numsteps.对于固定的模拟长度达到一定的接受比。在大多数情况下，0.65的默认目标接受率是好的。

从估计的MAP点开始调优，以获得更有效的调优。

[smp, tuneinfo] = tuneSampler (smp,“开始”，mappars）;

绘制调优过程中步长的演变，以确保步长调优已经收敛。显示达到的合格率。

图;情节(tuneinfo.StepSizeTuningInfo.StepSizeProfile);包含(“迭代”）;ylabel (“步长”）;Accratio = tuneinfo.stepsizetuninginfo.acceptanceratio.

Accratio = 0.6700

画出样本

从后验密度中抽取样本，使用几个独立的链。为链选择不同的初始点，随机分布在估计的MAP点周围。指定从马尔可夫链开始丢弃的老化样本数量和老化后生成的样本数量。

设置“VerbosityLevel”在对第一链采样期间打印详细信息的值。

NumChains = 4;链=细胞(NumChains, 1);燃烧= 500;NumSamples = 1000;为c = 1: NumChains如果(c == 1) = 1;其他的水平= 0;结束链{c} = drawSamples (smp,“开始”MAPpars + randn(大小(MAPpars)),．..'伯恩'，伯恩，“NumSamples”，numsamples，．..“VerbosityLevel”，等级，“NumPrint”, 300);结束

| ================================================================================== ||磨练|log pdf |步长|num步骤|ACC比率|发散|| ================================================================================== || 300 | -1.490356e+02 | 2.617e-01 | 11 | 9.467e-01 | 0 | | 600 | -1.493478e+02 | 2.199e-02 | 3 | 9.367e-01 | 0 | | 900 | -1.509868e+02 | 2.244e-01 | 5 | 9.422e-01 | 0 | | 1200 | -1.493611e+02 | 1.166e-01 | 15 | 9.300e-01 | 0 | | 1500 | -1.488397e+02 | 2.617e-01 | 11 | 9.320e-01 | 0 |

检查收敛性诊断

使用诊断计算标准MCMC诊断的方法。对于每个采样参数，该方法使用所有链来计算这些统计信息：

后平均估计（意思）
估计蒙特卡罗标准错误（MCSE)，即后验均值估计值的标准差
后验标准差的估计(SD）
后缘分布的第5和95分位数的估计(Q5和Q95）
后验均值估计的有效样本量(ess.）
Gelman-Rubin收敛统计量(小红帽）.根据经验，值小红帽少于1.1被解释为链条融合到所需分布的标志。如果小红帽对于任何变量都大于1.1，然后尝试用drawSamples方法。

显示诊断表和示例开头定义的采样参数的真值。由于先验分布对该数据集没有提供信息，所以真实值在(或接近)第5和95分位数之间。

truePars = [trueIntercept;trueBeta;log(trueNoiseSigma^2)]

DIAGS = 4x8表名称意味着MCSE SD Q5 Q95 ESS RHAT_________________________________________________________________________________________________________________0.43983 0.0065001 0.34398 -1.0192 0.12916 2800.5 1.0001 {'x4'} -0.063367 0.11.08645 0.14159 -0.28853 0.1838 2443.3 1.0004 Truepars = 2 3 0 0

看到样品

调查融合和混合等问题以确定绘制的样本是否代表了目标分布的合理随机实现。要检查输出，请使用第一链绘制样品的迹线图。

的drawSamples方法从马尔可夫链的开始就丢弃老化样本，以减少抽样起点的影响。此外，轨迹图看起来像高频噪声，样本之间没有任何可见的长期相关性。这一行为表明该链混合良好。

图;次要情节(2 2 1);情节(链{1}(:1));标题(sprintf (“拦截,链1”））;为p = 2:1+NumPredictors subplot(2,2,p);情节(链{1}(:p));标题(sprintf (“β(% d),链1”p - 1));结束子图（2,2,4）;绘图（链条{1}（：，结束））;标题(sprintf (“LogNoiseVariance链1”））;

将这些链组合成一个矩阵，并创建散点图和直方图来可视化1-D和2-D的边缘后验分布。

concatenatedsamples = Vertcat（链条{：}）;图;plotmatrix（concatenatedsamples）;标题（'所有链子合并'）;

用于计算后部分布的功能

的logPosterior函数返回线性模型的正态似然和正态先验乘积的对数。输入参数参数有格式[拦截; beta; lognoisevariance]．X和Y包含预测器和响应的值。

的普通人函数返回具有均值的多元正态概率密度的对数亩和标准偏差σ，指定为标量或列向量相同的长度P．第二个输出参数是相应的渐变。

功能[logpdf, gradlogpdf] = logPosterior(参数，X,Y，．..InterceptPriorMean InterceptPriorSigma,．..Betapriormean，Betapriorigma，．..lognoisevariarcememean，lognoisevariaresigma）%解包参数向量拦截=参数（1）;beta =参数（2：结束-1）;lognoisevariance =参数（结束）;%计算对数似然及其梯度sigma = sqrt（exp（lognoisevariance））;mu = x * beta +拦截;z =（y  -  mu）/ sigma;loglik = sum（-log（sigma） -  .5 * log（2 * pi） -  .5 * z. ^ 2）;GradIntercept1 = SUM（Z / SIGMA）;gradbeta1 = x'* z / sigma;gradLognoisevariance1 = Sum（ - 。5 + .5 *（Z. ^ 2））;%计算日志优先级和梯度[lpittercept，gradIntercept2] =正常项目（拦截，拦截ProRiormean，拦截产品）;[Lpbeta，Gradbeta2] =正常项目（Beta，Betapriormean，Betaprigma）;[lplognoisevar，gradLognoisevariance2] =正常项目（Lognoisevariance，Lognoisevariancememem，LognoisevariAceSigma）;LogPrior = LpIntercept + LPBeta + Lplognoisevar;%返回后验对数及其梯度logpdf = loglik + logprior;GragIntercept = GradIntercept1 + GradIntercept2;Gradbeta = Gradbeta1 + Gradbeta2;GradLognoisevariance = GradLognoisevariance1 + GradLognoisevariance2;gradlogpdf = [gradtercept; gradbeta; gradLognoisevariance];结束功能[logpdf,gradlogpdf] = normalPrior(P,Mu,Sigma) Z = (P - Mu)./Sigma;logpdf =总和(日志(σ)- 5 *日志(2 *π)- 5 * (z ^ 2));gradlogpdf = z /σ;结束