差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
dm:2012 [2012/07/26 08:33] – [残差をプロットする] wataludm:2012 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1
行 434: 行 434:
 </code> </code>
  
 +=== テスト用コード ===
 +
 +下記のコードを走らせて、エラーが出なければ、
 +
 +  * データの読み込み
 +  * 必要なライブラリのインストール
 +
 +が完了していて、準備は整っていることが確認できる。
 +
 +<code>
 +plot(software.debugging$t, software.debugging$Nt,
 +     type="l",
 +     xlab="Staff Days", ylab="Cumulative Faults")
 +library(gam)
 +software.debugging.gam <- gam(Nt~-1+bs(t),
 +                            data=software.debugging)
 +plot(software.debugging.gam,
 +     xlab="Staff Days", ylab="Cumulative Faults")
 +points(software.debugging$t, software.debugging$Nt,
 +       type="l")
 +n <- dim(software.debugging)[1]
 +software.debugging.diff <- data.frame(t=software.debugging$t[2:n], 
 +                                      t.diff=software.debugging$t[2:n]-software.debugging$t[1:(n-1)],
 +                                      Nt.diff=software.debugging$Nt[2:n]-software.debugging$Nt[1:(n-1)],
 +                                      Ct.diff=software.debugging$Ct[2:n]-software.debugging$Ct[1:(n-1)])
 +software.debugging.diff$dCt <- software.debugging.diff$Ct.diff/software.debugging.diff$t.diff
 +software.debugging.diff$dNt <- software.debugging.diff$Nt.diff/software.debugging.diff$t.diff
 +software.debugging.diff <- rbind(c(0,0,0,0,0,0),software.debugging.diff)
 +lambda.t <- function(theta, alpha, beta, j, data) {
 +  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j]))
 +  diff.2 <- (alpha+beta*sum(data$Ct.diff[1:(j-1)]*exp(theta*data$t[1:(j-1)])))
 +  return(diff.1*diff.2)
 +}
 +neg.log.lik <- function(x) {
 +  theta <- x[1]
 +  alpha <- x[2]
 +  beta  <- x[3]
 +  J <- dim(software.debugging.diff)[1]
 +  log.lik.temp <- 0
 +  for( j in c(2:J) ) {
 +    lambda.j <- lambda.t(theta,alpha,beta,j, software.debugging.diff)
 +    log.lik.temp <- log.lik.temp - lambda.j
 +    log.lik.temp <- log.lik.temp + software.debugging.diff$Nt.diff[j]*log(lambda.j)
 +  }
 +  return(-log.lik.temp)
 +}
 +fitted <- optim(c(0.001, 19, 0.005), neg.log.lik)
 +print(fitted)
 +</code>
 === まずは不具合発見数の成長を眺める === === まずは不具合発見数の成長を眺める ===
  
行 566: 行 615:
 points(software.debugging.diff$t, software.debugging.diff$dCt, pch=20) points(software.debugging.diff$t, software.debugging.diff$dCt, pch=20)
 </code> </code>
 +
  
 === モデルを当てはめる === === モデルを当てはめる ===
  
-  * このデータの場合、データ観測系列は1本なので、ロバストな推定方法は使わない。 +今回は、用いるモデルがRには用意されていないので、
-  * 先週の説明の通り、飽和しかけているので、曲線を強度関数に用いる (累積データに「飽和する曲線」を当てはめる場合、変曲点後まで観測できているか否かで、当てはまりがとても異なるので注意)+
  
-データの先頭に1行、0を入れておく。+  * モデルの強度関数と、強度関数を含む対数尤度関数を、Rの関数として記述する。 
 +  * 対数尤度関数を、optim()を用いて最大化する
  
-<code> +という手順で、パラメータの推定値を得る。
-software.debugging.diff <- rbind(c(0,0,0,0,0,0),software.debugging.diff) +
-</code> +
- +
-強度関数は，次の通り。 +
- +
-<code> +
-lambda.t <- function(theta, alpha, beta, j, data) { +
-  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j])) +
-  diff.2 <- (alpha+beta*sum(data$Ct.diff[1:(j-1)]*exp(theta*data$t[1:(j-1)]))) +
-  return(diff.1*diff.2) +
-} +
-#lambda.t(0.001704,19.1745,0.003045,10,software.debugging.diff) +
-</code> +
- +
-これを用いた対数尤度関数を、次のように定義する。 +
-<code> +
-neg.log.lik <- function(x) { +
-  theta <- x[1] +
-  alpha <- x[2] +
-  beta  <- x[3] +
-  J <- dim(software.debugging.diff)[1] +
-  log.lik.temp <- 0 +
-  for( j in c(2:J) ) { +
-    lambda.j <- lambda.t(theta,alpha,beta,j, software.debugging.diff) +
-    log.lik.temp <- log.lik.temp - lambda.j +
-    log.lik.temp <- log.lik.temp + software.debugging.diff$Nt.diff[j]*log(lambda.j) +
-  } +
-  return(-log.lik.temp) +
-} +
-</code> +
- +
-対数尤度の最大化は、上の関数の最小化と等しい。 +
- +
-Rには最小化の関数が幾つかある。nlminb(), optim(), nlm()など。 +
- +
-これを最小化するのは、結構、困難であるが、例えば初期値を +
- +
-<code> +
-fitted <- optim(c(0.001, 19, 0.005), neg.log.lik) +
-print(fitted) +
-</code> +
- +
-のように与えると、それなりの値に収束する。optim()の最初の引数は、最適化の初期値である。これは Cook and Lawless (2007)掲載の推定値を参考に設定した。 +
-以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。 +
- +
-=== モデルを当てはめる ===+
  
   * このデータの場合、データ観測系列は1本なので、ロバストな推定方法は使わない。   * このデータの場合、データ観測系列は1本なので、ロバストな推定方法は使わない。
行 680: 行 684:
  
 === 強度関数を変えてみる === === 強度関数を変えてみる ===
 +
 +ここでは、強度関数の変え方の一例を示す。
 +
 +  * 強度関数を変えるにはlambda.t()だけでなく、neg.log.lik()も変える必要が出てくることもある。今回は、パラメータを追加したので、neg.log.lik()の中で x[4] の扱いを追加し、4変数関数としなければならない。
 +  * さらに閾値を変えてみるとして、betaとbeta.2の切り替え点を300から500に変えるなど、してみている。
  
 過去の修正の影響がもっと軽い可能性はないかと、betaを2種類にしてみる。 過去の修正の影響がもっと軽い可能性はないかと、betaを2種類にしてみる。
 +まずは300日以内の影響をbetaとし、300日以降の影響をbeta.2とするように、強度関数を変更する。
  
 <code> <code>
行 695: 行 705:
 lambda.t(0.001704,19.1745,0.003045,0.003045,10,software.debugging.diff) lambda.t(0.001704,19.1745,0.003045,0.003045,10,software.debugging.diff)
 </code> </code>
 +
 +強度関数のパラメータがひとつ増えたので、対数尤度関数の引数xの次元もひとつ増える。
 +そのため、一行追加する。
  
 <code> <code>
行 715: 行 728:
 直近の改修よりも、以前の改修の方が、不具合の発生に大きく影響を与える、という母数の推定値が得られる。 直近の改修よりも、以前の改修の方が、不具合の発生に大きく影響を与える、という母数の推定値が得られる。
 これは、初期に多くの改修を行ったことと、無関係ではないかもしれない。 これは、初期に多くの改修を行ったことと、無関係ではないかもしれない。
 +
 +以下はbetaとbeta.2の切り替え時点を300から500に変更してみる場合の強度関数である。
  
 <code> <code>
行 727: 行 742:
 } }
 </code> </code>
 +
 +このような変更手段だけでも、かなりのバリエーションにはなるので、いろいろ試してみて欲しい。
  
 === モデルを選ぶ === === モデルを選ぶ ===
  
-AICは、モデルの予測精度を評価する指標として導出されており、次のように計算できる。+AICは、モデルの当てはまりの良さを評価する指標のひとつである。 
 + 
 +  * AICは予測精度の近似として導出されている。 
 +  * AICを基準にモデルを選ぶ、とはこの値が最小となるように、モデルを選ぶことである。 
 +  * AICを最小化するモデルが、必ずしも正解とは限らない。 
 + 
 +今回は、次のように計算される量をAICとして用いる。
  
 <code> <code>
行 737: 行 760:
 </code> </code>
  
-AICを基準にモデルを選ぶ、とはこの値が最小となるように、モデルを選ぶことである。+これをなるべく小さくするようなモデルを用いるのが望ましいが、対象についての知識・情報がある場合には、それに基づいて、最小ではないモデルを選択することも少なくない。
  
 === 残差をプロットする === === 残差をプロットする ===
行 811: 行 834:
 他にないか、頑張ってみて。 他にないか、頑張ってみて。
  
-=== テスト用コード === 
-<code> 
-software.debugging.diff <- data.frame(t=software.debugging$t[2:n],  
-                                      t.diff=software.debugging$t[2:n]-software.debugging$t[1:(n-1)], 
-                                      Nt.diff=software.debugging$Nt[2:n]-software.debugging$Nt[1:(n-1)], 
-                                      Ct.diff=software.debugging$Ct[2:n]-software.debugging$Ct[1:(n-1)]) 
-software.debugging.diff$dCt <- software.debugging.diff$Ct.diff/software.debugging.diff$t.diff 
-software.debugging.diff$dNt <- software.debugging.diff$Nt.diff/software.debugging.diff$t.diff 
-software.debugging.diff <- rbind(c(0,0,0,0,0,0),software.debugging.diff) 
-lambda.t <- function(theta, alpha, beta, j, data) { 
-  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j])) 
-  diff.2 <- (alpha+beta*sum(data$Ct.diff[1:(j-1)]*exp(theta*data$t[1:(j-1)]))) 
-  return(diff.1*diff.2) 
-} 
-neg.log.lik <- function(x) { 
-  theta <- x[1] 
-  alpha <- x[2] 
-  beta  <- x[3] 
-  J <- dim(software.debugging.diff)[1] 
-  log.lik.temp <- 0 
-  for( j in c(2:J) ) { 
-    lambda.j <- lambda.t(theta,alpha,beta,j, software.debugging.diff) 
-    log.lik.temp <- log.lik.temp - lambda.j 
-    log.lik.temp <- log.lik.temp + software.debugging.diff$Nt.diff[j]*log(lambda.j) 
-  } 
-  return(-log.lik.temp) 
-} 
-fitted <- optim(c(0.001, 19, 0.005), neg.log.lik) 
-print(fitted) 
-</code>