基于兩種統計模型的軟件缺陷預測

馬 由，湯 艷，解 斐

(中國電子科技集團第十五研究所 軟件測評中心,北京 100083)

0 引 言

本文采用兩種模型來實現軟件缺陷數預測，即一般的線性回歸模型和貝葉斯網絡模型，并基于LASSO嵌入式特征選擇算法選擇最優的因子集合。LASSO能夠在訓練模型時同時選擇因子，并通過L1算法避免過擬合問題[1]，最后分析線性回歸和貝葉斯網絡模型在軟件缺陷預測方面的優缺點。

在軟件缺陷數的預測中，因子是離散數據和連續數據結合的，采用線性回歸和貝葉斯網絡進行預測，能夠得出這種混合數據更適應于用哪種模型來預測。

1 相關概念

線性回歸是采用歷史數據來預測未來結果的最直接和有效的方法，一般的現實問題都可以用線性回歸模型來描述和預測，且結果往往很準確。

貝葉斯網絡是一個有向無環圖，圖中的節點代表一個獨立屬性，節點的邊代表兩個獨立屬性之間的依賴關系[2]。如圖1所示是一個貝葉斯網絡示例。

圖1 貝葉斯網絡

貝葉斯網絡結構確定時，其中所有屬性的聯合概率分布作為監督值，可以處理有監督學習的數據仿真及預測；當網絡結構不能事先明確，需要評價函數來學習得到網絡模型，評價函數常見的有AIC、BIC以及極大似然估計，都能夠基于假設近似的得到一個較為準確的模型[3]。當網絡結構確立后，就可以通過已知屬性預測未知屬性，例如通過吉布斯采樣采集已知屬性樣本，在網絡中計算未知屬性的后驗概率。

LASSO是一種常用的嵌入式特征選擇方法，用來解決特征值較多，LASSO方法有嵌入式特征選擇方法的優點，即能夠同時訓練模型和因子選擇，LASSO相較其它嵌入式特征因子,選擇的另一個顯著優勢是通過L1正則化形成了稀疏矩陣[4]，能夠最小化的對相關因子進行降維，且訓練結果較為準確，常用于復雜因子模型構建前的數據預處理。

2 研究過程

2.1 缺陷影響因子分析

缺陷影響因子是在軟件生命周期過程中作為輸入能夠影響輸出結果并引起缺陷的因素。本文所述的因子是從預測對象形成之前的階段中得到的，例如預測需求規格說明的缺陷數，其缺陷影響因子可能是需求分析人員的水平、用戶的表達程度、需求規格說明文檔的規模等，由此可以得知缺陷影響因子的值既有離散數也有連續數。

由于一個軟件缺陷在未確定原因之前，可能有多個影響因子[10]，因此借鑒故障樹分析的思路，采用構建無環有向圖的方式來分析軟件的缺陷影響因子。

因在軟件生命周期過程活動中引入的各種因素具有直接、間接兩種關系，沒有循環關系，因此可以構建有向無環圖。

具體采用兩種構建方式，一種是正向構建，即根據經驗通過原因推導缺陷，一種是逆向構建，通過常用的故障模式來逆向推導原因。構建有向無環圖默認兩個前提，一是缺陷和故障模式是一個抽象的概念，即在圖中只能作為一個節點；二是圖中的原因能夠量化，一般對定量的因子直接采集數據，對定性的因子則采用專家法進行量化。

正向構建時，是由原因推導出可能的結果，構建步驟如下：

從軟件生命周期活動中選取一個活動或一個產品作為起點v0；

分析起點與缺陷相關的能夠量化的因素，每個因素作為有向無環圖的起點w0i(i=1,2,3…n,n為因素個數)；

以v0為輸入，分析其輸出活動或產品，記為v1i(i=1,2,3…m,m為輸出個數)；

分析v1i，得到所有可能與缺陷相關且能夠量化的因素，記為w1j(j=1,2,3…l,l為因素個數)；

分析w0i與w1j，如果存在相關關系則繪出由w0i到w0j的路徑；

循環上述步驟構建有向無環圖V和W，直到目標節點，即缺陷節點。

上述步驟中，缺陷節點的輸入可以是直接影響因子，也可以是間接影響因子，隨應用經驗和分析粒度決定。上述步驟構建的W圖中所有的節點wij即缺陷影響因子，具體W圖如圖2所示。

圖2 缺陷影響因子關系有向無環W圖

逆向構建網絡步驟與正向構建思想相同，思路相反，這里不再贅述。為后續計算方便，圖2用表格的形式表示見表1。

表1 W圖的表格化形式

注：節點表示因子，前驅表示輸入，后繼表示輸出，邊表示與后繼的相關度。

2.2 貝葉斯網絡模型構建過程

將上節所構建的有向無環圖進行模型化處理即形成貝葉斯網絡模型，進而求解最優貝葉斯網絡作為預測模型。

求最優貝葉斯網絡是一個NP難題，常用的算法有貪心算法和約束消減算法，因缺陷預測前不確定影響因子的重要性，采用貪心算法結合LASSO特征選擇可以避免重要因子被剔除，貪心算法需要一個初始的網絡模型[11]。

假設有N個缺陷影響因子，構建網絡為B，因子參數為θ，將網絡B的結構記為

B=

(1)

其中，G為一個有向無環圖，如圖2所示，每一個節點對應一個影響因子，參數θ用來描述因子之間的關系，假設屬性xi在圖G中的父節點為πi，則將每個屬性的條件概率結構記為

θxi/yi=PB(xi/π)

(2)

假設G有n個節點，則網絡G的屬性x1,x2,x3…xn的聯合概率分布定義為

(3)

在式(3)中，πi為第i個類別，訓練集為D，則有

(4)

式中：Dπi,xi為第i個屬性在第i個分類里的個數，Dπi為第i個分類的個數。依據經驗僅僅比例估計條件概率結果較為脆弱，尤其當屬性較多訓練數據較少的情況下，因此將式(4)改為如下形式

(5)

式(5)采用m估計方法，其中m是等價樣本大小的參數，p為用戶指定的參數，一般依據m*p=1對p進行指定。

屬性量化值為連續型數值，則

(6)

(7)

為避免屬性和類別混淆，式(7)中πj為訓練集D中第j個類別，m為類別個數，n為屬性個數，當預測目標是一個值時，m=1。

2.3 模型學習及預測

貝葉斯網絡一般的學習算法是針對初始網絡模型定義一個評分函數，求評分函數的最大或最小值進而得到訓練模型[12]。本文確定訓練模型的思路是將初始貝葉斯網絡模型與LASSO特征選擇方法結合后進行學習，在得到一個系數稀疏解的同時將訓練模型確定下來的過程。

(1)確定模型

按照LASSO特征選擇方法，線性回歸以及概率函數以平方損失函數為評估模型，優化目標隨屬性的離散和連續特征不同而不同[11]。依據經驗，結合上節的初始網絡模型構建優化目標評估函數如下

(8)

式(8)中，θT為矩陣θ的轉置矩陣，λ為正則化參數，在LASSO中規定λ>0。

(2)求解模型參數

求解目標函數是求解使式(8)的值最小的θ的過程，在計算過程中最小二乘法是數學求解過程，但需要x滿秩，而在缺陷預測中并不能保證訓練集的每個因子都有量化值，所以采用梯度下降算法進行求解。

(9)

然后帶入梯度下降算法，梯度下降算法如下：

輸入：數據集D;

特征集A;

評估函數f(θ);//式(8)和式(9),y和x作為常量

步長step;//根據經驗取值，會影響迭代收斂速度

迭代閾值ξ

正則化參數λ;

過程：

初始化變量

f_change=f(θ);//兩次殘差的差值

f_current=f(θ);//當前殘差

循環求解

While D 非空 do

//梯度下降θ的值

θ=θ-step*f(θ′);//依據經驗f(θ′) 為f(θ) 的導數

//計算殘差差值

f_change=f_current-f(θ)；

//計算殘差

f_current=f(θ);

end while

將θ加入到中;

將該步的數據對y,x從D中刪除;

End while

(3)缺陷預測

3 實例及結果分析

本文采用某單位近五年積累的軟件研制過程的數據作為訓練和測試樣本，采集項目類型為辦公自動化項目需求規格說明的缺陷數，以及其它因子。生產這些數據的研發過程是嚴格按照GJB5000A-2008《軍用軟件研制能力成熟度模型》中的規定而制定和執行的，因此這些數據具有普遍性和代表性。具體如圖3所示。

圖3 需求缺陷影響因子集合

由于圖3中的缺陷因子在某單位的數據積累中都是定量的數據，可以直接采集到。

3.1 LASSO特征值選擇

對上述各屬性進行特征值選擇，采用R語言中的glmnet包進行，主要代碼如下：

#lasso特征選擇代碼

mydata<-read.csv(file="C:\Users\user17\Downloads\樣本數據.csv",header=TRUE)

x=as.matrix(mydata[, 1:8])

y=as.matrix(mydata[, 9])

la<-lars(x,y,type="lar")

plot(la)

summary(la)

運行結果如下：

LAR sequence

Computing X′X…

LARS Step 1： Variable 2 added

LARS Step 2： Variable 3 added

LARS Step 3： Variable 4 added

LARS Step 4： Variable 5 added

LARS Step 5： Variable 6 added

LARS Step 6： Variable 1 added

LARS Step 7： Variable 7 added

LARS Step 8： Variable 8 added

LARS Step 9： Variable 9 added

Computing residuals, RSS etc…

LARS/LAR

Call: lars(x=x,y=y, type="lar", trace=10)

Df Rss Cp

0 1 178.15 30.2213

1 2 165.24 16.2195

2 3 160.34 15.7000

3 4 160.11 13.6305

4 5 158.34 9.0796

5 6 156.87 6.6245

6 7 154.23 4.8069

7 8 148.40 8.0000

7 9 132.40 10.0000

依據運行結果，選擇CP值最小的步驟，即選擇屬性x2,x3,x4,x5,x6,x1,x7, 去掉x8,x9。

3.2 線性回歸預測

依據上一節結果嗎，去掉了x中與y相關性弱的節點，采用R語言進行線性回歸，代碼如下：

predata<-read.table("E:\mayou\test.csv",header=FALSE,sep=",")

y<-mydata[,8]

x1<-mydata[,1]

x2<-mydata[,2]

x3<-mydata[,3]

x4<-mydata[,4]

x5<-mydata[,5]

x6<-mydata[,6]

x7<-mydata[,7]

plot<-lm(y～x1+x2+x3+x4+x5+x6+x7,data=predata)

plot

#結果

Call:

lm(formula=y～x1+x2+x3+x4+x5+x6+x7,data=predata)

Coefficients:

(Intercept) x1 x2 x3 x4

42.437 42.4607 0.9408 -4.4549 -2.9723

x5 x6 x7

4.1206 -0.7707 -8.2829

#回歸診斷

plot(plot)

從運行結果可知，最后的模型為：

Y=2.4607x1+0.9408x2-4.4549x3-2.9723x4+4.1206x5-0.7707x6-8.2829x7， 依據上述預測模型，假設置信區間為0.95，將測試集帶入，代碼和結果如下：

predata<-read.table("E:\mayou\test.csv",header=FALSE,sep=",")

plot<-lm(v8～V1+V2+V3+V4+V5+V6+V7,data=predata)

newdata<-read.table("E:\mayou\test1.csv",header=FALSE,sep=",")

predict(plot,data.frame(newdata[,1:7]))

#結果

1 2 3 4

13.6510211 -0.6531836 13.6781614 26.7401970

5 6 7 8

7.7557199 23.2157111 31.3797127 22.2473188

9 10 11

5.0710729 24.7817223 2.4761034

與真實值的對比結果如下：

323027257131092243793326323518202621

誤差均值為2.8。

3.3 貝葉斯網絡預測

根據各因子的潛在關系，構建貝葉斯網絡模型，如圖4所示。

圖4 貝葉斯網絡模型

采集圖4中所有節點的數據，按照3∶1的比例分為訓練數據和測試數據集合。基于訓練數據，用R語言繪制貝葉斯網絡模型[12]，計算得到因子系數θ的值，主要代碼如下：

#安裝R貝葉斯庫

install.packages("D:/bnlearn_4.1.zip", repos = NULL, type = "source")

#引入庫

library(bnlearn)

#導入數據

mydata1<-read.table("D:/Rworkspace/sample.csv",header=FALSE,sep=",")

#離散化缺陷數

c8<-t(v8)

cc8<-sort(c8)

group<-c(seq(1,max(cc8),10),max(cc8))

t8<-cut(c8,breaks=group,labels=1:(length(group)-1))

t8

[1] 3 2 1 1 4 1 2 2 2 4 4 2 1 1 2 6 5 1 3 2 1 1 1 1 3 3 4 6 5 4 7 11 1 2 3

[36] 2 5 5 2 1 4 1 1

netdata<-cbind(predata[,1:7],t8)

#繪制貝葉斯節點

bayesnet<-empty.graph(names(mydata))

bayesnet<-set.arc(bayesnet,"V1","V2")

bayesnet<-set.arc(bayesnet,"V3","V2")

bayesnet<-set.arc(bayesnet,"V1","V3")

bayesnet<-set.arc(bayesnet,"V1","V4")

bayesnet<-set.arc(bayesnet,"V3","V4")

bayesnet<-set.arc(bayesnet,"V6","V4")

bayesnet<-set.arc(bayesnet,"V1","V8")

bayesnet<-set.arc(bayesnet,"V2","V8")

bayesnet<-set.arc(bayesnet,"V3","V8")

bayesnet<-set.arc(bayesnet,"V4","V8")

bayesnet<-set.arc(bayesnet,"V5","V8")

bayesnet<-set.arc(bayesnet,"V6","V8")

bayesnet<-set.arc(bayesnet,"V7","V8")

plot(bayesnet)

fitted<-bn.fit(bayesnet,mydata,method="mle")

fitted

pre<-predict(fitted,data=mydata,node='V8')

#顯示網絡

plot(bayesnet)

#用導入的數據訓練模型

fitted<-bn.fit(bayesnet,mydata1,method="mle")

#預測

newdata<-read.table("E:\mayou\test1.csv",header=FALSE,sep=",")

pre<-predict(fitted,data=newdata,node=′V8′)

運行結果中，需求規格缺陷節點的概率結果為：

Parameters of node V8 (Gaussian distribution)

Conditional density: V8|V1+V2+V3+V4+V5+V6+V7

Coefficients:

(Intercept) V1 V2

4.47062148 0.25553858 0.09860272

V3 V4 V5

-0.44597170 -0.33725753 0.40965399

V6 V7

-0.02004204 -0.81317686

Standard deviation of the residuals: 1.466314

由于缺陷數據是連續值，進行了離散化，因此需要將預測結果的離散值轉換為缺陷數，本文以最大閾值為預測值，與真實結果對比見表2。

表2 貝葉斯網絡預測結果與真實結果對比

貝葉斯網絡預測結果誤差均值為：3.5。

3.4 結果分析

由上述兩個實例結果可知，兩個模型都比較適合預測缺陷數量，但由于缺陷數是連續值，貝葉斯網絡預測模型的準確度要低一些，且計算較復雜。

4 結束語

由于在構建貝葉斯網絡時具有極大的主觀性，因此當數據本身已經是二次加工后的，例如定性數據轉換為定量數據這種情況，則貝葉斯網絡的預測準確度就會明顯降低。因此，當預測樣本數據經過二次加工后形成，且具有主觀性的情況下，采用線性模型預測結果可能更準確一些。

軟件缺陷由于多種原因導致，在實踐中由于各種原因量化的不準確性等問題，理論預測的結果往往不盡如人意，要解決類似的問題，應該從數據積累入手，逐步尋找規律提升預測的準確性。