一種加入類間因素的曲線聚類算法

許騰騰，王瑞，黃恒君

（蘭州財經大學 統計學院，甘肅 蘭州 730020）

隨著信息技術的不斷發展，數據獲取越來越便捷，數據采集的密集化程度也越來越高。隨之出現了一種具有函數特征的數據類型。如心理學研究中的腦電信號數據、生物技術中的基因微序列數據、化學計量中的光譜分析數據、經濟研究中的股票分時成交價數據、環境監測中的污染物濃度數據等，均隨著時間變化而表現出明顯的曲線特征。當前文獻中將這種數據類型稱為函數型數據(functional data)[1]。

一般而言，函數型數據的曲線形式無法直接獲取，通常僅能夠觀測到其離散樣本點，并針對離散數據進行傳統多元統計分析。當然，這種做法由于未能考慮到數據的函數特性(如連續、可導等)，同時需要處理高維問題，往往不能取得很好的分析效果[2]。為此，針對數據的曲線特征，人們提出了各種分析方法[3-4]，包括函數型主成分分析、函數型線性模型、函數型聚類分析等，將有限維的多元分析推廣到無限維的函數型數據上來。

聚類分析是數據探索、數據壓縮和展現的重要工具，本文討論函數型數據的聚類算法。目前，函數型數據聚類分析方法大致分為兩類[3]：一是原始數據法，該類方法直接針對離散樣本點進行聚類，屬于高維數據分析方法，文獻[5]對這種做法進行了綜述。二是投影方法，即以有限維的基底函數逼近曲線，將無限維的問題轉化為有限維問題展開分析。依據對基底函數所對應的系數處理方式不同，又可分為濾波法和自適應法。濾波法將基底函數所對應的系數設定為固定參數，分曲線擬合和聚類分析兩步展開：首先以有限維基底擬合曲線，對估計的參數執行傳統聚類算法。包括利用自組織映射(SOM)算法進行聚類和擬合函數型數據[6]；利用兩階段隨機過程分別完成數據降維和聚類[7]等。根據基底函數選擇利用B-樣條基底函數擬合數據并根據傳統聚類方法分析[8-9]，利用正交基函數進行聚類分析[10]等。自適應法是將基底函數所對應的系數作為隨機變量處理，將曲線擬合和聚類分析納入一個目標函數，采用類似EM的算法，同時進行優化。如利用機器學習和神經網絡模型SVR分析時空數據[11]、利用STM算法對時空數據進行聚類[12]以及時間序列數據[13]、經維度數據[14]等的聚類方法；基于K-medoids項目聚類的協同過濾推薦算法[15]；基于多元函數型主成分分析(FPCA)方法進行的改進混合模型同時進行曲線擬合與聚類分析[16]。在隨機過程中利用K-L散度度量，采用類似于EM算法進行聚類的算法[17]等。

盡管有眾多其他的算法[6,18]，目前的函數型聚類分析僅考慮了類內部的差異，而忽視了類間的差異性。對傳統離散數據的聚類研究表明[19]，同時考慮類內與類間差異有助于提升聚類效果。

受此啟發，本文提出一種加入類間因素的曲線聚類算法。本文的做法屬于濾波法，包括B-樣條基底擬合曲線、曲線距離確定、曲線聚類目標函數設定，以及加入類間因素的曲線聚類算法等。

1 加入類間因素的曲線聚類

根據前面的描述，本文討論的曲線聚類分析模型構建主要包含3個方面：1)由觀測離散數據生成函數型數據，這里采用B-樣條基底表述的方法；2)構造曲線之間的“距離”的表述，并通過用B-樣條基底系數，將曲線距離轉化為傳統歐氏距離；3)以所構造的距離作為親疏程度度量，并構建同時考慮類內差異和類間差異的目標函數，完成曲線聚類任務。

1.1 曲線生成

假 定 數 據 Yi=[yi1yi2···yim]T(i=1,2,···,n)由 如下模型生成：

式中： Xi(t)表 示待估計曲線，ε表示服從零均值同方差的獨立同分布隨機變量。進一步假定 Xi(t)可由一組基底 {φi1(t),φi2(t),···}表示，則有

稱這種做法為基底函數法，它是一種將離散數據轉化為曲線的常用平滑技術[3]。對待估計曲線 Xi(t)采取截斷處理，得到如式(3)的形式：

從而將無限維問題轉化為有限維估計方式。進一步假定[9]：

1) 對不同曲線 Xi(t)(i=1,2,···,n)采用一組相同的基底表述；

2) 基底函數設定為等距節點B-樣條基底。有

式中： L =K+M ， Bk,M(t)表 示第 k個內部節點數量為K 的 M 階 B-樣條基底函數。 BM(t)表示M階B-樣條基底函數。對于參數 αi，我們利用最小二乘法進行估計。

1.2 曲線距離

假定曲線 Xi(t)為 L2空 間的元素。則根據 L2范數定義，有曲線 Xi(t)和 Xj(t)的距離為

其中 ‖· ‖表示L2范數，由假定1)及式(4)知

結合式(6)，式(5)可轉化為

式中

其中，K 為 L ×L 實對稱矩陣， K中的每一個元素?Bi,Bj表示 L2空間的內積。但是類似于d2(i,j)=[αi-αj]T[αi-αj]這種形式的距離公式并不適用于非正交基底函數[9]，為將曲線距離用傳統距離公式表示，對 K 作楚列斯基(Cholesky)分解得 K =LLT，其中 L 為上三角矩陣，并令 bi=LTαi，式(7)可表示為

需要說明的是，式(8)完成了從曲線距離到一般距離的轉變，構成了將曲線聚類轉化為傳統多元聚類問題的基礎。利用式(8)，運用傳統聚類算法對 bi進 行聚類，得到 P 類 ，記為 i ∈ Gp(p ∈ 1,2,···,P)。

由 bi=LTαi得 到 B =AL ， 其中 A =[α1α2···αnp]T，B=[b1b2···bn]T。np表示第 Gp類中的曲線數量，p令 Xˉ (t0)表示隨機選取的一條曲線作為初始類中心，ˉ(Gp)(t)表示第 Gp類中的類中心。則有Xˉ(Gp)(t)=np-11TBL-1BM(t)。

1.3 改進的曲線聚類算法

聚類分析的目的是將同類型數據進行歸類，同時對不同類型的數據進行區分。文獻[19]針對傳統離散數據提出的K-means聚類擴展方法兼顧了類內、類間差異。具體來講，通過對數據集引入全局中心點實現類內差異最小化的同時類中心與全局中心點距離最大化。相比于K-means算法，這種做法提高了聚類效果[19]。

受此啟發，本文將K-means聚類分析擴展到函數型聚類分析上。本文的曲線聚類目標函數為

式中： Φ 表示待估參數矩陣(A或B)，U表示由uik構成的矩陣，其中 uik∈{0,1}，uik=1， Xi(t)表示曲線，(t0)表示隨機選取的一條曲線作為初始類中心，結合式(4)的曲線基底表述，得到目標函數：

根據前面關于曲線距離的描述將式(7)～(8)代入式(10)得到

式中： b?=LTα?， b0=LTα0， α?表示第k類類中心對應的參數， α0表示初始類中心曲線的參數。

目標函數確定后，式(11)中含有兩個未知參數 α 及 U 。通過固定一項求解另一項的步驟來求解式 (11)，即()

2) 固定U=U，求解函數FΦ,U。

針對1)，為使目標函數式(11)達到最小，當目標函數分子中曲線與對應類中心曲線距離小時uik=1，否則為 0，即

針對2)，假設 b?已知，對目標函數式(11)關于b?求偏導數：

得出

進一步化簡得到b?

即

在進行計算機編程時可以不斷對步驟1)、2)進行迭代，直至找出最優 U 和 Φ 。算法流程如下：

Input：X={X1,X2,···Xn},k

Initialize： Randomly choose an initialb0=b1,b2,···,bk

Repeat

Fixed Φ, use eq. (12) to solveU

Fixed U, use eq. (13) to solveΦ

Until convergence.

進一步，由 b?=LTα?， 求解出 b?可 得到參數 α?，并根據式(4)還原出類中心曲線。

2 算法效果模擬驗證與分析

為驗證本文曲線聚類算法的效果，利用模擬數據與文獻[9]中曲線聚類方法進行比較。模擬數據由兩組高斯分布生成兩類曲線構成。模擬過程中兩類高斯分布均值取0.5和1，方差取0.7和1。在確定類別的前提下比較本文算法與文獻[9]曲線聚類算法的聚類效果。聚類效果評價指標采用蘭德指數(Rand index)評價算法的性能[20]。同時分析兩組高斯分布的參數(均值和方差)對聚類的影響。分析結果顯示：同均值異方差情況下兩種曲線聚類方法聚類結果均存在一定的誤判，異均值異方差情況下二者聚類也存在誤判，異均值同方差情況下二者聚類未出現誤判。以下針對這一現象做出分析。

該部分采用R軟件進行數據模擬分析，每組包含 n條 數據，每條數據含有 m個數據點，則模擬數據中每組高斯分布要生成 m ×n個隨機數。為保證擬合結果的光滑，內部節點采用等距節點設置方式。針對高斯分布中的均值和方差分別在同均值異方差、同方差異均值、異均值異方差情況下分析本文的曲線聚類方法與已有曲線聚類方法的效果，并對相應結果進行分析。為便于表述，兩類模擬數據分別記為1類和2類，生成的區間長度設置為12。為便于展示，本文以圖1異均值異方差條件下兩種聚類方法比較為例。

圖1 模擬數據曲線聚類對比Fig.1 Comparison with simulated data of curve's clustering

圖1 表明：兩組高斯分布參數不同條件下，本文方法與文獻[9]相比，圖1(b)中1類曲線分布密集程度大于圖1(a)中1類曲線。為避免模擬次數少或其他原因對聚類效果的影響，對3種類型的數據分別模擬一萬次，比較兩種方法的平均錯判率，定義錯判率=abs(1類個數-n)/n，模擬驗證中m=12，n=50，錯判率下降比例=文獻[9]方法錯判率-本文方法錯判率。結果見表1。

表1、2表明：無論本文的曲線聚類還是文獻[9]中的曲線聚類方法，類中心的變化與高斯分布中均值有關，而聚類效果好壞與高斯分布的方差有關。對比表1、2中的同均值異方差和異均值異方差的錯判率及蘭德指數可以得出：當兩類高斯分布均值相同，方差不同時，兩種方法對應的蘭德指數相比于其他類型數據偏低。同時方差因素對聚類效果也會產生影響。綜合比較表1、2中的3類數據錯判率及蘭德指數，可以得到：對于曲線聚類分析，聚類效果會同時受數據總體均值和方差的影響，對比分析表1、2均值相同方差不同的情形，可以得到：均值對聚類的影響程度要大于方差，同時表1、2對兩種方法錯判率對比結果顯本文的方法能夠降低聚類錯判率從而提高聚類效果。

表1 3種類型模擬數據平均錯判率Table1 Average error rate of three types' simulated data

表2 3種類型模擬數據蘭德指數Table2 Rand index of three types' simulated data

3 NO2小時濃度曲線聚類效果分析

空氣質量，不僅關乎人類生存質量，同時也是衡量可持續發展能力和宜居程度的重要指標。NO2是一種重要的機動車尾氣污染物，其污染程度涉及人們生活出行的健康。近年來，空氣質量問題引起人們廣泛的關注，大氣污染監測數據成為人們了解空氣質量的客觀途徑，也構成空氣質量統計分析的數據基礎。

作為示例，通過實時網絡爬蟲手段[21]，采集蘭州市鐵路設計院空氣質量監測站(交通污染控制點)的NO2小時濃度數據，采用本文的曲線聚類算法展開大氣污染等級聚類分析，并與傳統曲線聚類結果進行比較。我們分析的樣本期為2013年6月1日—10月14日。

根據前面的方法，采用B-樣條基底函數進行曲線聚類分析。為保證擬合結果光滑，兩種聚類方法樣條基底階數M均設置為5，節點采用等距節點設置為11(文中采用廣義交叉驗證準則進行節點數量選擇)。考慮相同類中心下，與文獻[9]曲線聚類進行聚類效果對比，如圖2所示。

圖2表明，K=5時類中心聚類效果優于K=4，即隨著類中心個數的增加，兩種方法的聚類效果均有所提升，說明類中心個數的確定在曲線聚類中起到關鍵作用。但需要指出的是，本文方法的類中心分布曲線更為平滑，類間的類中心曲線分布更為分散，進一步說明本文提出的方法聚類效果優于已有聚類方法。此外，考慮到實際應用，可將圖2中的不同類別曲線看作空氣質量污染物等級劃分[20]。對比圖 2(a)、(c)與圖 2(b)、(d)可以發現，在空氣質量實時監測過程中，圖2(a)、(c)出現不同等級交叉情況，這對空氣質量等級劃分及應對會造成影響[22]。圖2(b)、(d)在進行空氣質量分析過程中能夠較好的對空氣質量進行聚類。另外，相比于針對離散數據的傳統K-means聚類分析[23]，本文方法能夠實時檢測NO2小時濃度變化趨勢，并依據該變化趨勢對污染物進行等級劃分。

為便于展示，本文以K=5的曲線聚類結果為例，結果見圖3。圖3表明，相比于已有曲線聚類算法，利用本文曲線聚類算法類內曲線分布集中，類間差異化明顯。這與圖2中兩種曲線聚類算法類中心比較結果相一致。說明本文方法具有較好的類間區分度。

為進一步驗證本文曲線聚類的聚類效果，對兩種方法的分類精確度采用公式：類間差異/(類內差異+類間差異)進行對比，見圖4。圖4表明，隨著類中心個數的增加，兩種曲線聚類算法聚類效果均有所提高。本文曲線聚類的聚類效果要好于文獻[9]的方法。通過與文獻[9]方法進行比較，本文方法在4類的聚類效果低于3類聚類效果，隨著類中心個數大于4類，聚類效果才逐步隨著類中心個數增加聚類效果不斷提升。說明本文方法存在一定的不穩定性。

圖2 曲線聚類類中心對比Fig.2 Comparison with curve cluster's center generated by different algorithms

圖3 NO2小時濃度數據曲線聚類對比Fig.3 Comparison with curve clustering of NO2 concentration

圖4 聚類效果對比結果Fig.4 Comparison with clustering effects

4 結束語

本文基于已有曲線聚類方法，針對聚類效果不明顯的問題，提出加入類間因素的擴展曲線聚類算法。加入類間因素能夠同時保證兩類數據類內差異較小和類間差異較大。模擬數據及實例應用表明，本文的曲線聚類算法有助于提高聚類效果。

需要說明的是，本文的目的是將同時考慮類內和類間差異的做法引入曲線聚類算法。但我們的做法屬于兩步法，即首先擬合曲線，然后進行聚類。這種做法很難達到兩部分的統一優化[24]。為此，后續的工作是，在同時考慮類內和類間差異的情況下，進行自適應算法研究，即將曲線擬合和聚類分析納入一個目標函數，同時進行優化。