「データ指向のソフトウェア品質マネジメント」の「テスト工程での不具合検出数の予測」をPythonで行う

 

「データ指向のソフトウェア品質マネジメント」(野中・小池・小室)の4.2節「テスト工程での不具合検出数の予測」をPythonで追ってみることにする。

今年のソフトウェア品質シンポジウムで著者の一人によるチュートリアルを受講したのだが、それがとても面白かったので購入した。

In [1]:
import math
import pandas as pd
from pandas import DataFrame, Series
import matplotlib.pyplot as plt
%matplotlib inline
 

まずは対象データを読み込む。これは書籍のサポートページからダウンロードできる。 ダウンロードに書籍記載のパスワードが必要なので生データをここには載せないが、 データ自体はFenton (2008)のものである。

ソースコード行数KLOCと不具合検出数Failures以外はカテゴリーデータ(順序尺度)で、D1は開発メンバーの経験、P5はステークホルダーの関与の度合い、P6は顧客とのコミュニケーション、P7は協力会社のマネジメントを表す。

In [2]:
data1 = pd.read_excel(u"4.2節/4-2-data.xls", index_col='Prj')
data1.columns
Out[2]:
Index([u'Failures', u'KLOC', u'D1', u'P5', u'P6', u'P7'], dtype='object')
 

ソースコード行数と不具合検出数のヒストグラムを描くと、いずれも右に裾が長い分布となっている。 また、ソースコード行数と不具合検出数には強い正の相関がある。

In [3]:
def scatter_corr_matrix(df):
    pd.scatter_matrix(df)
    plt.show()
    return df.corr()

scatter_corr_matrix(data1)
 
Out[3]:
 FailuresKLOC
Failures 1.000000 0.883083
KLOC 0.883083 1.000000
 

単回帰分析

まずはソースコード行数だけから不具合検出数を予測する単回帰分析を行ってみる。

「データ志向の~」では目的変数と説明変数が右に歪んだ分布となっている場合には、 それぞれの対数をとったうえで両変数を対数変換したモデル式を用いるのが典型的であるとして、 最初から対数変換した単回帰分析を行っているのだが、 ここでは対数変換しない場合とした場合で比べてみることにする。

In [4]:
import statsmodels.formula.api as smf

smf.ols(formula="Failures ~ KLOC", data=data1).fit().summary()
Out[4]:
OLS Regression Results
Dep. Variable: Failures R-squared: 0.780
Model: OLS Adj. R-squared: 0.772
Method: Least Squares F-statistic: 95.64
Date: Sat, 12 Nov 2016 Prob (F-statistic): 2.29e-10
Time: 23:10:02 Log-Likelihood: -200.39
No. Observations: 29 AIC: 404.8
Df Residuals: 27 BIC: 407.5
Df Model: 1    
Covariance Type: nonrobust    
  coefstd errtP>|t|[95.0% Conf. Int.]
Intercept 9.4969 65.505 0.145 0.886 -124.909 143.903
KLOC 11.4165 1.167 9.779 0.000 9.021 13.812
Omnibus: 0.160 Durbin-Watson: 2.363
Prob(Omnibus): 0.923 Jarque-Bera (JB): 0.091
Skew: -0.109 Prob(JB): 0.956
Kurtosis: 2.834 Cond. No. 78.8
In [5]:
smf.ols(formula="np.log(Failures) ~ np.log(KLOC)", data=data1).fit().summary()
Out[5]:
OLS Regression Results
Dep. Variable: np.log(Failures) R-squared: 0.779
Model: OLS Adj. R-squared: 0.771
Method: Least Squares F-statistic: 95.38
Date: Sat, 12 Nov 2016 Prob (F-statistic): 2.35e-10
Time: 23:10:02 Log-Likelihood: -29.481
No. Observations: 29 AIC: 62.96
Df Residuals: 27 BIC: 65.70
Df Model: 1    
Covariance Type: nonrobust    
  coefstd errtP>|t|[95.0% Conf. Int.]
Intercept 2.4606 0.325 7.578 0.000 1.794 3.127
np.log(KLOC) 0.9564 0.098 9.766 0.000 0.755 1.157
Omnibus: 5.037 Durbin-Watson: 2.194
Prob(Omnibus): 0.081 Jarque-Bera (JB): 1.856
Skew: -0.181 Prob(JB): 0.395
Kurtosis: 1.815 Cond. No. 9.02
 

Adj. R-squaredはあまり変わらないのだが、AICで見ると対数変換により改善している。

重回帰分析のための説明変数の検討

次に順序尺度を説明変数に加える検討をする。まずは順位尺度の列に順位を表す整数が入っているdata1-1のシートを読む。

In [6]:
data1_1 = pd.read_excel(u"4.2節/4-2-data.xls", index_col='Prj', sheetname='data1-1')
data1_1.describe()
 
C:\Users\ether\Anaconda2\lib\site-packages\numpy\lib\function_base.py:3834: RuntimeWarning: Invalid value encountered in percentile
  RuntimeWarning)
Out[6]:
 FailuresKLOCD1P5P6P7
count 31.000000 29.000000 31.000000 31.000000 31.000000 7.000000
mean 445.225806 39.379310 3.322581 4.032258 4.064516 2.857143
std 510.915695 40.680299 0.979357 0.752058 0.771815 0.899735
min 5.000000 1.000000 1.000000 3.000000 3.000000 2.000000
25% 100.000000 NaN 3.000000 3.500000 3.500000 NaN
50% 209.000000 NaN 4.000000 4.000000 4.000000 NaN
75% 662.500000 NaN 4.000000 5.000000 5.000000 NaN
max 1906.000000 155.000000 5.000000 5.000000 5.000000 4.000000
 

ソースコード行数がすでに説明変数に含まれているため、 不具合検出密度と各順序尺度の相関(スピアマンの順位相関係数)を計算する。

In [7]:
from scipy.stats import spearmanr
data1_1_d = data1_1.assign(Density = data1_1['Failures'] / data1_1['KLOC'])
data1_1_d.corr(method='spearman')
Out[7]:
 FailuresKLOCD1P5P6P7Density
Failures 1.000000 0.850419 -0.408264 -0.254397 -0.136753 0.018898 0.395567
KLOC 0.850419 1.000000 -0.114229 -0.055048 -0.031084 0.973329 -0.089946
D1 -0.408264 -0.114229 1.000000 0.488971 0.516152 0.632785 -0.692716
P5 -0.254397 -0.055048 0.488971 1.000000 0.645129 0.083666 -0.435407
P6 -0.136753 -0.031084 0.516152 0.645129 1.000000 -0.591960 -0.373876
P7 0.018898 0.973329 0.632785 0.083666 -0.591960 1.000000 -0.790569
Density 0.395567 -0.089946 -0.692716 -0.435407 -0.373876 -0.790569 1.000000
 

不具合検出密度と各順序尺度の関係を箱ひげ図で表す。

In [8]:
data1_1_d.boxplot(column='Density', by='P7', return_type='dict')
data1_1_d.boxplot(column='Density', by='D1', return_type='dict')
data1_1_d.boxplot(column='Density', by='P5', return_type='dict')
data1_1_d.boxplot(column='Density', by='P6', return_type='dict')
plt.show()
 
 
 
 
 

これを見ると不具合検出密度と相関が高いのはP7(協力会社のマネジメントが高いと不具合検出密度が低い)だが、「データ志向の~」ではP7はデータ件数が少ない(n=7)こと、ソースコード行数との相関係数も極めて大きいこと(多重共線性の疑い)を理由に退けている。

 

重回帰分析の説明変数をステップワイズで選択

P7以外の順序尺度をステップワイズで選択する。まずは全部を使う。

In [9]:
smf.ols(formula="np.log(Failures) ~ np.log(KLOC) + D1 + P5 + P6", data=data1).fit().summary()
Out[9]:
OLS Regression Results
Dep. Variable: np.log(Failures) R-squared: 0.897
Model: OLS Adj. R-squared: 0.847
Method: Least Squares F-statistic: 18.29
Date: Sat, 12 Nov 2016 Prob (F-statistic): 1.47e-07
Time: 23:10:03 Log-Likelihood: -18.504
No. Observations: 29 AIC: 57.01
Df Residuals: 19 BIC: 70.68
Df Model: 9    
Covariance Type: nonrobust    
  coefstd errtP>|t|[95.0% Conf. Int.]
Intercept 3.6092 0.788 4.582 0.000 1.961 5.258
D1[T.2_L] -0.0876 0.664 -0.132 0.896 -1.478 1.302
D1[T.3_M] -0.5328 0.712 -0.748 0.463 -2.023 0.957
D1[T.4_H] -1.1245 0.665 -1.690 0.107 -2.517 0.268
D1[T.5_VH] -1.2480 0.775 -1.611 0.124 -2.870 0.373
P5[T.4_H] -0.3889 0.430 -0.904 0.377 -1.290 0.512
P5[T.5_VH] 0.0721 0.682 0.106 0.917 -1.356 1.500
P6[T.4_H] 0.0671 0.338 0.198 0.845 -0.641 0.775
P6[T.5_VH] -0.3956 0.667 -0.593 0.560 -1.793 1.001
np.log(KLOC) 0.9235 0.089 10.430 0.000 0.738 1.109
Omnibus: 0.572 Durbin-Watson: 2.426
Prob(Omnibus): 0.751 Jarque-Bera (JB): 0.560
Skew: 0.296 Prob(JB): 0.756
Kurtosis: 2.664 Cond. No. 49.1
 

ここから1つずつ減らしていく。

In [10]:
formulas = [
    "np.log(Failures) ~ np.log(KLOC) + D1 + P5 + P6",
    "np.log(Failures) ~ np.log(KLOC) + D1 + P5",
    "np.log(Failures) ~ np.log(KLOC) + D1",
    "np.log(Failures) ~ np.log(KLOC)",
]
steps = [[res.llf, res.aic, res.rsquared, res.rsquared_adj, formula]
         for (res, formula) in [(smf.ols(formula=formula, data=data1).fit(), formula) for formula in formulas]]
DataFrame(steps, columns=['Log-Likelihood', 'AIC', 'R-squared', 'Adj. R-squared', 'formula'])
Out[10]:
 Log-LikelihoodAICR-squaredAdj. R-squaredformula
0 -18.504129 57.008257 0.896515 0.847496 np.log(Failures) ~ np.log(KLOC) + D1 + P5 + P6
1 -18.910947 53.821894 0.893571 0.858094 np.log(Failures) ~ np.log(KLOC) + D1 + P5
2 -19.531827 51.063655 0.888914 0.864765 np.log(Failures) ~ np.log(KLOC) + D1
3 -29.480676 62.961352 0.779383 0.771212 np.log(Failures) ~ np.log(KLOC)
 

AICおよびAdj. R-squaredで判断するとlog(KLOC)とD1を使うモデルが最も良い。

In [11]:
smf.ols(formula="np.log(Failures) ~ np.log(KLOC) + D1", data=data1).fit().summary()
Out[11]:
OLS Regression Results
Dep. Variable: np.log(Failures) R-squared: 0.889
Model: OLS Adj. R-squared: 0.865
Method: Least Squares F-statistic: 36.81
Date: Sat, 12 Nov 2016 Prob (F-statistic): 3.09e-10
Time: 23:10:03 Log-Likelihood: -19.532
No. Observations: 29 AIC: 51.06
Df Residuals: 23 BIC: 59.27
Df Model: 5    
Covariance Type: nonrobust    
  coefstd errtP>|t|[95.0% Conf. Int.]
Intercept 3.1819 0.622 5.117 0.000 1.895 4.468
D1[T.2_L] 0.0438 0.613 0.071 0.944 -1.225 1.312
D1[T.3_M] -0.2468 0.565 -0.436 0.667 -1.417 0.923
D1[T.4_H] -1.0533 0.557 -1.890 0.071 -2.206 0.100
D1[T.5_VH] -1.1696 0.660 -1.771 0.090 -2.536 0.197
np.log(KLOC) 0.9333 0.082 11.389 0.000 0.764 1.103
Omnibus: 0.773 Durbin-Watson: 2.459
Prob(Omnibus): 0.679 Jarque-Bera (JB): 0.831
Skew: 0.287 Prob(JB): 0.660
Kurtosis: 2.402 Cond. No. 44.1
 

このサマリーを見るとD1のダミー変数の中でも2_Lと3_MのP値が高い。そこで「データ指向の~」では3_M以下の水準を統合することを提案している。

In [12]:
d1_map = {'1_VL': '3_M', '2_L': '3_M', '3_M': '3_M', '4_H': '4_H', '5_VH': '5_VH'}
data1_d1_new = data1.assign(D1_new = data1['D1'].apply(lambda x: d1_map[x]))

ols_final = smf.ols(formula="np.log(Failures) ~ np.log(KLOC) + D1_new", data=data1_d1_new).fit()
ols_final.summary()
Out[12]:
OLS Regression Results
Dep. Variable: np.log(Failures) R-squared: 0.885
Model: OLS Adj. R-squared: 0.871
Method: Least Squares F-statistic: 64.22
Date: Sat, 12 Nov 2016 Prob (F-statistic): 6.93e-12
Time: 23:10:03 Log-Likelihood: -20.015
No. Observations: 29 AIC: 48.03
Df Residuals: 25 BIC: 53.50
Df Model: 3    
Covariance Type: nonrobust    
  coefstd errtP>|t|[95.0% Conf. Int.]
Intercept 3.1155 0.279 11.161 0.000 2.541 3.690
D1_new[T.4_H] -0.9215 0.201 -4.576 0.000 -1.336 -0.507
D1_new[T.5_VH] -1.0436 0.397 -2.630 0.014 -1.861 -0.226
np.log(KLOC) 0.9110 0.074 12.295 0.000 0.758 1.064
Omnibus: 0.826 Durbin-Watson: 2.461
Prob(Omnibus): 0.662 Jarque-Bera (JB): 0.841
Skew: 0.241 Prob(JB): 0.657
Kurtosis: 2.318 Cond. No. 15.0
 

データ加工により、Adj. R-squaredもAICも向上したことがわかる。

モデルの評価

予測モデルの評価として、「データ指向の~」ではまず予測値と実測値の相関係数を評価している。

In [13]:
scatter_corr_matrix(
    pd.concat(
        [data1.dropna(subset=['KLOC']).reset_index(drop=True),
         Series(np.exp(ols_final.predict()), name='predict')],
        axis=1)[['Failures', 'predict']])
 
Out[13]:
 Failurespredict
Failures 1.000000 0.886818
predict 0.886818 1.000000
 

上記では散布図も表示している。なお、「データ指向の~」(2015年1月13日第4刷)に掲載の散布図はX軸とY軸が逆のように思える(軸ラベル上Xが予測でYが実績なのだが、右上でなく左下の散布図になっている)。

相関係数が0.887であることから、うまく予測できているように見えると判断している。

このくだりにはちょっとよくわからないことろがあって、そもそも予測値と実績値の相関係数というのは予測モデルの評価の仕方として適当なのだろうかという疑問がある。 というのはモデルが誤っていて「本来の値の約半分の値を予測する」というものであった場合、予測モデルとしては役に立たないが、相関係数は半分でない場合と同様に高くなるはずだからだ。

次に予測誤差 MRE (Megnitude of Relative Error)と呼ばれる評価指標を出している。こちらのほうが直感的な定義に思える。MREは1つ1つの予測に対する値なので、全体で1つの指標値にするには予測誤差の平均とか中央値とかで見るようだ。

In [14]:
def mre(actual, predicted):
    df = DataFrame({'actual': actual, 'predicted': predicted})
    mre = df.eval("abs(actual - predicted) / actual")
    mre.plot.box()
    plt.show()
    print mre.describe().round(2)
    return mre.median()

mre(
    data1.dropna(subset=['KLOC'])['Failures'],
    np.exp(ols_final.predict()))
 
 
count    29.00
mean      0.44
std       0.30
min       0.07
25%       0.22
50%       0.34
75%       0.55
max       1.41
dtype: float64
Out[14]:
0.34027626468821626
 

中央値は34%程度(本文では30%と書いているが)である。

ポワソン回帰分析

この節のコラム記事では対数変換したうえでの重回帰分析(最小二乗法)ではなくポワソン回帰を行ったほうがよい予測モデルになるということが書かれている。本文には結果のみで、過程は省かれているが、これを行ってみたい。

Statsmodelsでポワソン回帰を行うにはglmを使う。コラムに掲載された結果と同じになるモデル式は下記だ。

In [16]:
import statsmodels.api as sm

poisson_model = smf.glm(
    formula='Failures ~ np.log(KLOC) + D1 + P5 + P6',
    data=data1,
    family=sm.families.Poisson()).fit()
poisson_model.summary()
Out[16]:
Generalized Linear Model Regression Results
Dep. Variable: Failures No. Observations: 29
Model: GLM Df Residuals: 19
Model Family: Poisson Df Model: 9
Link Function: log Scale: 1.0
Method: IRLS Log-Likelihood: -941.94
Date: Sat, 12 Nov 2016 Deviance: 1674.7
Time: 23:11:18 Pearson chi2: 1.59e+03
No. Iterations: 9    
  coefstd errzP>|z|[95.0% Conf. Int.]
Intercept 3.9355 0.064 61.320 0.000 3.810 4.061
D1[T.2_L] -0.3720 0.045 -8.224 0.000 -0.461 -0.283
D1[T.3_M] -0.7369 0.048 -15.207 0.000 -0.832 -0.642
D1[T.4_H] -0.9145 0.044 -20.781 0.000 -1.001 -0.828
D1[T.5_VH] -1.2531 0.088 -14.187 0.000 -1.426 -1.080
P5[T.4_H] -0.7771 0.035 -21.923 0.000 -0.847 -0.708
P5[T.5_VH] 0.1282 0.086 1.498 0.134 -0.040 0.296
P6[T.4_H] -0.1638 0.028 -5.910 0.000 -0.218 -0.109
P6[T.5_VH] -0.5596 0.083 -6.769 0.000 -0.722 -0.398
np.log(KLOC) 0.9393 0.013 70.423 0.000 0.913 0.965
In [17]:
scatter_corr_matrix(
    pd.concat(
        [data1.dropna(subset=['KLOC']).reset_index(drop=True),
         Series(poisson_model.predict(), name='predict')],
        axis=1)[['Failures', 'predict']])
 
Out[17]:
 Failurespredict
Failures 1.000000 0.950359
predict 0.950359 1.000000
 

予測値と実測値の相関係数は向上している。

一方MREを計算すると中央値では改善するが平均値では悪くなっており、またばらつきも大きくなっている。

In [18]:
mre(
    data1.dropna(subset=['KLOC'])['Failures'],
    poisson_model.predict())
 
 
count    29.00
mean      0.49
std       0.58
min       0.00
25%       0.12
50%       0.25
75%       0.52
max       1.83
dtype: float64
Out[18]:
0.25178105769280135
 

こうしてみてくると対数変換の有無、説明変数の選択、順序尺度の統合、回帰手法の選択といった多くの「ノブ」があって、評価尺度も1つではない。 モデリングが探索的な過程であると感じられる。

In [ ]: