Logistic回归模型的贝叶斯分析

打开脚本

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

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

贝叶斯推理

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

$P（\mathrm{parameters}）=\frac ；{P（\mathrm{data}parameters}）乘以P（\mathrm{parameters}）} ；{P（\mathrm{data}）}$

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

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

其中σ已知，θ的先验分布为

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

在这个公式中，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=正常值（10，西格玛，n，1）；mu=30；tau=20；θ=linspace（-40100500）；y1=标准PDF（平均值（x），θ，西格玛/sqrt（n））；y2=normpdf（θ、μ、τ）；后均值=tau^2*平均值（x）/（tau^2+西格玛^2/n）+西格玛^2*mu/n/（tau^2+西格玛^2/n）；postSD=sqrt（τ^2*σ^2/n/（τ^2+σ^2/n））；y3=正常PDF（θ、邮差、邮差）；图（θ，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]';%每种重量下每加仑行驶里程表现较差的汽车数量差=[120388141719151721]；

逻辑回归模型

逻辑回归是广义线性模型的特例，适用于这些数据，因为响应变量是二项式的。逻辑回归模型可以写成：

$$ 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;跟踪=切片样本（初始样本、nsamples、，“pdf”的帖子,“宽度”20 [2]);

取样器输出分析

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

子地块（2,1,1）绘图（迹线（：，1））标签(“拦截”）;次要情节(2,1,2)情节(跟踪(:,2))ylabel (“坡度”); xlabel(“样本编号”）;

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

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

movavg=过滤器（（1/50）*一（50,1），1，迹线）；子批次（2,1,1）绘图（movavg（：，1））xlabel(“样本数量”) ylabel (“拦截手段”）;次要情节(2,1,2)情节(movavg(:, 2))包含(“样本数量”) ylabel (“斜坡的方法”）;

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

由于沉降期代表的样本不能被合理地视为来自目标分布的随机实现，所以可能不建议在切片采样器输出开始时使用前50个左右的值。您可以只删除输出中的这些行，但是，也可以指定一个“老化”周期。当合适的老化长度已经知道时(可能是从以前的运行中)，这是很方便的。

跟踪=切片样本（初始样本、nsamples、，“pdf”的帖子,...“宽度”,[20 2],“燃烧”,50); 子地块（2,1,1）绘图（迹线（：，1））标签(“拦截”）;次要情节(2,1,2)情节(跟踪(:,2))ylabel (“坡度”）;

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

然而，还应该探讨跟踪图的第二个方面。虽然截距的轨迹看起来像高频噪声，但斜率的轨迹似乎有一个较低的频率分量，这表明在相邻迭代的值之间存在自相关。我们仍然可以从这个自相关的样本中计算平均值，但通常通过去除样本中的冗余来降低存储需求是很方便的。如果这消除了自相关性，我们也可以将其视为独立值的样本。例如，您可以通过只保留每10个值来稀释样本。

跟踪=切片样本（初始样本、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 = sqrt(1/nsamples) * [2;2);% 95%置信区间为iid正常实验室= {“用于截取的样本ACF”,“斜坡的ACF样本”}；对于i=1:2子批次（2,1，i）线手柄=阀杆（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]）；xlabel(“拦截”) 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））xlabel(“样本数量”) ylabel (“斜坡的方法”）;