尚　華，馮　牧，張貝貝

(1. 首都經(jīng)濟(jì)貿(mào)易大學(xué) 統(tǒng)計(jì)學(xué)院 北京 100070; 2. 中國(guó)科學(xué)技術(shù)大學(xué) 管理學(xué)院,安徽 合肥 230000)

?

基于Bayesian方法的參數(shù)估計(jì)和異常值檢測(cè)

尚華1，馮牧2，張貝貝1

(1. 首都經(jīng)濟(jì)貿(mào)易大學(xué) 統(tǒng)計(jì)學(xué)院 北京 100070; 2. 中國(guó)科學(xué)技術(shù)大學(xué) 管理學(xué)院,安徽 合肥 230000)

摘要：異常值檢測(cè)是當(dāng)前數(shù)據(jù)分析研究中的一個(gè)重要研究領(lǐng)域。模型中的異常值會(huì)直接影響建模、參數(shù)的估計(jì)、預(yù)測(cè)等問(wèn)題?；谀Ｐ偷漠惓Ｖ禉z測(cè)，傳統(tǒng)的做法是先對(duì)模型參數(shù)進(jìn)行估計(jì)，再進(jìn)行異常值檢測(cè)。而異常值的存在會(huì)影響參數(shù)估計(jì)，從而導(dǎo)致下一步異常值檢測(cè)的不可靠；反之異常值檢測(cè)也會(huì)影響參數(shù)估計(jì)。針對(duì)這些不足之處，提出了基于Bayesian方法的參數(shù)估計(jì)和異常值檢測(cè)，此方法可以將參數(shù)估計(jì)和異常值檢測(cè)同時(shí)實(shí)現(xiàn)，具體做法是在線性回歸模型中引入識(shí)別變量，基于Gibbs抽樣算法，給出識(shí)別變量后驗(yàn)概率的計(jì)算方法，通過(guò)比較這些識(shí)別變量的后驗(yàn)概率進(jìn)行異常值定位，同時(shí)給出參數(shù)的估算方法。通過(guò)大量的模擬實(shí)驗(yàn)，結(jié)果表明，與傳統(tǒng)方法相比，提出的方法對(duì)異常值更靈敏。

關(guān)鍵詞：線性回歸；識(shí)別變量；參數(shù)估計(jì)；異常值；Bayesian方法；Gibbs抽樣

0引言

對(duì)現(xiàn)代統(tǒng)計(jì)數(shù)據(jù)分析來(lái)說(shuō)，異常值是一個(gè)普遍存在的問(wèn)題。一般來(lái)說(shuō)，異常值是指那些不同于數(shù)據(jù)中大部分?jǐn)?shù)據(jù)的一個(gè)或多個(gè)觀察值。在線性回歸中，把偏離線性模式的觀測(cè)值定義為異常值。Hampel等[1]認(rèn)為，數(shù)據(jù)中包含10% 的異常值是很正常的。Hubert[2]和McCann-Welsh[3]認(rèn)為數(shù)據(jù)中包含0.25n(n為數(shù)據(jù)量)個(gè)異常值是一個(gè)污染上界。

在應(yīng)用背景下，回歸分析是一個(gè)重要的統(tǒng)計(jì)工具。在回歸技術(shù)中，一般采用最小二乘估計(jì)(ordinary least squares,OLS)，這是因?yàn)樗膫鹘y(tǒng)性和易計(jì)算性。然而由于異常值的出現(xiàn)，會(huì)使得用最小二乘估計(jì)的參數(shù)出現(xiàn)很大的偏差甚至錯(cuò)誤，若再用于預(yù)測(cè)，可能會(huì)產(chǎn)生誤導(dǎo)性的結(jié)論。另外，在很多實(shí)際情況下，可能對(duì)異常值本身感興趣，例如銀行詐騙、腫瘤監(jiān)測(cè)以及報(bào)警系統(tǒng)等。不僅響應(yīng)變量可能是異常值，解釋變量也可能是異常值(杠桿點(diǎn))。這2類異常值都會(huì)使得傳統(tǒng)的最小二乘估計(jì)不可信。這里我們只討論響應(yīng)變量是異常值的情況。

基于模型的異常值檢測(cè)一般分為2步：第1步對(duì)模型中的參數(shù)進(jìn)行估計(jì)；第2步對(duì)殘差進(jìn)行評(píng)價(jià)，過(guò)大殘差對(duì)應(yīng)的觀測(cè)值為異常值。為減小異常值對(duì)參數(shù)估計(jì)的影響，有很多穩(wěn)健估計(jì)方法，例如最小截?cái)喽朔ǖ墓烙?jì)[4](least trimmed squares，LTS)、M-估計(jì)(minimum)[5]和S-估計(jì)(solution)[6]等。但是，參數(shù)估計(jì)與異常值檢測(cè)之間相互影響[7-8]。若第1步參數(shù)估計(jì)得不穩(wěn)健，會(huì)使得第2步異常值檢測(cè)的正確率降低；反之，異常值檢測(cè)的精度也能影響參數(shù)估計(jì)的穩(wěn)健性。

本文要探討的是用Bayesian方法來(lái)同時(shí)進(jìn)行參數(shù)估計(jì)和異常值檢測(cè)。在本文中，用Bayesian的思想和原理[9-11]，并引入Gibbs抽樣算法[12]，可以把參數(shù)估計(jì)和異常值檢測(cè)同時(shí)實(shí)現(xiàn)。而不是先估計(jì)參數(shù)，再檢測(cè)異常值，避免了參數(shù)估計(jì)和異常值檢測(cè)相互影響的問(wèn)題。為了驗(yàn)證該方法的正確性，本文進(jìn)行了大量的模擬實(shí)驗(yàn)，驗(yàn)證了該方法的可行性和有效性。

1基于Bayesian方法的參數(shù)估計(jì)和異常值檢測(cè)

多元線性回歸模型一般形式為

其中，p(給定)為解釋變量的數(shù)目；βj(j=0,1,…,p)為回歸系數(shù)；εi～N(0,σ2),i.i.d.i=1,2,…,n；β=(β0,β1,…,βp)T；σ2為未知參數(shù)。

對(duì)每個(gè)觀測(cè)值，引入異常值識(shí)別變量

記δi為第i個(gè)異常程度的大小。并且假設(shè)：每個(gè)觀測(cè)值受到異常擾動(dòng)的先驗(yàn)概率都為α，即p(δi=1)=α。

1)根據(jù)共軛先驗(yàn)分布的選取準(zhǔn)則和實(shí)際應(yīng)用需要，取參數(shù)的先驗(yàn)分布分別為

其中，μ，ξ，α，β*，V，ν和λ為超參數(shù)。

根據(jù)以上假設(shè)，加入了異常值的多元線性模型為均值轉(zhuǎn)移模型(mean-shift)。

為判定觀測(cè)值中是否含有異常值以及確定它們的閾值，構(gòu)造如下Bayesian假設(shè)檢驗(yàn)問(wèn)題。

根據(jù)Bayesian假設(shè)檢驗(yàn)的原理，當(dāng)H1對(duì)應(yīng)的后驗(yàn)概率p(γi=1|Y)，Y=(y1,y2,…,yn)T大于H0對(duì)應(yīng)的后驗(yàn)概率p(γi=0|Y)即p(γi=1|Y)>0.5時(shí)，認(rèn)為H1成立，從而認(rèn)為第i個(gè)觀測(cè)值為異常值；否則，認(rèn)為第i個(gè)觀測(cè)值為正常值。這樣，問(wèn)題就歸結(jié)為計(jì)算每個(gè)觀測(cè)值為異常值的后驗(yàn)概率p(γi=1|Y)。

2基于Gibbs抽樣的后驗(yàn)概率值的計(jì)算和參數(shù)的估計(jì)

2.1參數(shù)的全條件分布

由于后驗(yàn)概率p(γi=1|Y)涉及的分布比較復(fù)雜，下面引入Gibbs抽樣算法來(lái)解決這些后驗(yàn)概率值的計(jì)算問(wèn)題。為此，根據(jù)Bayesian定理可得下列全條件分布。

Y的聯(lián)合概率密度函數(shù)為

(1)

(1)式中：β=(β0,β1,…,βp)T；δ=(δ1,δ2,…,δn)T；γ=(γ1,γ2,…,γn)T。

1)由Bayesian定理，在Y,σ2,δ,γ給定時(shí)，β的全條件分布為

(2)

這里Xi=(1,x1i,…,xpi)T。

2)在Y,β,δ,γ給定時(shí)，σ2的全條件分布為

(3)

(3)式中：v1=n+v，

3)在Y,β,σ2,δ,γ(-j)給定時(shí)，γj的全條件分布為

(4)

(4)式中：

4)在Y,β,σ2,δ(-j),γ給定時(shí)，δj的全條件分布為

(5)

2.2識(shí)別變量后驗(yàn)概率值的計(jì)算

設(shè)β(r),(σ2)(r),δ(r),γ(r)，r=1,2,…,R為用Gibbs抽樣算法從上述全條件分布中抽取的樣本，則異常值的識(shí)別變量后驗(yàn)概率值的公式為

(6)

2.3參數(shù)β的估計(jì)

同理，參數(shù)β的估計(jì)值為

(7)

3線性模型異常值檢測(cè)的Bayesian方法的實(shí)施具體步驟

第1步確定先驗(yàn)分布中的超參數(shù)。例在本文中給出這些超參數(shù)的一組具體取值如下

第2步由Bayesian估計(jì)方法和超參數(shù)的取值，確定Gibbs抽樣的初值。

向量的上角標(biāo)(i,k)的含義為：該向量的第1個(gè)分量到第k-1個(gè)分量是第i+1次抽樣的樣本，第k個(gè)分量到最后一個(gè)分量為第i次抽樣抽取的樣本。例：(γ)(s-1,j)=((γ1)(s),…(γj-1)(s),(γj)(s-1),…,(γn)(s-1))T。

重復(fù)上述抽樣過(guò)程直到Markov鏈達(dá)到穩(wěn)定，取穩(wěn)定之后的R個(gè)Gibbs樣本

第4步按照(1)—(3)式計(jì)算識(shí)別變量的后驗(yàn)概率值，并按判定規(guī)則對(duì)異常值判定。

第5步按照(4)式估計(jì)參數(shù)β的大小。

4算例與分析

4.1算例

從三元正態(tài)分布β～N(β*,V-1)中隨機(jī)抽取一個(gè)向量為(-0.018,2.029,-1.028)T。

將上述向量只取整數(shù)為(0,2,-1)T。

故考慮模型

經(jīng)模擬產(chǎn)生100個(gè)數(shù)據(jù)。其中，(x1i,x2i)T，i=1,2,…,100來(lái)自于二元標(biāo)準(zhǔn)正態(tài)分布。

下面用3種方案進(jìn)行模擬和實(shí)驗(yàn)。

方案1在第19個(gè)觀測(cè)值上加上一個(gè)大小為-5的異常擾動(dòng)。

方案2在第20和79個(gè)觀測(cè)值上分別加上大小為的2.5，-4的異常擾動(dòng)。

方案3在第19至23個(gè)觀測(cè)值上分別加上大小為-6，4，9，-7，8的異常擾動(dòng)。

4.2異常值檢測(cè)

用Bayesian方法異常值識(shí)別變量的后驗(yàn)概率值如圖1所示。由圖1a可以看出，第19個(gè)觀測(cè)值為異常值的后驗(yàn)概率大于0.5，判定為異常值；由圖1b可以看出，第20個(gè)和第79個(gè)觀測(cè)值為異常值；由圖1c看出，從第19到23個(gè)觀測(cè)值均被檢測(cè)出為異常值。

圖1　基于方案1-3的識(shí)別變量后驗(yàn)概率值Fig.1　Posterior probabilities of these classification variables based on the schemes of 1-3

下面分別用3種穩(wěn)健的估計(jì)方法：LTS估計(jì)、M-估計(jì)和S-估計(jì)先估計(jì)參數(shù)，再進(jìn)行異常值判斷(用3sigma準(zhǔn)則)。對(duì)方案1，3種方法與Bayesian估計(jì)同樣都能檢測(cè)出第19個(gè)觀測(cè)值為異常值。對(duì)方案2，3種方法都能檢測(cè)出第79個(gè)觀測(cè)值為異常值，但是無(wú)法檢測(cè)出第20個(gè)觀測(cè)值為異常值。由此說(shuō)明，Bayesian方法對(duì)異常值更加敏感。對(duì)方案3，由于遮蔽現(xiàn)象，3種方法都能檢測(cè)出第19，21，22，23個(gè)觀測(cè)值為異常值，但都無(wú)法檢測(cè)出第20個(gè)觀測(cè)值為異常值。說(shuō)明Bayesian對(duì)異常值的遮蔽現(xiàn)象有一定作用，能有效地檢測(cè)出連續(xù)幾個(gè)在一起的異常值。

4.3參數(shù)估計(jì)

模型中參數(shù)β的真實(shí)值為(0,2,-1)T。3種方案參數(shù)估計(jì)的結(jié)果分別對(duì)應(yīng)于表1-3。從3個(gè)表中的數(shù)據(jù)可以看出，Bayesian方法估計(jì)出的參數(shù)值準(zhǔn)確度相對(duì)較高。

表1　基于方案1的4種方法的參數(shù)估計(jì)值

故用Bayesian方法來(lái)同時(shí)進(jìn)行參數(shù)估計(jì)和異常值檢測(cè)是行之有效的。

表2　基于方案2的4種方法的參數(shù)估計(jì)值

表3　基于方案3的4種方法的參數(shù)估計(jì)值

5結(jié)束語(yǔ)

異常值檢測(cè)是當(dāng)前數(shù)據(jù)分析研究中的一個(gè)熱點(diǎn)問(wèn)題。本文是在線性模型的基礎(chǔ)上，引入識(shí)別變量，基于Bayesian方法并結(jié)合Gibbs抽樣算法，給出識(shí)別變量的后驗(yàn)概率值的計(jì)算方法和參數(shù)的估算方法；同時(shí)估計(jì)模型參數(shù)和進(jìn)行異常值檢測(cè)。進(jìn)行了大量的模擬實(shí)驗(yàn)，結(jié)果表明，該方法對(duì)于解決線性模型數(shù)據(jù)中異常值檢測(cè)和參數(shù)估計(jì)是可行和有效的。

參考文獻(xiàn)：

[1]HAMPEL F R, RONCHETTI E M, ROUSSEEUW P J, et al. The Approach based on Influence Functions.[M]New York: John Wiley and Sons, 1986.

[2]HUBERT M, ROUSSEEUW P, VAN A S. High-breakdown robust multivariate methods[J]. Statistical science, 2008,23(1),92-119.

[3]MCCANN L, WELSCH R E. Robust Variable Selection Using Least Angle Regression and Elemental Set Sampling[J].Computational Statistics & Data Analysis, 2007, 52(1),249-257.

[4]ROUSSEEUW P J, VAN D K. Computing LTS regression for large data sets [J]. Data Mining and Knowledge Discovery, 2006,(12),29-45.

[5]HUBER P J.Robust Statistics[M].New York:Wiley,1981.

[6]DAVIES P L. Asymptotic behavior of S-estimates of multivariate location parameters and dispersion matrices[J]. Ann, Statist, 1987,(15),1269-1292.

[7]MARONNA R A, MARTIN D R, YOHAI V J. Robust Statistics:Theory and Methods,Chichester[M]: New York:Wiley, 2006.

[8]SHE Y, OWEN A B. Outlier detection using nonconver

penalized regression[J].Journal of the American Statistical Association, 2011,106(494),626-639.

[9]BERGER J O. Statistical decision theory and Bayesian analysis[M]. New York:Wiley,1985.

[10] 茆詩(shī)松．貝葉斯統(tǒng)計(jì)學(xué)[M].北京：中國(guó)統(tǒng)計(jì)出版社，1999.

MAO Shisong. Bayesian statistics[M]. Beijing: Statistics Press of China, 1999.

[11] 吳喜之.現(xiàn)代貝葉斯統(tǒng)計(jì)學(xué)[M].北京：中國(guó)統(tǒng)計(jì)出版社，2000.

WU Xizhi, Modern Bayesian statistics [M]. Beijing: Statistics Press of China, 2000.

[12] CHRISTIAN P R. Monte carlo statistical methods [M]. Berlin:Springer, 2004.

Parameter estimation and outliers detection based on Bayesian method

SHANG Hua1, FENG Mu2, ZHANG Beibei1

(1. College of Statistics, Capital University of Economics and Business, Beijing, 100070, P.R. China;2. College of Management, University of Science and Technology of China, Hefei, 230000, P.R. China)

Abstract:Outliers detection is an important research field in the current data analysis. Outliers in the data will affect the modeling, estimating parameters, forecasting and other issues directly. The conventional methods of outliers detection based on the model are to estimate the model parameters firstly, and then detect the abnormal value. The presence of outliers affects the parameter estimation, which results the in unreliability of outlier detection consequently; On the contrary, the presence of outliers will affect the parameter estimation. In this paper, we propose a new outliers detecting method based on Bayesian method, which can estimate parameters and detect outliers simultaneously. This method is introducing classification variables into linear regression model. Using Gibbs sampling a procedure for computing the posterior probabilities of classification variables and obtaining the estimation of parameters is designed. The outliers can be detected by comparing the posterior probabilities of these classification variables. A large number of simulation experiments illustrate that the proposed method is more sensitive to outliers compared with traditional methods.

Keywords：linear regression; classification variables；parameter estimation; outlier; Bayesian method; Gibbs sampling

DOI：10.3979/j.issn.1673-825X.2016.01.021

收稿日期：2015-03-12

修訂日期：2015-10-12通訊作者：尚華hnshanghua@tom.com

基金項(xiàng)目：國(guó)家自然科學(xué)基金(11426159)；首都經(jīng)濟(jì)貿(mào)易大學(xué)研究生科技創(chuàng)新項(xiàng)目(12013120061)

Foundation Items：The National Natural Science Foundation of China(11426159); The Postgraduate Technology Innovation Project of Capital University of Economics and Business (12013120061)

中圖分類號(hào)：TP391

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1673-825X(2016)01-0138-05

作者簡(jiǎn)介：

尚華(1981-)，女，河南新鄉(xiāng)人，講師，在讀博士，研究方向?yàn)閿?shù)理統(tǒng)計(jì)，數(shù)據(jù)挖掘。E-mail:hnshanghua@tom.com。

馮牧(1989-)，女，湖南岳陽(yáng)人，在讀博士，研究方向?yàn)闀r(shí)間序列，數(shù)據(jù)挖掘。

張貝貝(1983-)，女，河南濟(jì)源人，講師，博士，研究方向?yàn)闀r(shí)間序列，數(shù)據(jù)挖掘。

(編輯：張誠(chéng))