差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
dm:2011 [2011/07/21 08:50] – [今回の解析の流れ] wataludm:2011 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1
行 282: 行 282:
  
 ==== 演習2 ==== ==== 演習2 ====
 +=== 日程 ===
 +|出題|2011.07.21|
 +|〆切|2011.08.04|
 +
 +
 === 今回の解析の流れ === === 今回の解析の流れ ===
  
行 361: 行 366:
 今回は二つのパッケージを用いるので、インストールしておく。 今回は二つのパッケージを用いるので、インストールしておく。
 <code> <code>
 +install.packages(pkgs=c("gam"), 
 +                 repos='http://cran.md.tsukuba.ac.jp/',
 +                 dependencies=TRUE)
 library(gam) library(gam)
 </code> </code>
行 369: 行 377:
  
 <code> <code>
-plot(software.testing$t, software.testing$Nt,+plot(software.debugging$t, software.debugging$Nt,
      type="l",      type="l",
      xlab="Staff Days", ylab="Cumulative Faults")      xlab="Staff Days", ylab="Cumulative Faults")
行 378: 行 386:
 <code> <code>
 library(gam) library(gam)
-software.testing.gam <- gam(Nt~-1+bs(t), +software.debugging.gam <- gam(Nt~-1+bs(t), 
-                            data=software.testing) +                            data=software.debugging) 
-plot(software.testing.gam,+plot(software.debugging.gam,
      xlab="Staff Days", ylab="Cumulative Faults")      xlab="Staff Days", ylab="Cumulative Faults")
-points(software.testing$t, software.testing$Nt,+points(software.debugging$t, software.debugging$Nt,
        type="l")        type="l")
 </code> </code>
行 391: 行 399:
  
 <code> <code>
-n <- dim(software.testing)[1] +n <- dim(software.debugging)[1] 
-software.testing.diff <- data.frame(t=software.testing$t[2:n],  +software.debugging.diff <- data.frame(t=software.debugging$t[2:n],  
-                                    t.diff=software.testing$t[2:n]-software.testing$t[1:(n-1)], +                                    t.diff=software.debugging$t[2:n]-software.debugging$t[1:(n-1)], 
-                                    Nt.diff=software.testing$Nt[2:n]-software.testing$Nt[1:(n-1)], +                                    Nt.diff=software.debugging$Nt[2:n]-software.debugging$Nt[1:(n-1)], 
-                                    Ct.diff=software.testing$Ct[2:n]-software.testing$Ct[1:(n-1)]) +                                    Ct.diff=software.debugging$Ct[2:n]-software.debugging$Ct[1:(n-1)]) 
-software.testing.diff$dCt <- software.testing.diff$Ct.diff/software.testing.diff$t.diff +software.debugging.diff$dCt <- software.debugging.diff$Ct.diff/software.debugging.diff$t.diff 
-software.testing.diff$dNt <- software.testing.diff$Nt.diff/software.testing.diff$t.diff+software.debugging.diff$dNt <- software.debugging.diff$Nt.diff/software.debugging.diff$t.diff
 </code> </code>
  
行 403: 行 411:
  
 <code> <code>
-software.testing.gam <- gam(Nt.diff~-1+bs(t), +software.debugging.gam <- gam(Nt.diff~-1+bs(t), 
-                            data=software.testing.diff) +                            data=software.debugging.diff) 
-plot(software.testing.gam,+plot(software.debugging.gam,
      xlab="Staff Days", ylab="1st Order Diff of Cumulative Faults")      xlab="Staff Days", ylab="1st Order Diff of Cumulative Faults")
-points(software.testing.diff$t,software.testing.diff$Nt.diff,+points(software.debugging.diff$t,software.debugging.diff$Nt.diff,
      pch=20)      pch=20)
 </code> </code>
行 415: 行 423:
 二階差分も確認してみる。 二階差分も確認してみる。
 <code> <code>
-n <- dim(software.testing.diff)[1] +n <- dim(software.debugging.diff)[1] 
-software.testing.diff.2 <- data.frame(t=software.testing.diff$t[2:n],  +software.debugging.diff.2 <- data.frame(t=software.debugging.diff$t[2:n],  
-                                    t.diff=software.testing.diff$t[2:n]-software.testing.diff$t[1:(n-1)], +                                    t.diff=software.debugging.diff$t[2:n]-software.debugging.diff$t[1:(n-1)], 
-                                    Nt.diff.2=software.testing.diff$Nt[2:n]-software.testing.diff$Nt[1:(n-1)], +                                    Nt.diff.2=software.debugging.diff$Nt[2:n]-software.debugging.diff$Nt[1:(n-1)], 
-                                    Ct.diff.2=software.testing.diff$Ct[2:n]-software.testing.diff$Ct[1:(n-1)]) +                                    Ct.diff.2=software.debugging.diff$Ct[2:n]-software.debugging.diff$Ct[1:(n-1)]) 
-software.testing.diff.2$dCt <- software.testing.diff.2$Ct.diff.2/software.testing.diff.2$t +software.debugging.diff.2$dCt <- software.debugging.diff.2$Ct.diff.2/software.debugging.diff.2$t 
-software.testing.diff.2$dNt <- software.testing.diff.2$Nt.diff.2/software.testing.diff.2$t+software.debugging.diff.2$dNt <- software.debugging.diff.2$Nt.diff.2/software.debugging.diff.2$t
 </code> </code>
  
行 427: 行 435:
  
 <code> <code>
-software.testing.gam <- gam(Nt.diff.2~-1+bs(t), +software.debugging.gam <- gam(Nt.diff.2~-1+bs(t), 
-                            data=software.testing.diff.2) +                            data=software.debugging.diff.2) 
-plot(software.testing.gam,+plot(software.debugging.gam,
      xlab="Staff Days", ylab="2nd Order Diff of Cumulative Faults")      xlab="Staff Days", ylab="2nd Order Diff of Cumulative Faults")
-points(software.testing.diff.2$t,software.testing.diff.2$Nt.diff.2,+points(software.debugging.diff.2$t,software.debugging.diff.2$Nt.diff.2,
      pch=20)      pch=20)
 </code> </code>
行 438: 行 446:
  
 <code> <code>
-plot(software.testing.diff.2$t,software.testing.diff.2$Nt.diff.2,+plot(software.debugging.diff.2$t,software.debugging.diff.2$Nt.diff.2,
      pch=20)      pch=20)
 </code> </code>
行 447: 行 455:
  
 <code> <code>
-software.testing.gam <- gam(dNt~-1+bs(t), +software.debugging.gam <- gam(dNt~-1+bs(t), 
-                            data=software.testing.diff.2) +                            data=software.debugging.diff.2) 
-plot(software.testing.gam,+plot(software.debugging.gam,
      xlab="Staff Days", ylab="2nd Order Diff of Cumulative Faults per Days")      xlab="Staff Days", ylab="2nd Order Diff of Cumulative Faults per Days")
-points(software.testing.diff.2$t,software.testing.diff.2$dNt,+points(software.debugging.diff.2$t,software.debugging.diff.2$dNt,
      pch=20)      pch=20)
 </code> </code>
行 458: 行 466:
  
 <code> <code>
-plot(software.testing.diff.2$t[software.testing.diff.2$dNt>=0], +plot(software.debugging.diff.2$t[software.debugging.diff.2$dNt>=0], 
-     software.testing.diff.2$dNt[software.testing.diff.2$dNt>=0], +     software.debugging.diff.2$dNt[software.debugging.diff.2$dNt>=0], 
-     xlim=c(min(software.testing.diff.2$t), +     xlim=c(min(software.debugging.diff.2$t), 
-            max(software.testing.diff.2$t)), +            max(software.debugging.diff.2$t)), 
-     ylim=c(min(software.testing.diff.2$dNt), +     ylim=c(min(software.debugging.diff.2$dNt), 
-            max(software.testing.diff.2$dNt)),+            max(software.debugging.diff.2$dNt)),
      pch=20,col="blue")      pch=20,col="blue")
-points(software.testing.diff.2$t[software.testing.diff.2$dNt<0], +points(software.debugging.diff.2$t[software.debugging.diff.2$dNt<0], 
-       software.testing.diff.2$dNt[software.testing.diff.2$dNt<0],+       software.debugging.diff.2$dNt[software.debugging.diff.2$dNt<0],
        pch=20,col="red")        pch=20,col="red")
 </code> </code>
行 477: 行 485:
  
 <code> <code>
-plot(software.testing.diff$t,software.testing.diff$Ct.diff,+plot(software.debugging.diff$t,software.debugging.diff$Ct.diff,
      pch=20)      pch=20)
 </code> </code>
行 484: 行 492:
  
 <code> <code>
-plot(software.testing.diff$t,software.testing.diff$dCt,+plot(software.debugging.diff$t,software.debugging.diff$dCt,
      pch=20)      pch=20)
 </code> </code>
行 491: 行 499:
  
 <code> <code>
-software.testing.gam <- gam(dCt~-1+bs(t), data=software.testing.diff) +software.debugging.gam <- gam(dCt~-1+bs(t), data=software.debugging.diff) 
-plot(software.testing.gam,+plot(software.debugging.gam,
      xlab="Staff Days", ylab="Correction Lines Per Staff Day")      xlab="Staff Days", ylab="Correction Lines Per Staff Day")
-points(software.testing.diff$t, software.testing.diff$dCt, pch=20)+points(software.debugging.diff$t, software.debugging.diff$dCt, pch=20)
 </code> </code>
  
行 505: 行 513:
  
 <code> <code>
-software.testing.diff <- rbind(c(0,0,0,0,0,0),software.testing.diff)+software.debugging.diff <- rbind(c(0,0,0,0,0,0),software.debugging.diff)
 </code> </code>
  
行 516: 行 524:
   return(diff.1*diff.2)   return(diff.1*diff.2)
 } }
-#lambda.t(0.001704,19.1745,0.003045,10,software.testing.diff)+#lambda.t(0.001704,19.1745,0.003045,10,software.debugging.diff)
 </code> </code>
  
行 525: 行 533:
   alpha <- x[2]   alpha <- x[2]
   beta  <- x[3]   beta  <- x[3]
-  J <- dim(software.testing.diff)[1]+  J <- dim(software.debugging.diff)[1]
   log.lik.temp <- 0   log.lik.temp <- 0
   for( j in c(2:J) ) {   for( j in c(2:J) ) {
-    lambda.j <- lambda.t(theta,alpha,beta,j, software.testing.diff)+    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 - lambda.j
-    log.lik.temp <- log.lik.temp + software.testing.diff$Nt.diff[j]*log(lambda.j)+    log.lik.temp <- log.lik.temp + software.debugging.diff$Nt.diff[j]*log(lambda.j)
   }   }
   return(-log.lik.temp)   return(-log.lik.temp)
行 550: 行 558:
 以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。 以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。
  
 +=== モデルを当てはめる ===
 +
 +  * このデータの場合、データ観測系列は1本なので、ロバストな推定方法は使わない。
 +  * 先週の説明の通り、飽和しかけているので、曲線を強度関数に用いる (累積データに「飽和する曲線」を当てはめる場合、変曲点後まで観測できているか否かで、当てはまりがとても異なるので注意)
 +
 +データの先頭に1行、0を入れておく。
 +
 +<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)掲載の推定値を参考に設定した。
 +以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。
 +
 +== 補足 ==
 +optimという関数は、目的関数と初期値を与えると最適化してくれる関数、で、幾つかの有名な最適化アルゴリズムを利用できて便利なのですが、探索範囲を制限することができない。
 +
 +今回の演習で用いているのはNelder-Meed法という、Simplex法に似たアルゴリズムで、初期値から出発して、探索していきます。その際に、大きな歩幅で探索すると、warningが表示されるように計算きない点に到達する。そのため、今回のような「計算可能な探索範囲」が狭い最適化問題では、目的関数を評価できない旨のwarningが沢山、ログに残る。ただし、初期値の近傍の局所最適解に到達している場合には、その解が表示されるので、warningがある旨のメッセージではなく、最終的な結果が出力されているか否かで、最適化が行えたかどうか、判断して欲しい。
  
 +(質問を受けたので補足)
 === 強度関数を変えてみる === === 強度関数を変えてみる ===
  
行 565: 行 631:
   return(diff.1*diff.2)   return(diff.1*diff.2)
 } }
-lambda.t(0.001704,19.1745,0.003045,0.003045,10,software.testing.diff)+lambda.t(0.001704,19.1745,0.003045,0.003045,10,software.debugging.diff)
 </code> </code>
  
行 574: 行 640:
   beta  <- x[3]   beta  <- x[3]
   beta.2 <- x[4]   beta.2 <- x[4]
-  J <- dim(software.testing.diff)[1]+  J <- dim(software.debugging.diff)[1]
   log.lik.temp <- 0   log.lik.temp <- 0
   for( j in c(2:J) ) {   for( j in c(2:J) ) {
-    lambda.j <- lambda.t(theta,alpha,beta,beta.2,j, software.testing.diff)+    lambda.j <- lambda.t(theta,alpha,beta,beta.2,j, software.debugging.diff)
     log.lik.temp <- log.lik.temp - lambda.j     log.lik.temp <- log.lik.temp - lambda.j
-    log.lik.temp <- log.lik.temp + software.testing.diff$Nt.diff[j]*log(lambda.j)+    log.lik.temp <- log.lik.temp + software.debugging.diff$Nt.diff[j]*log(lambda.j)
   }   }
   return(-log.lik.temp)   return(-log.lik.temp)
行 612: 行 678:
 時点jと時点j-1の間の強度は 時点jと時点j-1の間の強度は
 <code> <code>
-lambda.t(fitted$par[1], fitted$par[2], fitted$par[3], fitted$par[4],j,software.testing.diff)+lambda.t(fitted$par[1], fitted$par[2], fitted$par[3], fitted$par[4],j,software.debugging.diff)
 </code> </code>
 で求められる。テストデータに付与するには、 で求められる。テストデータに付与するには、
 <code> <code>
 fitted.int <- function(theta, alpha, beta, beta.2, data) { fitted.int <- function(theta, alpha, beta, beta.2, data) {
-  lambda <- software.testing.diff$t+  lambda <- software.debugging.diff$t
   lambda <- 0   lambda <- 0
-  J <- dim(software.testing.diff)[1]+  J <- dim(software.debugging.diff)[1]
   for( j in c(2:J) ) {   for( j in c(2:J) ) {
     lambda[j] <- lambda.t(theta, alpha, beta, beta.2, j, data)     lambda[j] <- lambda.t(theta, alpha, beta, beta.2, j, data)
行 632: 行 698:
 <code> <code>
 fitted.diff <- fitted.int(fitted$par[1], fitted$par[2], fitted$par[3], fitted$par[4], fitted.diff <- fitted.int(fitted$par[1], fitted$par[2], fitted$par[3], fitted$par[4],
-                   software.testing.diff)+                   software.debugging.diff)
 </code> </code>
  
行 672: 行 738:
 今のところ、これが一番良いモデル。 今のところ、これが一番良いモデル。
 他にないか、頑張ってみて。 他にないか、頑張ってみて。
 +
 +=== テスト用コード ===
 +<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>
 +