范朝霞，趙 明，楊劍鋒

(1.貴州大學(xué) 理學(xué)院，貴州 貴陽 550025;2.耶夫勒大學(xué) 技術(shù)與可持續(xù)發(fā)展學(xué)院，瑞典 80176;3.貴州理工學(xué)院 信息工程學(xué)院，貴州 貴陽550003)

基于CU變換的非齊次泊松過程的統(tǒng)計(jì)驗(yàn)證模型

范朝霞1，趙 明2，楊劍鋒3*

(1.貴州大學(xué) 理學(xué)院，貴州 貴陽 550025;2.耶夫勒大學(xué) 技術(shù)與可持續(xù)發(fā)展學(xué)院，瑞典 80176;3.貴州理工學(xué)院 信息工程學(xué)院，貴州 貴陽550003)

開源軟件中故障輸入數(shù)據(jù)的非齊次泊松過程(NHPP)的檢驗(yàn)和均值函數(shù)的估計(jì)是排隊(duì)論分析中重要的子問題。本文依據(jù)均值函數(shù)的連續(xù)性和緩慢變化性，把NHPP轉(zhuǎn)化為小區(qū)間泊松性的分段常數(shù)非齊次泊松過程(PC NHPP)進(jìn)行研究，主要運(yùn)用條件均勻(CU)變換法對(duì)數(shù)據(jù)調(diào)整變換后，驗(yàn)證泊松性。本文還運(yùn)用Bugzilla的錯(cuò)誤報(bào)告數(shù)據(jù)對(duì)模型實(shí)踐分析，結(jié)果顯示CU變換后的驗(yàn)證效果可以接受。最后運(yùn)用R語言擬合參數(shù)，計(jì)算出故障輸入過程參數(shù)，為排隊(duì)論模型構(gòu)造和求解奠定基礎(chǔ)。

NHPP檢驗(yàn)；CU變換；K-S檢驗(yàn)； R語言

針對(duì)軟件排錯(cuò)過程中故障輸入過程的統(tǒng)計(jì)分析方法[1]，國內(nèi)目前應(yīng)用廣泛的是假設(shè)錯(cuò)誤輸入過程是NHPP，然后根據(jù)三類常見累積強(qiáng)度函數(shù)模型：線性模型、冪率模型、對(duì)數(shù)線性模型等，擬合比較最優(yōu)的累積強(qiáng)度函數(shù)模型，求解參數(shù)，運(yùn)用到可靠性分析[2]中。但是這些忽略了一個(gè)重要的問題，即研究都是提前假設(shè)了NHPP，沒有對(duì)其合理性給出證明。因?yàn)镹HPP的強(qiáng)度函數(shù)λ(t)形式變化的多樣性和不確定性，NHPP的驗(yàn)證也成了一個(gè)難題。國外針對(duì)NHPP的假設(shè)檢驗(yàn)有一定程度的研究，Brown[3]提出了用PC NHPP近似NHPP的驗(yàn)證方法，根據(jù)λ(t)對(duì)于時(shí)間變化的連續(xù)性和緩慢性，在一定小的區(qū)間上可以認(rèn)為λ(t)不變，然后在每個(gè)充分小子區(qū)間上驗(yàn)證泊松性。Kim[4]總結(jié)了四種NHPP的檢驗(yàn)方法，對(duì)小區(qū)間的泊松過程做恰當(dāng)?shù)霓D(zhuǎn)換，變成我們易處理的分布類型，然后進(jìn)行K-S檢驗(yàn)，最終驗(yàn)證數(shù)據(jù)的NHPP性。

參數(shù)估計(jì)是在NHPP得到驗(yàn)證后，對(duì)擬合模型的估算參數(shù)[5-6]，求解均值函數(shù)和強(qiáng)度函數(shù)，進(jìn)而為排隊(duì)論模型分析做準(zhǔn)備?？煽啃灾谐Ｓ玫氖亲訁^(qū)間上的分段函數(shù)擬合[7]，但是鑒于簡潔性，本文在整區(qū)間上運(yùn)用最小二乘估計(jì)。

1 常見NHPP模型

1.1 Goel-Okumoto模型

G-O模型是NHPP類軟件可靠性最基本的模型[8]，該模型的適用環(huán)境較理想化，但是模型簡單，便于求解，故而應(yīng)用廣泛。該模型假設(shè)故障被檢測到的概率b是恒定不變的。

強(qiáng)度函數(shù)為：λ(t)=abe-bt，表示的是單位時(shí)間故障輸入率的變化趨勢。

均值函數(shù)為：m(t)=a(1-e-bt)，其表示的實(shí)際意義就是到時(shí)間t為止的累積頻數(shù)，故而本文中用均值函數(shù)來擬合原數(shù)據(jù)的累積頻數(shù)。

1.2 Yamada Delayed S-Shaped模型

Y-D模型是G-O模型的進(jìn)一步推導(dǎo)，該模型適當(dāng)?shù)胤艑捔瞬糠掷碚摷僭O(shè)，更符合實(shí)際情況。模型假設(shè)故障被檢測到的概率b是與時(shí)間有關(guān)的函數(shù)，且檢測率隨著時(shí)間的增加、系統(tǒng)內(nèi)故障的減少而越來越大。

強(qiáng)度函數(shù)為：λ(t)=ab2te-bt，

均值函數(shù)為：m(t)=a[1-(1+bt)e-bt]，

a>0，b>0

1.3 Inflected S-Shaped模型

I-S模型也是假設(shè)故障檢測率b是一個(gè)與時(shí)間有關(guān)的函數(shù)，只是函數(shù)的形式發(fā)生了變化，自變量時(shí)間t與檢測率函數(shù)b(t)不再是多項(xiàng)式關(guān)系，變成了負(fù)指數(shù)形式的變化趨勢。

2 基于CU變換的NHPP統(tǒng)計(jì)檢驗(yàn)?zāi)Ｐ?/h2>

在呼叫中心系統(tǒng)中，一天24小時(shí)的呼入率是不一樣的。因此，在以往的排隊(duì)論分析中，均假設(shè)每小時(shí)的呼入率是常數(shù)，在每個(gè)時(shí)間段上建立排隊(duì)論模型，綜合分析對(duì)應(yīng)指標(biāo)。類似于常用的分段常數(shù)的理念，直接把NHPP的強(qiáng)度函數(shù)轉(zhuǎn)化為PC NHPP進(jìn)行檢驗(yàn)分析。

2.1 CU(conditional uniform)轉(zhuǎn)換K-S檢驗(yàn)

定理1 泊松分布事件發(fā)生時(shí)刻條件均勻原理：在[0，T]時(shí)間段內(nèi)，已知事件發(fā)生了n次的前提條件下，各事件發(fā)生的時(shí)刻{Ti，i∈1∶n}在不考慮順序依存關(guān)系發(fā)生的情況下，可看做相互獨(dú)立的U[0，T]隨機(jī)變量[9]。

本文的CU變換思想就源自泊松分布事件發(fā)生時(shí)刻條件均勻原理。當(dāng)不考慮Ti的順序時(shí)，{Ti，i∈1∶n}～U[0，T]，變形為{Ti/T，i∈1∶n}～U[0，1]，則經(jīng)驗(yàn)累積分布函數(shù)Fn(x)：

(1)

同時(shí)，由均勻分布知對(duì)應(yīng)的理論累積分布函數(shù)F(x)：

F(x)=x，0≤x≤1

2.2 Log(logarithmic)變換K-S檢驗(yàn)

定理2 泊松分布事件發(fā)生的時(shí)間間隔Xn，n=1，2，…服從參數(shù)為λ的指數(shù)分布，且相互獨(dú)立。

由定理2知，欲檢驗(yàn)小區(qū)間上的泊松性質(zhì)，也可通過檢驗(yàn)事件發(fā)生的時(shí)間間隔是不是服從指數(shù)分布來檢驗(yàn)數(shù)據(jù)的泊松性。Brown構(gòu)造了率1指數(shù)分布隨機(jī)變量的Log變換：

1{Xjlog，n≤x}理論累積分布函數(shù)F(x)：

F(x)=1-e-λx，x≥0

2.3 K-S檢驗(yàn)

1）掃描矢量化陜西省民政廳提供的1998年鄉(xiāng)級(jí)行政區(qū)域界線協(xié)議書附圖，獲取全省鄉(xiāng)級(jí)行政區(qū)域界線、界址點(diǎn)、界樁點(diǎn)和三交點(diǎn)等原始界線矢量數(shù)據(jù)，并對(duì)矢量數(shù)據(jù)進(jìn)行坐標(biāo)轉(zhuǎn)換，即1954年北京坐標(biāo)系、1956黃海高程系轉(zhuǎn)換至2000國家大地坐標(biāo)系、1985國家高程基準(zhǔn)。

定義 Kolmogorov分布函數(shù)：

K-S檢驗(yàn)是檢驗(yàn)小樣本數(shù)據(jù)的泊松性。K-S檢驗(yàn)分為單樣本K-S檢驗(yàn)和兩樣本K-S檢驗(yàn)，單樣本是檢驗(yàn)一組樣本數(shù)據(jù)和已知的概率分布類型的擬合優(yōu)度的。

由定理3，統(tǒng)計(jì)量的漸近分布為：

3 實(shí)例驗(yàn)證

3.1 數(shù)據(jù)來源與清洗

Bugzilla失效數(shù)據(jù)是由Mozilla公司開發(fā)的錯(cuò)誤追蹤系統(tǒng)Bugzilla(http://www.bugzilla.org/)在版本升級(jí)和內(nèi)測期間由于系統(tǒng)失效而產(chǎn)生的一系列被詳細(xì)記錄的故障數(shù)據(jù)組成。本文數(shù)據(jù)即2010年11月Bugzilla的第四版本正式內(nèi)測上線開始記錄的4.0版本歷史故障數(shù)據(jù)。統(tǒng)計(jì)每月接收到的故障頻數(shù)。在Bugzilla V4.0中，故障輸入時(shí)間共延續(xù)45個(gè)單位，其中12個(gè)空數(shù)據(jù)，且都集中于序列尾部，可參照?qǐng)D1。這就是統(tǒng)計(jì)中的拖尾現(xiàn)象。本文采取連續(xù)兩次遇零截尾，即Freq連續(xù)兩次為零時(shí)，截去后面的數(shù)據(jù)。

圖1 故障輸入過程頻數(shù)

3.2 NHPP的檢驗(yàn)

由于K-S檢驗(yàn)需要隨機(jī)變量是連續(xù)的，故而針對(duì)單個(gè)到達(dá)的離散變量，可以考慮其到達(dá)時(shí)刻、到達(dá)時(shí)間間隔等連續(xù)型變量，然后采用適合于小樣本檢驗(yàn)的K-S檢驗(yàn)準(zhǔn)確推導(dǎo)結(jié)論。

3.2.1 CU-KS檢驗(yàn)

CU變換基于定理1對(duì)子區(qū)間樣本數(shù)據(jù)進(jìn)行轉(zhuǎn)換。針對(duì)故障輸入時(shí)間數(shù)據(jù)。記到達(dá)時(shí)刻Opened_Day為Ti，區(qū)間長度T為31天，則根據(jù)CU變換理論，此處有Ti/T～U[0，1]。針對(duì)第四個(gè)子區(qū)間的CU變換數(shù)據(jù)整理如表1：

表1 第四個(gè)子區(qū)間CU變換數(shù)據(jù)

運(yùn)用統(tǒng)計(jì)軟件中的R軟件結(jié)合K-S檢驗(yàn)理論檢驗(yàn)數(shù)據(jù)的泊松性質(zhì)，運(yùn)行結(jié)論如下：

>ks.test(Temp，"punif")

>D = 0.12054， p-value = 0.9453

顯然，子區(qū)間上的泊松性是明顯的。但樣本量是19，在原數(shù)據(jù)中是相對(duì)較多的，故而檢驗(yàn)效果相對(duì)較好，也證明了在足夠數(shù)據(jù)量的情況下，小區(qū)間的泊松性質(zhì)驗(yàn)證很好。

但是在表2中有數(shù)據(jù)量Freq為1或0的情況，此時(shí)CU變換后小區(qū)間的泊松性質(zhì)并不樂觀。結(jié)果說明子區(qū)間上小樣本情況下的泊松性質(zhì)相對(duì)于多樣本的表現(xiàn)較差。同時(shí)需要注意，多樣本的小區(qū)間數(shù)據(jù)也出現(xiàn)有泊松檢驗(yàn)效果較差的現(xiàn)象，分析發(fā)現(xiàn)，故障會(huì)有偶然聚集出現(xiàn)的特殊情況，這對(duì)統(tǒng)計(jì)量的均勻性檢驗(yàn)很不利，故而出現(xiàn)了少有的特例情況。

表2 CU變換后各子區(qū)間上的p值和D值

但是整體來講，E(p)=0.522，E(D)=0.397，整體檢驗(yàn)效果還是達(dá)到了的，已經(jīng)滿足了不能拒絕原假設(shè)PC NHPP的條件。

3.2.2 Log-KS檢驗(yàn)

Log變換是基于定理2對(duì)子區(qū)間數(shù)據(jù)進(jìn)行的對(duì)數(shù)變換。第四個(gè)子區(qū)間Log變換后數(shù)據(jù)如表3所示，根據(jù)Log變換理論可知，數(shù)據(jù)表中的序列0，是因?yàn)榍昂髢蓚€(gè)時(shí)間點(diǎn)一致使得真數(shù)為1造成的。無窮小量(Inf)則是因?yàn)檎鏀?shù)趨于0所致。

表3 第四個(gè)子區(qū)間Log變換數(shù)據(jù)

針對(duì)第四子區(qū)間變換數(shù)據(jù)進(jìn)行K-S率1指數(shù)檢驗(yàn)，R程序運(yùn)行結(jié)果如下：

>ks.test(xlog，"pexp")

>D = 0.31579， p-value = 0.04521

顯然的，同一組預(yù)處理的數(shù)據(jù)，Log變換后檢驗(yàn)效果沒有CU檢驗(yàn)效果明顯。這也就暗示我們Log變換的數(shù)據(jù)預(yù)處理更復(fù)雜，需要更多的細(xì)節(jié)調(diào)整。

Log變換后的整體結(jié)論是E(p)=0.378，E(D)=0.327，雖然通過檢驗(yàn)，但Log變換后的數(shù)據(jù)的率1指數(shù)檢驗(yàn)效果沒有CU變換好。這與數(shù)據(jù)選取的時(shí)間不是絕對(duì)連續(xù)，使得對(duì)數(shù)變換對(duì)結(jié)果的影響較大等因素有關(guān)。Log變換時(shí)，應(yīng)該注意盡可能使得每個(gè)樣本數(shù)據(jù)都不一樣，避免0值和無窮小值出現(xiàn)，影響檢驗(yàn)結(jié)果。鑒于驗(yàn)證NHPP性已經(jīng)達(dá)到，此處Log變換的數(shù)據(jù)處理方法不做更深研究。

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

鑒于故障輸入過程N(yùn)HPP性質(zhì)的檢驗(yàn)成立，此處用常用的NHPP類軟件可靠性模型對(duì)數(shù)據(jù)進(jìn)行擬合分析。用Yamada Delayed S-Shaped(Y-D)模型擬合分析：

>nls(X$CFreq ～ a*(1-(1+b*X$time)*exp(-b*X$time))，data=X[，c(1，3)]，start = list(a=200，b=0.5))

用R里的nls函數(shù)對(duì)均值函數(shù)的參數(shù)進(jìn)行估計(jì)，得a=252，b=0.2044，且兩個(gè)參數(shù)是顯著性不等于0。此時(shí)的標(biāo)準(zhǔn)殘差為3.811，迭代次數(shù)為6，容差為7個(gè)分?jǐn)?shù)位。

均值函數(shù)：

m(t)=252·[1-(1+0.2044·t)·e-0.2044t]

圖2 Y-D模型擬合

由圖2可以看出，Y-D模型擬合了故障輸入過程累積頻數(shù)的趨勢。分析K-S擬合檢驗(yàn)結(jié)果，p值為0.9794，接近1，且D值相對(duì)較小，接受Y-D擬合。

>ks.test(X$CFreq，f(X$time))

>D = 0.11111， p-value = 0.9794

4 結(jié)論

本文運(yùn)用CU變換和Log變換對(duì)NHPP檢驗(yàn)，效果顯著，但是針對(duì)數(shù)據(jù)的預(yù)處理問題，本文根據(jù)數(shù)據(jù)類型和R編程難易程度選擇了大區(qū)間以月為時(shí)間單位，小區(qū)間以天為時(shí)間單位的處理方法，故而出現(xiàn)了小區(qū)間中有0值出現(xiàn)的情況，所幸K-S檢驗(yàn)過關(guān)。

CU變換和Log變化是用了泊松分布的兩個(gè)特性延伸推導(dǎo)出的NHPP驗(yàn)證方法，其實(shí)還有直接驗(yàn)證法(頻數(shù)統(tǒng)計(jì)法)，也稱為標(biāo)準(zhǔn)泊松檢驗(yàn)。Lewis[11]也曾提出了比較復(fù)雜的Lewis變換對(duì)PP過程進(jìn)行驗(yàn)證，并給出了推導(dǎo)，所以選擇合適的、有效的數(shù)據(jù)變換方法針對(duì)NHPP的檢驗(yàn)都會(huì)有一定的效益。

[1] Dohi T， Matsuoka T， Osaki S. An Infinite Server Queuing Model for Assessment of the Software Reliability[J]. Electronics and Communications in Japan， 2002， 85(3): 43-51.

[2] Huang C Y， Hung T Y. Software reliability analysis and assessment using queueing models with multiple change-points[J]. Computers & Mathematics with Applications， 2010， 60(7): 2015-2030.

[3] Brown L， Zhao L. Statistical Analysis of a Telephone Call Center: A Queueing-Science Perspective[J]. Journal of the American Statistical Association， 2005， 100(March):36-50.

[4] Kim S H， Whitt W. Choosing arrival process models for service systems: Tests of a nonhomogeneous Poisson process[J]. Naval Research Logistics， 2014， 61(1):66-90.

[5] Massey W A， Parker G A， Whitt W. Estimating the parameters of a nonhomogeneous Poisson process with linear rate[J]. Telecommunication Systems， 1996， 5(4):361-388.

[6] 茆詩松. 高等數(shù)理統(tǒng)計(jì)[M]. 北京：高等教育出版社， 1998.

[7] 徐仁佐，劉蓮君，潘志宏，等. NHPP模型擬合質(zhì)量的改進(jìn)[J]. 自然科學(xué)進(jìn)展:國家重點(diǎn)實(shí)驗(yàn)室通訊， 1991(6):535-542.

[8] 楊劍鋒. 復(fù)雜數(shù)據(jù)下的軟件可靠性分析方法[D].貴陽：貴州大學(xué)，2014.

[9] 張波，張景肖.應(yīng)用隨機(jī)過程[M].北京：清華大學(xué)出版社，2007.

[10] 朱力行. Kolmogorov統(tǒng)計(jì)量的精確分布及其在Bootstrap逼近中的應(yīng)用(英文)[J]. Journal of Mathematical Research with Applications， 1991， 11(2):163-164.

[11] Lewis P A W. Some results on tests for Poisson processes[J]. Biometrika， 1965， 36(52):67-77.

(責(zé)任編輯：曾 晶)

A Test Model of NHPP: Based on CU Transformation

FAN Chaoxia1，ZHAO Ming2，YANG Jianfeng3*

(1.College of Science， Guizhou University， Guiyang 550025，China; 2.Faculty of Technology and Sustainable Development，University of Gavle， Sweden 80176;3. College of Information Engineering， Guizhou Institute of Technology， Guiyang 550003，China)

The test of NHPP and the estimation of the mean function about the BUG input data in opened source software are important sub- problems in queuing theory analysis. In this paper， the test of NHPP was converted to the test of Piecewise- Constant NHPP as the continuity and slowness of the mean function， and finally problem converted to the test of the PP on subintervals firstly， and then combined to verify the NHPP on interval. The verification of PP on subinterval is mainly through the CU transformation of the data reasonablely. Also Bugzilla 's error report data was used to analyze the model， and the result shows acceptable. Finally， the parameters of the BUG input process are calculated using the R language fitting parameters， which lays the foundation for the construction and solution of the following queuing theory model.

NHPP test; CU transformation; K- S test; R language

1000-5269(2016)06-0010-04

10.15958/j.cnki.gdxbzrb.2016.06.03

2016-10-17

貴州省科學(xué)技術(shù)基金計(jì)劃(黔科合J字[2015]2064號(hào))；高層次人才科研啟動(dòng)經(jīng)費(fèi)項(xiàng)目(XJGC20150106)

范朝霞(1990-)，女，在讀碩士，研究方向：應(yīng)用統(tǒng)計(jì)，Email：917855385@qq.com.

*通訊作者： 楊劍鋒，Email：jfyang1@163.com.

O211.6

A