基于两种统计模型的软件缺陷预测
2020-04-24马由，汤艳，解斐
马 由，汤 艳，解 斐
(中国电子科技集团第十五研究所 软件测评中心,北京 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 结束语
由于在构建贝叶斯网络时具有极大的主观性，因此当数据本身已经是二次加工后的，例如定性数据转换为定量数据这种情况，则贝叶斯网络的预测准确度就会明显降低。因此，当预测样本数据经过二次加工后形成，且具有主观性的情况下，采用线性模型预测结果可能更准确一些。
软件缺陷由于多种原因导致，在实践中由于各种原因量化的不准确性等问题，理论预测的结果往往不尽如人意，要解决类似的问题，应该从数据积累入手，逐步寻找规律提升预测的准确性。
