差分
このページの2つのバージョン間の差分を表示します。
両方とも前のリビジョン前のリビジョン次のリビジョン | 前のリビジョン最新のリビジョン両方とも次のリビジョン | ||
dm:2011 [2011/07/21 01:50] – watalu | dm:2011 [2011/08/02 10:28] – [補足] watalu | ||
---|---|---|---|
行 282: | 行 282: | ||
==== 演習2 ==== | ==== 演習2 ==== | ||
+ | === 日程 === | ||
+ | |出題|2011.07.21| | ||
+ | |〆切|2011.08.04| | ||
+ | |||
+ | |||
+ | === 今回の解析の流れ === | ||
+ | |||
+ | - 解析対象はSoftware Debugging Data (Dalal and McIntosh, 1994; Cook and Lawless, 2007, Sec. 1.2.2 and Appendix C) | ||
+ | - まずはデータを、テスト時間 vs 累積発見不具合数、テスト時間 vs 累積発見不具合数の一階差分、テスト時間 vs 累積発見不具合数の二階差分、テスト時間 vs 累積改修行数、テスト時間 vs 累積改修行数の一階差分、など図示しながら、吟味していき、観測期間が変曲点を通り過ぎている(観測期間の終わりの方で累積のグラフの傾きが減少傾向にある)ことを確認する | ||
+ | - 飽和する成長曲線モデルを当てはめる、今回は観測系列が1系列のみなので、非定常ポアソン過程を用いる (ロバスト分散などは考えないし、検定なども行わず、強度関数のモデルを工夫するだけ、という宣言) | ||
+ | - 当てはめた累積強度関数とデータを重ねて描く、マルチンゲール残差、あるいはポアソン残差などの残差プロットを描く、などしつつAICも計算しておく (これらの残差プロットを、Cook and Lawlessでは「特にモデルに瑕疵があるとは思えない」としていることを、記しておく) | ||
+ | - 上のモデルをベースに、モデルを少しずつ複雑にしてみて、どのようなモデルが適当か、調べていく | ||
+ | |||
+ | 最初に当てはめる強度関数は | ||
+ | |||
+ | < | ||
+ | |||
+ | である。この初項は飽和する曲線で、第二項は直近の改修が大きな影響を、古い改修ほど小さな影響を与えるように指数平滑化したような関数であり、改修行数を計数過程の強度関数に反映させるよう提案されたモデルである。 | ||
+ | |||
+ | これは、期間 < | ||
+ | |||
+ | < | ||
+ | |||
+ | と一定になる、と書き直すことができる。すると、この期間内の期待発見件数は | ||
+ | |||
+ | < | ||
+ | |||
+ | となる。これを < | ||
+ | |||
+ | < | ||
+ | L\left(\theta, | ||
+ | </ | ||
+ | |||
+ | であることから、 | ||
+ | |||
+ | < | ||
+ | \log L = \sum_{j=1}^k \left(-\mu_j+N_j \log \mu_j\right) | ||
+ | </ | ||
+ | |||
+ | を未知母数 <jsm> \theta, \alpha, \beta </ | ||
+ | |||
+ | < | ||
+ | AIC\left(\hat{\theta}, | ||
+ | </ | ||
+ | |||
+ | ただし< | ||
+ | |||
+ | < | ||
+ | N_j-\hat{\mu}_j=N_j-\left(e^{-\hat{\theta} t_{j-1}}-e^{-\hat{\theta} t_j}\right)\left(\hat{\alpha} + \hat{\beta} \sum_{l=0}^{j-1} c_l e^{\hat{\theta} t_l}\right) | ||
+ | </ | ||
+ | |||
+ | もしくは、 | ||
+ | |||
+ | < | ||
+ | \left(N_j-\hat{\mu}_j\right)/ | ||
+ | </ | ||
+ | |||
+ | 今回の解析では、Cook and Lawless (2007)では、上の強度関数をこのデータに対して当てはめているが、これを少しだけ精緻化することを試みてもらう。 | ||
+ | |||
+ | 下記では、 | ||
+ | |||
+ | < | ||
+ | \lambda\left(t|H\left(t\right)\right)=\theta e^{-\theta t}\left\{\alpha+\beta\sum_{l=1}^{d} c_l e^{\theta t_l}+\beta_2\sum_{l=d+1}^{j-1} c_l e^{\theta t_l}\right\} | ||
+ | </ | ||
+ | |||
+ | と、改修履歴に依存する部分を < | ||
+ | |||
+ | これを更に拡張するなど、できたら試みて欲しいとは思っている。 | ||
+ | |||
=== 準備 === | === 準備 === | ||
行 297: | 行 366: | ||
今回は二つのパッケージを用いるので、インストールしておく。 | 今回は二つのパッケージを用いるので、インストールしておく。 | ||
< | < | ||
+ | install.packages(pkgs=c(" | ||
+ | | ||
+ | | ||
library(gam) | library(gam) | ||
</ | </ | ||
行 305: | 行 377: | ||
< | < | ||
- | plot(software.testing$t, software.testing$Nt, | + | plot(software.debugging$t, software.debugging$Nt, |
| | ||
| | ||
行 314: | 行 386: | ||
< | < | ||
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, |
| | ||
- | points(software.testing$t, software.testing$Nt, | + | points(software.debugging$t, software.debugging$Nt, |
| | ||
</ | </ | ||
行 327: | 行 399: | ||
< | < | ||
- | 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: | + | t.diff=software.debugging$t[2: |
- | Nt.diff=software.testing$Nt[2: | + | Nt.diff=software.debugging$Nt[2: |
- | Ct.diff=software.testing$Ct[2: | + | Ct.diff=software.debugging$Ct[2: |
- | software.testing.diff$dCt <- software.testing.diff$Ct.diff/ | + | software.debugging.diff$dCt <- software.debugging.diff$Ct.diff/ |
- | software.testing.diff$dNt <- software.testing.diff$Nt.diff/ | + | software.debugging.diff$dNt <- software.debugging.diff$Nt.diff/ |
</ | </ | ||
行 339: | 行 411: | ||
< | < | ||
- | 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, |
| | ||
- | points(software.testing.diff$t, | + | points(software.debugging.diff$t, |
| | ||
</ | </ | ||
行 351: | 行 423: | ||
二階差分も確認してみる。 | 二階差分も確認してみる。 | ||
< | < | ||
- | n <- dim(software.testing.diff)[1] | + | n <- dim(software.debugging.diff)[1] |
- | software.testing.diff.2 <- data.frame(t=software.testing.diff$t[2: | + | software.debugging.diff.2 <- data.frame(t=software.debugging.diff$t[2: |
- | t.diff=software.testing.diff$t[2: | + | t.diff=software.debugging.diff$t[2: |
- | Nt.diff.2=software.testing.diff$Nt[2: | + | Nt.diff.2=software.debugging.diff$Nt[2: |
- | Ct.diff.2=software.testing.diff$Ct[2: | + | Ct.diff.2=software.debugging.diff$Ct[2: |
- | software.testing.diff.2$dCt <- software.testing.diff.2$Ct.diff.2/ | + | software.debugging.diff.2$dCt <- software.debugging.diff.2$Ct.diff.2/ |
- | software.testing.diff.2$dNt <- software.testing.diff.2$Nt.diff.2/ | + | software.debugging.diff.2$dNt <- software.debugging.diff.2$Nt.diff.2/ |
</ | </ | ||
行 363: | 行 435: | ||
< | < | ||
- | 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, |
| | ||
- | points(software.testing.diff.2$t, | + | points(software.debugging.diff.2$t, |
| | ||
</ | </ | ||
行 374: | 行 446: | ||
< | < | ||
- | plot(software.testing.diff.2$t, | + | plot(software.debugging.diff.2$t, |
| | ||
</ | </ | ||
行 383: | 行 455: | ||
< | < | ||
- | 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, |
| | ||
- | points(software.testing.diff.2$t, | + | points(software.debugging.diff.2$t, |
| | ||
</ | </ | ||
行 394: | 行 466: | ||
< | < | ||
- | plot(software.testing.diff.2$t[software.testing.diff.2$dNt> | + | plot(software.debugging.diff.2$t[software.debugging.diff.2$dNt> |
- | | + | |
- | | + | |
- | max(software.testing.diff.2$t)), | + | max(software.debugging.diff.2$t)), |
- | | + | |
- | max(software.testing.diff.2$dNt)), | + | max(software.debugging.diff.2$dNt)), |
| | ||
- | points(software.testing.diff.2$t[software.testing.diff.2$dNt< | + | points(software.debugging.diff.2$t[software.debugging.diff.2$dNt< |
- | | + | |
| | ||
</ | </ | ||
行 413: | 行 485: | ||
< | < | ||
- | plot(software.testing.diff$t, | + | plot(software.debugging.diff$t, |
| | ||
</ | </ | ||
行 420: | 行 492: | ||
< | < | ||
- | plot(software.testing.diff$t, | + | plot(software.debugging.diff$t, |
| | ||
</ | </ | ||
行 427: | 行 499: | ||
< | < | ||
- | software.testing.gam <- gam(dCt~-1+bs(t), | + | software.debugging.gam <- gam(dCt~-1+bs(t), |
- | plot(software.testing.gam, | + | plot(software.debugging.gam, |
| | ||
- | points(software.testing.diff$t, software.testing.diff$dCt, pch=20) | + | points(software.debugging.diff$t, software.debugging.diff$dCt, pch=20) |
</ | </ | ||
行 441: | 行 513: | ||
< | < | ||
- | software.testing.diff <- rbind(c(0, | + | software.debugging.diff <- rbind(c(0, |
</ | </ | ||
行 452: | 行 524: | ||
return(diff.1*diff.2) | return(diff.1*diff.2) | ||
} | } | ||
- | lambda.t(0.001704, | + | #lambda.t(0.001704, |
</ | </ | ||
行 461: | 行 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, | + | lambda.j <- lambda.t(theta, |
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) | ||
行 483: | 行 555: | ||
</ | </ | ||
- | のように与えると、それなりの値に収束する。 | + | のように与えると、それなりの値に収束する。optim()の最初の引数は、最適化の初期値である。これは Cook and Lawless (2007)掲載の推定値を参考に設定した。 |
+ | 以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。 | ||
+ | === モデルを当てはめる === | ||
+ | |||
+ | * このデータの場合、データ観測系列は1本なので、ロバストな推定方法は使わない。 | ||
+ | * 先週の説明の通り、飽和しかけているので、曲線を強度関数に用いる (累積データに「飽和する曲線」を当てはめる場合、変曲点後まで観測できているか否かで、当てはまりがとても異なるので注意) | ||
+ | |||
+ | データの先頭に1行、0を入れておく。 | ||
+ | |||
+ | < | ||
+ | software.debugging.diff <- rbind(c(0, | ||
+ | </ | ||
+ | |||
+ | 強度関数は,次の通り。 | ||
+ | |||
+ | < | ||
+ | lambda.t <- function(theta, | ||
+ | diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j])) | ||
+ | diff.2 <- (alpha+beta*sum(data$Ct.diff[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, | ||
+ | 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) | ||
+ | } | ||
+ | </ | ||
+ | |||
+ | 対数尤度の最大化は、上の関数の最小化と等しい。 | ||
+ | |||
+ | Rには最小化の関数が幾つかある。nlminb(), | ||
+ | |||
+ | これを最小化するのは、結構、困難であるが、例えば初期値を | ||
+ | |||
+ | < | ||
+ | fitted <- optim(c(0.001, | ||
+ | print(fitted) | ||
+ | </ | ||
+ | |||
+ | のように与えると、それなりの値に収束する。optim()の最初の引数は、最適化の初期値である。これは Cook and Lawless (2007)掲載の推定値を参考に設定した。 | ||
+ | 以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。 | ||
+ | |||
+ | == 補足 == | ||
+ | optimという関数は、目的関数と初期値を与えると最適化してくれる関数、で、幾つかの有名な最適化アルゴリズムを利用できて便利なのですが、探索範囲を制限することができない。 | ||
+ | |||
+ | 今回の演習で用いているのはNelder-Meed法という、Simplex法に似たアルゴリズムで、初期値から出発して、探索していきます。その際に、大きな歩幅で探索すると、warningが表示されるように計算きない点に到達する。そのため、今回のような「計算可能な探索範囲」が狭い最適化問題では、目的関数を評価できない旨のwarningが沢山、ログに残る。ただし、初期値の近傍の局所最適解に到達している場合には、その解が表示されるので、warningがある旨のメッセージではなく、最終的な結果が出力されているか否かで、最適化が行えたかどうか、判断して欲しい。 | ||
+ | |||
+ | (質問を受けたので補足) | ||
=== 強度関数を変えてみる === | === 強度関数を変えてみる === | ||
行 499: | 行 631: | ||
return(diff.1*diff.2) | return(diff.1*diff.2) | ||
} | } | ||
- | lambda.t(0.001704, | + | lambda.t(0.001704, |
</ | </ | ||
行 508: | 行 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, | + | lambda.j <- lambda.t(theta, |
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) | ||
行 546: | 行 678: | ||
時点jと時点j-1の間の強度は | 時点jと時点j-1の間の強度は | ||
< | < | ||
- | lambda.t(fitted$par[1], | + | lambda.t(fitted$par[1], |
</ | </ | ||
で求められる。テストデータに付与するには、 | で求められる。テストデータに付与するには、 | ||
< | < | ||
fitted.int <- function(theta, | fitted.int <- function(theta, | ||
- | 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, | lambda[j] <- lambda.t(theta, | ||
行 566: | 行 698: | ||
< | < | ||
fitted.diff <- fitted.int(fitted$par[1], | fitted.diff <- fitted.int(fitted$par[1], | ||
- | | + | |
</ | </ | ||
行 606: | 行 738: | ||
今のところ、これが一番良いモデル。 | 今のところ、これが一番良いモデル。 | ||
他にないか、頑張ってみて。 | 他にないか、頑張ってみて。 | ||
+ | |||
+ | === テスト用コード === | ||
+ | < | ||
+ | software.debugging.diff <- data.frame(t=software.debugging$t[2: | ||
+ | t.diff=software.debugging$t[2: | ||
+ | Nt.diff=software.debugging$Nt[2: | ||
+ | Ct.diff=software.debugging$Ct[2: | ||
+ | software.debugging.diff$dCt <- software.debugging.diff$Ct.diff/ | ||
+ | software.debugging.diff$dNt <- software.debugging.diff$Nt.diff/ | ||
+ | software.debugging.diff <- rbind(c(0, | ||
+ | lambda.t <- function(theta, | ||
+ | diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j])) | ||
+ | diff.2 <- (alpha+beta*sum(data$Ct.diff[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, | ||
+ | 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, | ||
+ | print(fitted) | ||
+ | </ | ||