Logistic回归模型的贝叶斯分析

打开脚本

这个例子展示了如何为逻辑回归模型做出贝叶斯推理slicesample．

统计推断通常基于最大似然估计(MLE)。MLE选择的参数使数据的可能性最大化，直观上很有吸引力。在极大误差估计中，假设参数是未知但固定的，并具有一定的置信度进行估计。在贝叶斯统计中，将未知参数的不确定性用概率量化，将未知参数视为随机变量。

贝叶斯推理

贝叶斯推理是将关于模型或模型参数的先验知识结合起来分析统计模型的过程。这种推论的根源是贝叶斯定理:

$$ P (\ mathrm{参数|数据})= \压裂# xA;{P (\ mathrm{数据|参数})\乘以P (\ mathrm{参数})}& # xA;{P (\ mathrm{数据})}& # xA;\propto \mathrm{likelihood} \times \mathrm{prior}$$

例如，假设我们有正常的观察结果

$X|\theta \sim N(\theta， \sigma^2)$

已知的先验分布是什么

$\theta \sim N(\mu， \tau^2)$

在这个公式中，有时称为超参数的mu和tau也是已知的。如果我们观察n的样本X，我们可以得到的后验分布

$$ \θ| X \ sim N \离开(\压裂{\τ^ 2}{\σ^ 2 / N + \ t ^ 2} \酒吧X + & # xA;\压裂{\σ^ 2 / n}{{\σ^ 2 / n + \ t ^ 2}} \μ,& # xA;\压裂{(\σ^ 2 / n) \ \τ^ 2}{\σ^ 2 / n + & # xA;\τ^ 2}\ $$

下图显示了的先验、似然和后验。

rng (0,“旋风”);n = 20;σ= 50;x = normrnd(10,σ,n, 1);μ= 30;τ= 20;Theta = linspace(- 40,100,500);日元= normpdf(意味着(x),θ,σ/√(n));y2 = normpdf(θ,μ,τ);postMean =τ^ 2 *意味着(x) /(τ)^ 2 +σ^ 2 / n) +σ^ 2 *亩/ n /(τ)^ 2 +σ^ 2 / n); postSD = sqrt(tau^2*sigma^2/n/(tau^2+sigma^2/n)); y3 = normpdf(theta, postMean,postSD); plot(theta,y1,“- - -”θ,y2,“——”,θ,y3,“-”。)传说(“可能性”，“之前”，“后”)包含(‘\θ）

汽车试验数据

在一些简单的问题中，如前面的正态均值推理例子，很容易求出封闭形式的后验分布。但在涉及非共轭先验的一般问题中，后验分布很难或不可能解析计算。我们将考虑逻辑回归作为一个例子。这个例子涉及到一个实验，帮助建模不同重量的汽车在里程测试中失败的比例。这些数据包括观察到的重量、测试的汽车数量和失败的汽车数量。我们将使用权重的转换版本来减少回归参数估计中的相关性。

一套汽车砝码Weight = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]';重量=(重量- 2800)/ 1000;%重定向和缩放%每种重量测试的汽车数量Total = [48 42 31 34 31 21 23 23 21 16 17 21]';%在每种重量下mpg性能较差的汽车数量Poor = [1 2 0 3 8 8 14 17 19 15 17 21]';

逻辑回归模型

逻辑回归是广义线性模型的一种特殊情况，适合于这些数据，因为响应变量是二项的。逻辑回归模型可以写成:

$$ P (\ mathrm{失败})= \压裂{e ^ {Xb}} {1 + e ^ {Xb}} $$

其中X为设计矩阵，b为包含模型参数的向量。在MATLAB®中，我们可以将这个方程写成:

logitp = @ (b, x) exp (b(1) +(2)。* x) / (1 + exp (b(1) +(2)。* x));

如果你有一些先验知识或一些非信息性的先验是可用的，你可以指定模型参数的先验概率分布。例如，在本例中，我们使用正常先验进行拦截b1和斜率b2,即

Prior1 = @(b1) normpdf(b1,0,20);拦截前百分比Prior2 = @(b2) normpdf(b2,0,20);斜率先验%

根据贝叶斯定理，模型参数的联合后验分布与似然与先验的乘积成正比。

Post = @(b) prod(binopdf(poor,total,logitp(b,weight))).．.%的可能性* prior1(b(1)) * prior2(b(2));%先验

注意，在这个模型中，后验的归一化常数在分析上是难以处理的。然而，即使不知道归一化常数，如果你知道模型参数的近似范围，你也可以可视化后验分布。

B1 = linspace(-2.5， -1, 50);B2 = linspace(3, 5.5, 50);simpost = 0 (50,50);为i = 1:长度(b1)为j = 1:长度(b2) simpost后(i, j) = ((b1(我),b2 (j)]);结束；结束；网格(b2, b1, simpost)包含(“坡”) ylabel (“拦截”) zlabel (“后验密度”)视图(-110,30)

这个后验在参数空间中沿着对角线被拉长，这表明，在我们看完数据后，我们相信参数是相关的。这很有趣，因为在我们收集任何数据之前，我们假设它们是独立的。这种相关性来自于我们的先验分布和似然函数的结合。

片抽样

在贝叶斯数据分析中，常用蒙特卡罗方法来总结后验分布。其思想是，即使你不能解析地计算后验分布，你可以从分布中生成一个随机样本，并使用这些随机值来估计后验分布或衍生的统计数据，如后验平均值、中位数、标准差等。切片抽样是一种设计用于从具有任意密度函数的分布中抽样的算法，该分布已知的密度函数不超过比例常数——这正是从一个归一化常数未知的复杂后验分布中抽样所需要的。该算法不生成独立样本，而是生成一个平稳分布为目标分布的马尔可夫序列。因此，切片采样器是一种马尔可夫链蒙特卡罗(MCMC)算法。然而，它不同于其他著名的MCMC算法，因为只需要指定缩放后验——不需要建议或边际分布。

本例展示了如何使用切片采样器作为里程测试逻辑回归模型贝叶斯分析的一部分，包括从模型参数的后验分布生成随机样本，分析采样器的输出，并对模型参数进行推断。第一步是生成一个随机样本。

Initial = [1 1];nsamples = 1000;跟踪= slicesample (nsamples初始,“pdf”的帖子,“宽度”20 [2]);

采样器输出分析

在从切片采样器中获得随机样本后，重要的是研究收敛和混合等问题，以确定样本是否可以合理地视为来自目标后验分布的一组随机实现。查看边缘轨迹图是检查输出的最简单方法。

次要情节(2,1,1)情节(跟踪(:1))ylabel (“拦截”);次要情节(2,1,2)情节(跟踪(:,2))ylabel (“坡”);包含(的样本数量);

从这些图中可以明显看出，参数起始值的影响需要一段时间才能消失(可能是50个左右的样本)，然后过程才开始看起来平稳。

在检查收敛性时，使用移动窗口来计算统计数据，如样本平均值、中位数或样本标准差，也是很有帮助的。这产生了一个比原始样本轨迹更平滑的图形，并且可以更容易地识别和理解任何非平稳性。

Movavg = filter((1/50)*ones(50,1)， 1, trace);次要情节(2,1,1)情节(movavg(: 1)包含(样品的数量) ylabel (“拦截手段”);次要情节(2,1,2)情节(movavg(:, 2))包含(样品的数量) ylabel (“坡度的含义”);

因为这些是50次迭代窗口的移动平均值，所以前50个值与图的其余部分是不可比较的。然而，每个图的其余部分似乎证实了参数后验均值在大约100次迭代后收敛到平稳。很明显，这两个参数是相互相关的，与前面的后验密度图一致。

由于磨合期代表的样本不能被合理地视为来自目标分布的随机实现，所以最好不要在片采样器输出的开始使用前50个左右的值。您可以只删除输出的那些行，但是，也可以指定一个“老化”周期。如果已经知道了合适的老化时间(可能是通过以前的测试)，这就很方便了。

跟踪= slicesample (nsamples初始,“pdf”的帖子,.．.“宽度”20 [2],“燃烧”, 50);次要情节(2,1,1)情节(跟踪(:1))ylabel (“拦截”);次要情节(2,1,2)情节(跟踪(:,2))ylabel (“坡”);

这些痕迹图似乎没有显示出任何非平稳性，表明老化期已经完成了它的工作。

然而，痕迹情节的第二个方面也应该被探索。虽然截距的轨迹看起来像高频噪声，但斜率的轨迹似乎有一个较低的频率成分，表明相邻迭代的值之间存在自相关性。我们仍然可以从这个自相关样本中计算均值，但通常通过删除样本中的冗余来减少存储需求是很方便的。如果这消除了自相关，它也允许我们将其视为独立值的样本。例如，您可以通过只保留每10个值来简化示例。

跟踪= slicesample (nsamples初始,“pdf”的帖子,“宽度”20 [2],.．.“燃烧”, 50岁,“薄”10);

为了检查这种细化的效果，从迹中估计样本自相关函数并使用它们来检查样本是否快速混合是有用的。

F = fft(去趋势跟踪,“不变”));F = F .* conj(F);ACF =传输线(F);: ACF = ACF (21);%保留延迟高达20。ACF = real([ACF(1:21,1) ./ ACF(1,1)].．.ACF(一21,2)。/ ACF(1、2)]);%正常化。Bounds =√(1/nsamples) * [2;2);% 95% CI为iid正常实验室= {“拦截ACF样本”，“斜坡ACF样本”}；为i = 1:2 subplot(2,1,i) lineHandles = stem(0:20, ACF(:，i)，“填充”，“r-o”);lineHandles。MarkerSize = 4;网格(“上”)包含(“滞后”{我})ylabel(实验室)在情节([0.5 - 0.5;20 20]， [bounds([1 1]) bounds([2 2])]，“- b”);图([0 20]，[0 0]，“- k”);持有从一个=轴;轴([1](1:3));结束

第一个滞后点的自相关值对于截距参数是显著的，对于斜率参数更是如此。为了进一步降低相关性，我们可以使用更大的细化参数重复采样。然而，出于本例的目的，我们将继续使用当前的示例。

模型参数的推断

正如预期的那样，样本的直方图模拟了后验密度的图。

次要情节(1,1,1)hist3(跟踪,(25、25));包含(“拦截”) ylabel (“坡”) zlabel (“后验密度”)视图(-110,30)

可以使用直方图或核平滑密度估计来总结后验样本的边缘分布特性。

次要情节(2,1,1)嘘(跟踪(:1))包含(“拦截”);次要情节(2,1,2)ksdensity(跟踪(:,2))包含(“坡”);

您还可以计算描述性统计数据，如后验均值或随机样本百分位数。为了确定样本大小是否足够大，以达到所需的精度，监视作为样本数量函数的迹的所需统计量是有帮助的。

csum = cumsum(跟踪);次要情节(2,1,1)情节(csum(: 1)“。/ (1:nsamples))包含(样品的数量) ylabel (“拦截手段”);次要情节(2,1,2)情节(csum(:, 2)’。/ (1:nsamples))包含(样品的数量) ylabel (“坡度的含义”);

在这种情况下，1000的样本量似乎足以为后验均值估计提供良好的精度。

bHat =意味着(跟踪)

bHat = -1.6931 4.2569

总结

统计和机器学习工具箱™提供了多种功能，允许您轻松指定可能性和先验。它们可以组合得到后验分布。的slicesample函数使您能够在MATLAB中使用马尔可夫链蒙特卡洛模拟进行贝叶斯分析。它甚至可以用于后验分布的问题，用标准随机数生成器难以抽样。