趙斯穎

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

0 引言

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

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

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

1 模型描述

1.1 邏輯回歸模型

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

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

1.2 馬爾可夫模型

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

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

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

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

2 實(shí)證分析

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 正態(tài)分布

2.3 馬爾可夫鏈蒙特卡羅

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

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

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

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

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

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

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

圖6軌跡圖

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

3 參數(shù)可視化

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

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

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

圖7 參數(shù)分布

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

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

圖8 睡眠概率分布

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

4 結(jié)語

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