Logistic回归模型的贝叶斯分析

打开脚本

这个例子展示了如何使用逻辑回归模型进行贝叶斯推理slicesample．

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

贝叶斯推理

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

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

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

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

已知的先验分布是什么

$\theta \sim N(\mu， \tau²)$

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

$$ \θ| 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;= 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,“-”。)传说(“可能性”，“之前”，“后”)包含(‘\θ）

汽车试验数据

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

一套汽车重量重量= [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]';%每种重量下每加仑行驶里程表现较差的汽车数量Poor = [1 2 0 3 8 8 14 17 19 15 17 21]';

逻辑回归模型

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

$$ 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,即

if (ref (b1, 1) = 1 and ref (b1, 1) = 1 and ref (b1, 1) = 1 and ref (b1, 1) = 1 and% prior for interceptstickline (b2,0, 0)， colorred, linethick0;斜率前%

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

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。(1) / (1,1)．..ACF(一21,2)。/ ACF(1、2)]);%正常化。Bounds = sqrt(1/nsamples) * [2;2);% 95%置信区间为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 (“斜坡的方法”）;