基于MCMC算法的Logistic回归模型的参数估计——以智能睡眠手环为例

赵斯颖

（西安财经大学统计学院，西安 710000）

0 引言

近年来，随着马尔可夫链蒙特卡罗（MCMC）方法的普及，越来越多的学者采用这种算法来模拟复杂的多元分布。在实际工作、学习中经常会遇到数据缺失、信息不全或者非标准化的多维随机变量的抽样问题，一般的抽样方法很难对其进行处理，而MCMC方法就可以用来解决这类采样难的问题。MCMC作为一种随机抽样方法，被广泛应用于机器学习、深度学习和自然语言处理，其中包括Gibbs抽样法和M-H算法。使用MCMC方法估计模型中的参数，就是构造一个以参数的后验分布（目标分布）为平稳分布的马尔可夫链，然后利用从平稳分布中提取的样本点计算蒙特卡罗积分。

国内对模型的研究主要集中在其应用方面，对模型参数估计的研究相对较少，窦剑军等采用Gibbs抽样法对多元Logistic回归模型的参数进行估计，最终得到的参数后验分布为正态分布。蒋捷等利用MCMC算法中的Gibbs抽样来推断参数的后验均值，以此来估计模型的参数。王纯杰等采用M-H算法和切片Gibbs算法，计算Logistic回归模型的后验分布，结果表明这两种方法都具有可行性。刘贞以贝叶斯分析为背景，讨论了回归模型系数变点相关问题，并采用MCMC方法对得到的贝叶斯估计进行实证模拟。

在经典统计中，通常用最大似然估计和最小二乘法来估计逻辑回归模型的参数，在本文中，我们通过MCMC方法对Logistic回归模型做参数估计，在对参数的先验分布取正态分布时，推导出后验联合分布，并通过Python与MCMC相结合的方法对实际数据的拟合过程进行演示，对模型参数进行迭代，得到参数估计结果落在95%置信区间内，说明MCMC估计参数的结果是可靠的，从而证明模型的有效性。

1 模型描述

1.1 逻辑回归模型

逻辑回归（binomial logistic regression）是一个假设样本服从伯努利分布，利用梯度下降法和极大似然估计求解的二分类模型。其看似是一种回归算法，其实是分类算法，即已知许多特征，通过这些特征预测出类别的标签。要确定它属于哪个类别，就需要设置一个阈值，如果大于此阈值，则将其分配给一个类别；如果小于此阈值，则将其分配给另一个类别。二分类问题只有两个标签，记为条件概率（|），这里随机变量的值是实数，而随机变量的值通常是0或1，它的模型表示如下：

上式中∈R是输入，∈{}0,1是输出，∈R和∈R是参数，模型中的分类和值的正负是相关联的。

1.2 马尔可夫模型

马尔可夫模型是一种基于马尔可夫性质的统计模型，其假设在随机过程中给定某一时刻状态的概率分布，这个概率分布只与它前面的状态有关，所以马尔可夫链被广泛应用在时间序列模型中。

用数字描述就是假设某一时间序列是…X,X,X,X,X…，在+1时刻的条件概率只与时刻有关，这也是马尔可夫链模型一个非常有用的性质，即：

所以只要得到两个状态之间的转换概率，就能计算出对应的概率分布的样本。

马尔可夫模型还包含隐马尔可夫（HMM）、马尔可夫链蒙特卡罗（MCMC）和马尔可夫随机场（MRF），其中MCMC和MRF被广泛应用于模型参数的近似估计。

2 实证分析

2.1 睡眠数据分析

本文采集了智能手环中储存了两个多月的数据，将智能手环中的数据导出为csv文件格式，以便对数据进行建模。这里将智能手环的数据分为睡眠和清醒两种状态，手环主要根据摆动的频率、次数或者佩戴者心律的变化来判断睡眠和清醒状态。

通过绘制睡眠数据的散点图，我们可以更好地了解睡眠和清醒状态的时间分布，如图1所示，笔者通常在22点以后进入睡眠状态，而仅有少部分入睡时间在22点以前。根据实际情况，0点以后入睡的概率很高（接近1），而21点以前入睡的概率很小，所以散点图只描绘了从清醒到睡眠这段时间的状态，即从21点到0点的数据分布。

图1睡眠数据分布

图2 的清醒数据分布显示，笔者的清醒时间大都集中在早晨6点以后，在5点至6点之间，大部分是处于睡眠状态，比较符合实际情况。与上面睡眠数据分布同理，这里只展示从睡眠状态转换到清醒状态的时间段，即早上5点至8点的数据分布，而8点以后清醒的概率非常高（接近1）。

图2 清醒数据分布

通过分析可知，在本例中一共包含入睡、清醒两个状态和一个时间维度，对于这种问题，通常可以使用逻辑函数对其进行建模：

2.2 参数的先验分布

我们对逻辑函数的参数α和取不同的值来观察曲线的形状，结构如图3所示。

图3 不同α和β参数的逻辑曲线

由图3可知，当=0时，逻辑函数是一条型曲线，在曲线的中心位置增长速度较快；当＞0时，逻辑函数曲线向右偏移，由此只要找到合适的和参数就能很好地使用逻辑函数来拟合我们的数据。

因为参数的先验分布事先是未知的，所以需要使用贝叶斯方法来找到合适的和。通常假设这两个参数服从两个正态分布，正态分布有两个重要参数（均值）和（标准差），和的取值不同得到的正态分布曲线也不一样，的取值可以是正数也可以是负数；而的值通常取正数。这里用另一个参数（精度）来替代，的取值是的倒数，所以正态分布的标准差和精度成反比，标准差越大，精度越低；标准差越小，精度越高。正太分布概率密度函数的定义如下：

描绘出均值和精度完全不同的三组正态分布曲线，如图4所示。可以看出，均值决定了图形的中心对称位置；精度决定了图形的宽窄，取值越大，曲线越窄，分布精度越高；越小，曲线越宽，其分布精度也越低。

图4 正态分布

2.3 马尔可夫链蒙特卡罗

蒙特卡罗法又称随机抽样法，是在概率模型的基础上，通过计算机反复模拟得到近似解。而马尔可夫蒙特卡罗（简称MCMC）是以马尔可夫链为概率模型的蒙特卡罗法，它的核心是从一个已知概率密度函数的概率分布中进行采样，来构造最接近真实数据的概率分布。

由于参数和我们无法直接计算，所以需要通过生成和的数千个样本来构建最接近真实分布的近似值，逻辑函数的两个参数和服从正态先验概率，通过选择随机值，我们可以得到参数可能值的范围，如图5所示，越趋向于中心先验概率就越高。

图5 正态先验的参数空间

正态分布的期望值是：

正态分布的标准差是：

参数和服从正态分布，而正态分布的参数和具体数值未知，我们设定初始值=0、=0.05，利用MCMC对和进行抽样，最大化和的似然度。

有多种方法可用来诊断MCMC的收敛性，这里主要介绍两种检验方法。一种是Gelman检验，在给定置信水平为1的前提下，用每条链的轨迹区间1对应的长度除以整个轨迹区间长度作为诊断指数，当该指数接近1时就认为链收敛。另一种是轨迹图检验法，通过查看采样值路径来判断马尔可夫链蒙特卡罗链在迭代过程中是否达到稳定状态，收敛的马尔可夫链形似白噪声序列，在一个水平上没有趋势和周期的小幅波动。本文通过Pymc3提供的praceplot函数来绘制所有样本的轨迹，结果如图6所示。

图6轨迹图

图6 的轨迹图直观地展示了MCMC迭代过程中的采样值，将迭代次数设定为14000，可以看到马氏链迭代轨迹较平稳，并没有出现明显的波动周期性和异常值。

3 参数可视化

逻辑函数描述了从清醒到入睡的变化过程，但逻辑函数α和的真实值未知，通常假设α和来自于由和定义的正态分布，利用MCMC算法分别对参数α和采样和,使得到的和最能接近真实数据的分布，对清醒和睡眠建模为一个伯努利变量，清醒用1表示，入睡用0表示。特定时间下入睡的伯努利变量由逻辑函数来定义：

这里(t)为一个有独立时间变量的逻辑函数，所以上面的公式也可以写成：

MCMC的目标是使用现有的数据找到α和参数，并假定它们来自于正态先验。本文使用Python中功能强大的贝叶斯推理库Pymc3，该库包含了马尔可夫链蒙特卡罗和其他推理算法，通过代码创建模型并执行MCMC，为和抽取5000个样本，指定的抽样算法是Metropolic Hastings，我们将数据输入模型，并设定模型数据为伯努利变量，模型为数据找到最有可能的参数和，跟踪变量trace包含了所有后验分布的样本，据此可以对这些样本进行可视化，以探索它们在采样过程中的变化。随着样本的增加，MCMC算法收敛于最可能的值。在trace中返回的值都是参数的所有样本，我们可以利用直方图对这些值进行可视化，结果如图7所示。

图7 参数分布

从图7可以看到，α、的均值比较接近0，从前面描绘的逻辑函数曲线图来看，趋向于0时曲线趋于平坦，这说明清醒和睡眠的时间点存在重合的情况，这和实际情况相符，而α不等于0说明逻辑函数曲线发生了偏移，这也就意味着模型找到了一个从清醒到睡眠的转折点，从睡眠数据的分布图上看这个时间点应该是在晚上10点左右。

为了更加清楚地了解夜晚的睡眠状况，本文描绘出5000个样本数据的睡眠概率分布，如图8所示。

图8 睡眠概率分布

从图8可以看到，晚上9点半睡眠的后验概率均值大概在0.00至0.08之间，这说明笔者在晚上9点半就入睡的可能性比较小，这符合实际情况，大多数成年人一般不会在9点半就入睡，不过10点半的睡眠概率突然增加到了0.74，这也比较符合现实情况。从之前的分析可知，睡眠概率大于0.5的转折点是在晚上10点至10点半之间，也就是大部分成年人在这个时间点都会入睡，这应该比较符合常理。

4 结语

MCMC算法对整个统计学领域的研究有着深远的影响，通过MCMC随机模拟我们可以找到一个合适的模型来解决某些高维复杂的问题，特别是参数后验分布的计算。本文将MCMC方法运用于具体的实例中，以正态分布函数为建议分布函数，采用M-H抽样算法，运用Matlab对Logistic模型进行参数估计，通过构建睡眠模型让笔者了解了自己的睡眠方式。有了以上的经验，还可以根据其他数据信息，比如星期几、日常活动如何影响我们的睡眠，本文对此假设并未深入研究，但这是使用贝叶斯方法分析真实数据的良好开端。