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

打开脚本

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

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

$$ P(\θ| y) = \压裂{P (y | \θ)(\θ)}{P (y)} $$

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

线性回归模型

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

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

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

$$ mu_i = + x_i^T， $$

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

$美元\α\ sim \ mathcal {N} (\ alpha_0 \ sigma_ \α^ 2)美元$ ，

$美元\ beta_i \ sim \ mathcal {N} (\ beta_ {0}, \ sigma_{\β}^ 2)美元$ ．

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

噪声方差是一个(平方)尺度参数，只能是正的。这样就可以更容易、更自然地将它的对数视为自由参数，它是无界的。在噪声方差的对数之前指定一个正态分布:

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

写出自由参数的后验密度的对数 $$ \θ=(\α;\β;日志\ \σ^ 2)美元$ 作为

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

忽略常数项，取最后两项的和 $g(\θ)美元$ ．要使用HMC，请创建一个计算函数句柄 $g(\θ)美元$ 及其梯度 ${\部分g} /{\部分\θ}$ 对于任意值 $\θ美元$ ．用于计算的函数 $g(\θ)美元$ 位于脚本的末尾。

创建数据集

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

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

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

NumData = 100;rng (“默认”）%的再现性X =兰特(NumData NumPredictors);= X*trueBeta + trueIntercept;y = normrnd(μ,trueNoiseSigma);

定义后验概率密度

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

InterceptPriorMean = 0;InterceptPriorSigma = 10;BetaPriorMean = 0;BetaPriorSigma = 10;LogNoiseVarianceMean = 0;LogNoiseVarianceSigma = 2;

保存一个函数logPosterior在MATLAB®路径上返回先验和似然乘积的对数，以及这个对数的梯度。的logPosterior函数在本例的最后定义。然后，调用带参数的函数来定义logpdf的输入参数hmcSampler函数。

logpdf = @(参数)logPosterior(参数,X, y,．..InterceptPriorMean InterceptPriorSigma,．..BetaPriorMean BetaPriorSigma,．..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] = estimateMAP (smp,“VerbosityLevel”, 0);MAPIntercept = mapservers (1) MAPBeta = mapservers (2:end-1) mapservers (2:end-1)

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

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

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

调整取样器

为获得有效的采样，选择合适的采样参数值是非常重要的。找到好的值的最好方法是自动调优MassVector，StepSize,NumSteps参数使用tuneSampler方法。使用该方法:

1.调优MassVector取样器。

2.调优StepSize和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);结束

|==================================================================================| | ITER日志PDF | |步长| NUM | | ACC比例不同的步骤  | |==================================================================================| | 300 e + 02 | -1.490356 | 2.617 e-01 | 11 | 9.467 e-01 | 0 | | 600 | -1.493478 e + 02 | 2.199 e-02 | 3 | 9.367 e-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)]

诊断接头= 4×8表名的意思是MCSE SD Q5 Q95 ESS小红帽  ______ _________ _________ _______ ________ _______ ______ ______ {' x1} 2.3792 0.0055101 0.27821 1.9116 2.8372 2549.3 1.0009{“x2”}2.5533 0.0062076 0.32836 2.0223 3.0959 2798.1 1.0004 {x3的}-0.43983 0.0065001 0.34398 -1.0192 0.12916 2800.5 1.0001{“x4”}-0.063367 0.0028645 0.14159 -0.288530.1838 2443.3 1.0004 truePars = 2 3 00

看到样品

研究收敛和混合等问题，以确定抽取的样本是否代表了目标分布的一组合理的随机实现。要检查输出，请使用第一个链绘制样本的跟踪图。

的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函数返回线性模型的正态似然和正态先验乘积的对数。输入参数参数有格式[拦截;β;LogNoiseVariance]．X和Y分别包含预测器和响应的值。

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

函数[logpdf, gradlogpdf] = logPosterior(参数，X,Y，．..InterceptPriorMean InterceptPriorSigma,．..BetaPriorMean BetaPriorSigma,．..LogNoiseVarianceMean LogNoiseVarianceSigma)%解包参数向量拦截=参数(1);β=参数(2:end-1);LogNoiseVariance =参数(结束);%计算对数似然及其梯度σ=√exp (LogNoiseVariance));= X* +截距;Z = (Y - Mu)/;loglik =总和(日志(σ)- 5 *日志(2 *π)- 5 * z ^ 2);gradIntercept1 =总和(Z /σ);gradBeta1 = X ' * Z /σ;gradLogNoiseVariance1 =(-求和。5 + 5 * (z ^ 2));%计算日志优先级和梯度[LPIntercept, gradIntercept2] = normalPrior(Intercept,InterceptPriorMean,InterceptPriorSigma);[LPBeta, gradBeta2] = normalPrior(Beta,BetaPriorMean,BetaPriorSigma);[LPLogNoiseVar, gradLogNoiseVariance2] = normalPrior(LogNoiseVariance,LogNoiseVarianceMean,LogNoiseVarianceSigma);logprior = LPIntercept + LPBeta + LPLogNoiseVar;%返回后验对数及其梯度Logpdf = loglik + logprior;gradIntercept = gradIntercept1 + gradIntercept2;gradBeta = gradBeta1 + gradBeta2;gradLogNoiseVariance = gradLogNoiseVariance1 + gradlognoisvariance2;gradlogpdf = [gradIntercept; gradBeta gradLogNoiseVariance];结束函数[logpdf,gradlogpdf] = normalPrior(P,Mu,Sigma) Z = (P - Mu)./Sigma;logpdf =总和(日志(σ)- 5 *日志(2 *π)- 5 * (z ^ 2));gradlogpdf = z /σ;结束