基于MCMC算法的Logistic回歸模型的參數(shù)估計
——以智能睡眠手環(huán)為例

趙斯穎

（西安財經(jīng)大學(xué)統(tǒng)計學(xué)院，西安 710000）

0 引言

近年來，隨著馬爾可夫鏈蒙特卡羅（MCMC）方法的普及，越來越多的學(xué)者采用這種算法來模擬復(fù)雜的多元分布。在實際工作、學(xué)習(xí)中經(jīng)常會遇到數(shù)據(jù)缺失、信息不全或者非標(biāo)準(zhǔn)化的多維隨機變量的抽樣問題，一般的抽樣方法很難對其進行處理，而MCMC方法就可以用來解決這類采樣難的問題。MCMC作為一種隨機抽樣方法，被廣泛應(yīng)用于機器學(xué)習(xí)、深度學(xué)習(xí)和自然語言處理，其中包括Gibbs抽樣法和M-H算法。使用MCMC方法估計模型中的參數(shù)，就是構(gòu)造一個以參數(shù)的后驗分布（目標(biāo)分布）為平穩(wěn)分布的馬爾可夫鏈，然后利用從平穩(wěn)分布中提取的樣本點計算蒙特卡羅積分。

國內(nèi)對模型的研究主要集中在其應(yīng)用方面，對模型參數(shù)估計的研究相對較少，竇劍軍等采用Gibbs抽樣法對多元Logistic回歸模型的參數(shù)進行估計，最終得到的參數(shù)后驗分布為正態(tài)分布。蔣捷等利用MCMC算法中的Gibbs抽樣來推斷參數(shù)的后驗均值，以此來估計模型的參數(shù)。王純杰等采用M-H算法和切片Gibbs算法，計算Logistic回歸模型的后驗分布，結(jié)果表明這兩種方法都具有可行性。劉貞以貝葉斯分析為背景，討論了回歸模型系數(shù)變點相關(guān)問題，并采用MCMC方法對得到的貝葉斯估計進行實證模擬。

在經(jīng)典統(tǒng)計中，通常用最大似然估計和最小二乘法來估計邏輯回歸模型的參數(shù)，在本文中，我們通過MCMC方法對Logistic回歸模型做參數(shù)估計，在對參數(shù)的先驗分布取正態(tài)分布時，推導(dǎo)出后驗聯(lián)合分布，并通過Python與MCMC相結(jié)合的方法對實際數(shù)據(jù)的擬合過程進行演示，對模型參數(shù)進行迭代，得到參數(shù)估計結(jié)果落在95%置信區(qū)間內(nèi)，說明MCMC估計參數(shù)的結(jié)果是可靠的，從而證明模型的有效性。

1 模型描述

1.1 邏輯回歸模型

邏輯回歸（binomial logistic regression）是一個假設(shè)樣本服從伯努利分布，利用梯度下降法和極大似然估計求解的二分類模型。其看似是一種回歸算法，其實是分類算法，即已知許多特征，通過這些特征預(yù)測出類別的標(biāo)簽。要確定它屬于哪個類別，就需要設(shè)置一個閾值，如果大于此閾值，則將其分配給一個類別；如果小于此閾值，則將其分配給另一個類別。二分類問題只有兩個標(biāo)簽，記為條件概率（|），這里隨機變量的值是實數(shù)，而隨機變量的值通常是0或1，它的模型表示如下：

上式中∈R是輸入，∈{}0,1是輸出，∈R和∈R是參數(shù)，模型中的分類和值的正負是相關(guān)聯(lián)的。

1.2 馬爾可夫模型

馬爾可夫模型是一種基于馬爾可夫性質(zhì)的統(tǒng)計模型，其假設(shè)在隨機過程中給定某一時刻狀態(tài)的概率分布，這個概率分布只與它前面的狀態(tài)有關(guān)，所以馬爾可夫鏈被廣泛應(yīng)用在時間序列模型中。

用數(shù)字描述就是假設(shè)某一時間序列是…X,X,X,X,X…，在+1時刻的條件概率只與時刻有關(guān)，這也是馬爾可夫鏈模型一個非常有用的性質(zhì)，即：

所以只要得到兩個狀態(tài)之間的轉(zhuǎn)換概率，就能計算出對應(yīng)的概率分布的樣本。

馬爾可夫模型還包含隱馬爾可夫（HMM）、馬爾可夫鏈蒙特卡羅（MCMC）和馬爾可夫隨機場（MRF），其中MCMC和MRF被廣泛應(yīng)用于模型參數(shù)的近似估計。

2 實證分析

2.1 睡眠數(shù)據(jù)分析

本文采集了智能手環(huán)中儲存了兩個多月的數(shù)據(jù)，將智能手環(huán)中的數(shù)據(jù)導(dǎo)出為csv文件格式，以便對數(shù)據(jù)進行建模。這里將智能手環(huán)的數(shù)據(jù)分為睡眠和清醒兩種狀態(tài)，手環(huán)主要根據(jù)擺動的頻率、次數(shù)或者佩戴者心律的變化來判斷睡眠和清醒狀態(tài)。

通過繪制睡眠數(shù)據(jù)的散點圖，我們可以更好地了解睡眠和清醒狀態(tài)的時間分布，如圖1所示，筆者通常在22點以后進入睡眠狀態(tài)，而僅有少部分入睡時間在22點以前。根據(jù)實際情況，0點以后入睡的概率很高（接近1），而21點以前入睡的概率很小，所以散點圖只描繪了從清醒到睡眠這段時間的狀態(tài)，即從21點到0點的數(shù)據(jù)分布。

圖1睡眠數(shù)據(jù)分布

圖2 的清醒數(shù)據(jù)分布顯示，筆者的清醒時間大都集中在早晨6點以后，在5點至6點之間，大部分是處于睡眠狀態(tài)，比較符合實際情況。與上面睡眠數(shù)據(jù)分布同理，這里只展示從睡眠狀態(tài)轉(zhuǎn)換到清醒狀態(tài)的時間段，即早上5點至8點的數(shù)據(jù)分布，而8點以后清醒的概率非常高（接近1）。

圖2 清醒數(shù)據(jù)分布

通過分析可知，在本例中一共包含入睡、清醒兩個狀態(tài)和一個時間維度，對于這種問題，通常可以使用邏輯函數(shù)對其進行建模：

2.2 參數(shù)的先驗分布

我們對邏輯函數(shù)的參數(shù)α和取不同的值來觀察曲線的形狀，結(jié)構(gòu)如圖3所示。

圖3 不同α和β參數(shù)的邏輯曲線

由圖3可知，當(dāng)=0時，邏輯函數(shù)是一條型曲線，在曲線的中心位置增長速度較快；當(dāng)＞0時，邏輯函數(shù)曲線向右偏移，由此只要找到合適的和參數(shù)就能很好地使用邏輯函數(shù)來擬合我們的數(shù)據(jù)。

因為參數(shù)的先驗分布事先是未知的，所以需要使用貝葉斯方法來找到合適的和。通常假設(shè)這兩個參數(shù)服從兩個正態(tài)分布，正態(tài)分布有兩個重要參數(shù)（均值）和（標(biāo)準(zhǔn)差），和的取值不同得到的正態(tài)分布曲線也不一樣，的取值可以是正數(shù)也可以是負數(shù)；而的值通常取正數(shù)。這里用另一個參數(shù)（精度）來替代，的取值是的倒數(shù)，所以正態(tài)分布的標(biāo)準(zhǔn)差和精度成反比，標(biāo)準(zhǔn)差越大，精度越低；標(biāo)準(zhǔn)差越小，精度越高。正太分布概率密度函數(shù)的定義如下：

描繪出均值和精度完全不同的三組正態(tài)分布曲線，如圖4所示。可以看出，均值決定了圖形的中心對稱位置；精度決定了圖形的寬窄，取值越大，曲線越窄，分布精度越高；越小，曲線越寬，其分布精度也越低。

圖4 正態(tài)分布

2.3 馬爾可夫鏈蒙特卡羅

蒙特卡羅法又稱隨機抽樣法，是在概率模型的基礎(chǔ)上，通過計算機反復(fù)模擬得到近似解。而馬爾可夫蒙特卡羅（簡稱MCMC）是以馬爾可夫鏈為概率模型的蒙特卡羅法，它的核心是從一個已知概率密度函數(shù)的概率分布中進行采樣，來構(gòu)造最接近真實數(shù)據(jù)的概率分布。

由于參數(shù)和我們無法直接計算，所以需要通過生成和的數(shù)千個樣本來構(gòu)建最接近真實分布的近似值，邏輯函數(shù)的兩個參數(shù)和服從正態(tài)先驗概率，通過選擇隨機值，我們可以得到參數(shù)可能值的范圍，如圖5所示，越趨向于中心先驗概率就越高。

圖5 正態(tài)先驗的參數(shù)空間

正態(tài)分布的期望值是：

正態(tài)分布的標(biāo)準(zhǔn)差是：

參數(shù)和服從正態(tài)分布，而正態(tài)分布的參數(shù)和具體數(shù)值未知，我們設(shè)定初始值=0、=0.05，利用MCMC對和進行抽樣，最大化和的似然度。

有多種方法可用來診斷MCMC的收斂性，這里主要介紹兩種檢驗方法。一種是Gelman檢驗，在給定置信水平為1的前提下，用每條鏈的軌跡區(qū)間1對應(yīng)的長度除以整個軌跡區(qū)間長度作為診斷指數(shù)，當(dāng)該指數(shù)接近1時就認為鏈?zhǔn)諗俊Ａ硪环N是軌跡圖檢驗法，通過查看采樣值路徑來判斷馬爾可夫鏈蒙特卡羅鏈在迭代過程中是否達到穩(wěn)定狀態(tài)，收斂的馬爾可夫鏈形似白噪聲序列，在一個水平上沒有趨勢和周期的小幅波動。本文通過Pymc3提供的praceplot函數(shù)來繪制所有樣本的軌跡，結(jié)果如圖6所示。

圖6軌跡圖

圖6 的軌跡圖直觀地展示了MCMC迭代過程中的采樣值，將迭代次數(shù)設(shè)定為14000，可以看到馬氏鏈迭代軌跡較平穩(wěn)，并沒有出現(xiàn)明顯的波動周期性和異常值。

3 參數(shù)可視化

邏輯函數(shù)描述了從清醒到入睡的變化過程，但邏輯函數(shù)α和的真實值未知，通常假設(shè)α和來自于由和定義的正態(tài)分布，利用MCMC算法分別對參數(shù)α和采樣和,使得到的和最能接近真實數(shù)據(jù)的分布，對清醒和睡眠建模為一個伯努利變量，清醒用1表示，入睡用0表示。特定時間下入睡的伯努利變量由邏輯函數(shù)來定義：

這里(t)為一個有獨立時間變量的邏輯函數(shù)，所以上面的公式也可以寫成：

MCMC的目標(biāo)是使用現(xiàn)有的數(shù)據(jù)找到α和參數(shù)，并假定它們來自于正態(tài)先驗。本文使用Python中功能強大的貝葉斯推理庫Pymc3，該庫包含了馬爾可夫鏈蒙特卡羅和其他推理算法，通過代碼創(chuàng)建模型并執(zhí)行MCMC，為和抽取5000個樣本，指定的抽樣算法是Metropolic Hastings，我們將數(shù)據(jù)輸入模型，并設(shè)定模型數(shù)據(jù)為伯努利變量，模型為數(shù)據(jù)找到最有可能的參數(shù)和，跟蹤變量trace包含了所有后驗分布的樣本，據(jù)此可以對這些樣本進行可視化，以探索它們在采樣過程中的變化。隨著樣本的增加，MCMC算法收斂于最可能的值。在trace中返回的值都是參數(shù)的所有樣本，我們可以利用直方圖對這些值進行可視化，結(jié)果如圖7所示。

圖7 參數(shù)分布

從圖7可以看到，α、的均值比較接近0，從前面描繪的邏輯函數(shù)曲線圖來看，趨向于0時曲線趨于平坦，這說明清醒和睡眠的時間點存在重合的情況，這和實際情況相符，而α不等于0說明邏輯函數(shù)曲線發(fā)生了偏移，這也就意味著模型找到了一個從清醒到睡眠的轉(zhuǎn)折點，從睡眠數(shù)據(jù)的分布圖上看這個時間點應(yīng)該是在晚上10點左右。

為了更加清楚地了解夜晚的睡眠狀況，本文描繪出5000個樣本數(shù)據(jù)的睡眠概率分布，如圖8所示。

圖8 睡眠概率分布

從圖8可以看到，晚上9點半睡眠的后驗概率均值大概在0.00至0.08之間，這說明筆者在晚上9點半就入睡的可能性比較小，這符合實際情況，大多數(shù)成年人一般不會在9點半就入睡，不過10點半的睡眠概率突然增加到了0.74，這也比較符合現(xiàn)實情況。從之前的分析可知，睡眠概率大于0.5的轉(zhuǎn)折點是在晚上10點至10點半之間，也就是大部分成年人在這個時間點都會入睡，這應(yīng)該比較符合常理。

4 結(jié)語

MCMC算法對整個統(tǒng)計學(xué)領(lǐng)域的研究有著深遠的影響，通過MCMC隨機模擬我們可以找到一個合適的模型來解決某些高維復(fù)雜的問題，特別是參數(shù)后驗分布的計算。本文將MCMC方法運用于具體的實例中，以正態(tài)分布函數(shù)為建議分布函數(shù)，采用M-H抽樣算法，運用Matlab對Logistic模型進行參數(shù)估計，通過構(gòu)建睡眠模型讓筆者了解了自己的睡眠方式。有了以上的經(jīng)驗，還可以根據(jù)其他數(shù)據(jù)信息，比如星期幾、日常活動如何影響我們的睡眠，本文對此假設(shè)并未深入研究，但這是使用貝葉斯方法分析真實數(shù)據(jù)的良好開端。