差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
dm:2011 [2011/07/20 21:39] – wataludm:2011 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1
行 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では「特にモデルに瑕疵があるとは思えない」としていることを、記しておく)
 +  - 上のモデルをベースに、モデルを少しずつ複雑にしてみて、どのようなモデルが適当か、調べていく
 +
 +最初に当てはめる強度関数は
 +
 +<jsmath>\lambda\left(t|H\left(t\right)\right)=\alpha\theta e^{-\theta t}+\beta\theta e^{-\theta t} \int_0^t \frac{dC\left(u\right)}{e^{-\theta u}}</jsmath>
 +
 +である。この初項は飽和する曲線で、第二項は直近の改修が大きな影響を、古い改修ほど小さな影響を与えるように指数平滑化したような関数であり、改修行数を計数過程の強度関数に反映させるよう提案されたモデルである。
 +
 +これは、期間 <jsm>[t_{j-1}, t_{j})</jsm> 内の改修がすべて <jsm> t_{j} </jsm> に追加されたとして、この期間内の強度関数は
 +
 +<jsmath>\lambda\left(t|H\left(t\right)\right)=\theta e^{-\theta t}\left\{\alpha+\beta\sum_{l=1}^{j-1} c_l e^{\theta t_l}\right\}</jsmath>
 +
 +と一定になる、と書き直すことができる。すると、この期間内の期待発見件数は
 +
 +<jsmath>E\left[N_j|H_j\right] = \int_{t_{j-1}}^{t_j} \lambda\left(t|H\left(t\right)\right)dt = \left(e^{-\theta t_{j-1}}-e^{-\theta t_j}\right)\left(\alpha + \beta \sum_{l=0}^{j-1} c_l e^{\theta t_l}\right)</jsmath>
 +
 +となる。これを <jsm>\mu_j</jsm> とおくと、ポアソン過程の尤度関数が
 +
 +<jsmath>
 +L\left(\theta, \alpha, \beta\right)=\prod_{j=1}^k e^{-\mu_j}{\mu_j}^N_j
 +</jsmath>
 +
 +であることから、
 +
 +<jsmath>
 +\log L = \sum_{j=1}^k \left(-\mu_j+N_j \log \mu_j\right)
 +</jsmath>
 +
 +を未知母数 <jsm> \theta, \alpha, \beta </jsm> について最大化すれば、良いことがわかる。その他、モデル選択に必要なAICは
 +
 +<jsmath>
 +AIC\left(\hat{\theta}, \hat{\alpha}, \hat{\beta}\right)=-2\log L\left(\hat{\theta}, \hat{\alpha}, \hat{\beta}\right)+2k
 +</jsmath>
 +
 +ただし<jsm>k</jsm>は母数の数である。残差は
 +
 +<jsmath>
 +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)
 +</jsmath>
 +
 +もしくは、
 +
 +<jsmath>
 +\left(N_j-\hat{\mu}_j\right)/\sqrt{\hat{\mu}_j}
 +</jsmath>
 +
 +今回の解析では、Cook and Lawless (2007)では、上の強度関数をこのデータに対して当てはめているが、これを少しだけ精緻化することを試みてもらう。
 +
 +下記では、
 +
 +<jsmath>
 +\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\}
 +</jsmath>
 +
 +と、改修履歴に依存する部分を <jsm>l=1,\ldots,d</jsm> と <jsm>l=d+1,\ldots,j-1</jsm> のように前半と後半に分け、それぞれに母数 <jsm>\beta, \beta_2</jsm> と付与してみる例を与えてある。<jsm>d</jsm>の与え方は、現時点から300時間遡った時点まで、と最初から300時間、の二種類の与え方は例示している。
 +
 +これを更に拡張するなど、できたら試みて欲しいとは思っている。
 +
  
 === 準備 === === 準備 ===
行 287: 行 356:
 データの読み込み データの読み込み
 <code> <code>
-testing.data <- read.table("http://stat.inf.uec.ac.jp/library/dm.2011/software.txt",+software.debugging <- read.table("http://stat.inf.uec.ac.jp/library/dm.2011/software.txt",
                            sep=",", header=T)                            sep=",", header=T)
 </code> </code>
  
-=== 生データ ===+このデータは3つのフィールドを持つ。 
 +|t|Nt|Ct| 
 +|テスト時間(人時間)|累積不具合発見数|累積改修行数|
  
 +今回は二つのパッケージを用いるので、インストールしておく。
 <code> <code>
-t Nt Ct +install.packages(pkgs=c("gam"),  
-0.0 0 0 +                 repos='http://cran.md.tsukuba.ac.jp/', 
-4.8 0 16012 +                 dependencies=TRUE) 
-6.0 0 16012 +library(gam)
-14.3 7 32027 +
-22.8 7 48042 +
-32.1 7 58854 +
-41.4 7 69669 +
-51.2 11 80483 +
-60.6 12 91295 +
-70.0 13 102110 +
-79.9 15 112925 +
-91.3 20 120367 +
-97.0 21 127812 +
-107.7 22 135257 +
-119.1 28 142702 +
-127.6 40 150147 +
-135.1 44 152806 +
-142.8 46 155464 +
-148.9 48 158123 +
-156.6 52 167081 +
-163.9 52 167704 +
-169.7 59 174626 +
-170.1 59 174626 +
-170.6 59 174626 +
-174.7 63 181548 +
-179.6 68 188473 +
-185.5 71 194626 +
-194.0 88 200782 +
-200.3 93 206937 +
-207.2 97 213093 +
-211.9 98 219248 +
-217.0 105 221355 +
-223.5 113 223462 +
-227.0 113 225568 +
-234.1 122 227675 +
-241.6 129 229784 +
-250.7 141 233557 +
-259.8 155 237330 +
-268.3 166 241103 +
-277.2 178 244879 +
-285.5 186 247946 +
-294.2 190 251016 +
-295.7 190 251016 +
-298.0 190 254086 +
-305.2 195 257155 +
-312.3 201 260225 +
-318.2 209 260705 +
-328.9 224 261188 +
-334.8 231 261669 +
-342.7 243 262889 +
-350.5 252 263629 +
-356.3 259 264367 +
-360.6 271 265107 +
-365.7 277 265845 +
-374.9 282 266585 +
-386.5 290 267325 +
-396.5 300 268607 +
-408.0 310 269891 +
-417.3 312 271175 +
-424.9 321 272457 +
-434.2 326 273741 +
-442.7 339 275025 +
-451.4 346 276556 +
-456.1 347 278087 +
-460.8 351 279618 +
-466.0 356 281149 +
-472.3 359 283592 +
-476.4 362 286036 +
-480.9 367 288480 +
-486.8 374 290923 +
-495.8 376 293367 +
-505.7 380 295811 +
-516.0 392 298254 +
-526.2 399 300698 +
-527.3 401 300698 +
-535.8 405 303142 +
-546.3 415 304063 +
-556.1 425 305009 +
-568.1 440 305956 +
-577.2 457 306902 +
-578.3 457 306902 +
-587.2 467 307849 +
-595.5 473 308795 +
-605.6 480 309742 +
-613.9 491 310688 +
-621.6 496 311635 +
-623.4 496 311635 +
-636.3 502 311750 +
-649.7 517 311866 +
-663.9 527 312467 +
-675.1 540 313069 +
-677.4 543 313069 +
-677.9 544 313069 +
-688.4 553 313671 +
-698.1 561 314273 +
-710.5 573 314783 +
-720.9 581 315294 +
-731.6 584 315805 +
-732.7 585 315805 +
-733.6 585 315805 +
-746.7 586 316316 +
-761.0 598 316827 +
-776.5 612 318476 +
-793.5 621 320125 +
-807.2 636 321774 +
-811.8 639 321774 +
-812.5 639 321774 +
-829.0 648 323423 +
-844.4 658 325072 +
-860.5 666 326179 +
-876.7 674 327286 +
-892.0 679 328393 +
-895.5 686 328393 +
-910.8 690 329500 +
-925.1 701 330608 +
-938.3 710 330435 +
-952.0 720 330263 +
-965.0 729 330091 +
-967.7 729 330091 +
-968.6 731 330091 +
-981.3 740 329919 +
-1013.9 759 330036 +
-1030.1 776 330326 +
-1044.0 781 330616 +
-1047.0 782 330616 +
-1059.7 783 330906 +
-1072.6 787 331196 +
-1085.7 793 331486 +
-1098.4 796 331577 +
-1112.4 797 331669 +
-1113.5 798 331669 +
-1114.1 798 331669 +
-1128.0 802 331760 +
-1139.1 805 331852 +
-1163.2 823 332167 +
-1174.3 827 332391 +
-1184.6 832 332615 +
-1198.3 834 332839 +
-1210.3 836 333053 +
-1221.1 839 333267 +
-1230.5 842 333481 +
-1231.6 842 333481 +
-1240.9 844 333695 +
-1249.5 845 333909 +
-1262.2 849 335920 +
-1271.3 851 337932 +
-1279.8 854 339943 +
-1287.4 855 341955 +
-1295.1 859 341967 +
-1304.8 860 341979 +
-1305.8 865 342073 +
-1313.3 867 342168 +
-1314.4 867 342168 +
-1320.0 867 342262 +
-1325.3 867 342357 +
-1330.6 870 342357 +
-1334.2 870 342358 +
-1336.7 870 342358+
 </code> </code>
 +
 +=== まずは不具合発見数の成長を眺める ===
 +
 +累積不具合発見数のグラフは次の一行で描ける。
 +
 +<code>
 +plot(software.debugging$t, software.debugging$Nt,
 +     type="l",
 +     xlab="Staff Days", ylab="Cumulative Faults")
 +</code>
 +
 +特にモデルは用いていないが、平滑化した曲線も付与してみる。
 +
 +<code>
 +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")
 +</code>
 +
 +飽和しているように見えるのは、平滑化の特性なので、ここでは無視する。
 +
 +次に、不具合発見数、改修行数とも、時点毎の差分を取ってみる。
 +
 +<code>
 +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
 +</code>
 +
 +時点間の不具合発見数のグラフを描く。
 +
 +<code>
 +software.debugging.gam <- gam(Nt.diff~-1+bs(t),
 +                            data=software.debugging.diff)
 +plot(software.debugging.gam,
 +     xlab="Staff Days", ylab="1st Order Diff of Cumulative Faults")
 +points(software.debugging.diff$t,software.debugging.diff$Nt.diff,
 +     pch=20)
 +</code>
 +
 +一階差分を平滑化してみると、減少しているように見える。
 +
 +二階差分も確認してみる。
 +<code>
 +n <- dim(software.debugging.diff)[1]
 +software.debugging.diff.2 <- data.frame(t=software.debugging.diff$t[2:n], 
 +                                    t.diff=software.debugging.diff$t[2:n]-software.debugging.diff$t[1:(n-1)],
 +                                    Nt.diff.2=software.debugging.diff$Nt[2:n]-software.debugging.diff$Nt[1:(n-1)],
 +                                    Ct.diff.2=software.debugging.diff$Ct[2:n]-software.debugging.diff$Ct[1:(n-1)])
 +software.debugging.diff.2$dCt <- software.debugging.diff.2$Ct.diff.2/software.debugging.diff.2$t
 +software.debugging.diff.2$dNt <- software.debugging.diff.2$Nt.diff.2/software.debugging.diff.2$t
 +</code>
 +
 +グラフを描いてみる。
 +
 +<code>
 +software.debugging.gam <- gam(Nt.diff.2~-1+bs(t),
 +                            data=software.debugging.diff.2)
 +plot(software.debugging.gam,
 +     xlab="Staff Days", ylab="2nd Order Diff of Cumulative Faults")
 +points(software.debugging.diff.2$t,software.debugging.diff.2$Nt.diff.2,
 +     pch=20)
 +</code>
 +
 +平滑化しないと、こんな感じ。
 +
 +<code>
 +plot(software.debugging.diff.2$t,software.debugging.diff.2$Nt.diff.2,
 +     pch=20)
 +</code>
 +
 +フラットかもしれない。
 +
 +時間あたりに直してみる。
 +
 +<code>
 +software.debugging.gam <- gam(dNt~-1+bs(t),
 +                            data=software.debugging.diff.2)
 +plot(software.debugging.gam,
 +     xlab="Staff Days", ylab="2nd Order Diff of Cumulative Faults per Days")
 +points(software.debugging.diff.2$t,software.debugging.diff.2$dNt,
 +     pch=20)
 +</code>
 +
 +二階差分の非負値を青、負値を赤でプロットしてみると、徐々に負に向かっているようにも見える。
 +
 +<code>
 +plot(software.debugging.diff.2$t[software.debugging.diff.2$dNt>=0],
 +     software.debugging.diff.2$dNt[software.debugging.diff.2$dNt>=0],
 +     xlim=c(min(software.debugging.diff.2$t),
 +            max(software.debugging.diff.2$t)),
 +     ylim=c(min(software.debugging.diff.2$dNt),
 +            max(software.debugging.diff.2$dNt)),
 +     pch=20,col="blue")
 +points(software.debugging.diff.2$t[software.debugging.diff.2$dNt<0],
 +       software.debugging.diff.2$dNt[software.debugging.diff.2$dNt<0],
 +       pch=20,col="red")
 +</code>
 +
 +変曲点があるようにも見えるので、今回は飽和すると思って、分析を続ける。
 +
 +=== 次に改修行数の成長を眺める ===
 +
 +時点間の改修行数のグラフも描いてみる。
 +
 +<code>
 +plot(software.debugging.diff$t,software.debugging.diff$Ct.diff,
 +     pch=20)
 +</code>
 +
 +更に、改修行数を改修に要した期間で割ったグラフも描いてみる。
 +
 +<code>
 +plot(software.debugging.diff$t,software.debugging.diff$dCt,
 +     pch=20)
 +</code>
 +
 +いずれにせよ、初期に大規模な改修が行われているが、徐々に改修行数は減ってきているのが見て取れるはず。
 +
 +<code>
 +software.debugging.gam <- gam(dCt~-1+bs(t), data=software.debugging.diff)
 +plot(software.debugging.gam,
 +     xlab="Staff Days", ylab="Correction Lines Per Staff Day")
 +points(software.debugging.diff$t, software.debugging.diff$dCt, pch=20)
 +</code>
 +
 +=== モデルを当てはめる ===
 +
 +  * このデータの場合、データ観測系列は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)掲載の推定値を参考に設定した。
 +以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。
 +
 +=== モデルを当てはめる ===
 +
 +  * このデータの場合、データ観測系列は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がある旨のメッセージではなく、最終的な結果が出力されているか否かで、最適化が行えたかどうか、判断して欲しい。
 +
 +(質問を受けたので補足)
 +=== 強度関数を変えてみる ===
 +
 +過去の修正の影響がもっと軽い可能性はないかと、betaを2種類にしてみる。
 +
 +<code>
 +lambda.t <- function(theta, alpha, beta, beta.2, j, data) {
 +  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j]))
 +  data.t <- data$t[1:(j-1)]
 +  beta.s <- data.t
 +  beta.s[(data$t[j]-data.t)<=300] <- beta
 +  beta.s[(data$t[j]-data.t)>300] <- beta.2
 +  diff.2 <- (alpha+sum(beta.s*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,0.003045,10,software.debugging.diff)
 +</code>
 +
 +<code>
 +neg.log.lik <- function(x) {
 +  theta <- x[1]
 +  alpha <- x[2]
 +  beta  <- x[3]
 +  beta.2 <- x[4]
 +  J <- dim(software.debugging.diff)[1]
 +  log.lik.temp <- 0
 +  for( j in c(2:J) ) {
 +    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 + software.debugging.diff$Nt.diff[j]*log(lambda.j)
 +  }
 +  return(-log.lik.temp)
 +}
 +</code>
 +
 +直近の改修よりも、以前の改修の方が、不具合の発生に大きく影響を与える、という母数の推定値が得られる。
 +これは、初期に多くの改修を行ったことと、無関係ではないかもしれない。
 +
 +<code>
 +lambda.t <- function(theta, alpha, beta, beta.2, j, data) {
 +  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j]))
 +  data.t <- data$t[1:(j-1)]
 +  beta.s <- data.t
 +  beta.s[(data$t[j]-data.t)<=500] <- beta
 +  beta.s[(data$t[j]-data.t)>500] <- beta.2
 +  diff.2 <- (alpha+sum(beta.s*data$Ct.diff[1:(j-1)]*exp(theta*data$t[1:(j-1)])))
 +  return(diff.1*diff.2)
 +}
 +</code>
 +
 +=== モデルを選ぶ ===
 +
 +AICは、モデルの予測精度を評価する指標として導出されており、次のように計算できる。
 +<code>
 +fitted <- optim(c(0.001,19,0.005,0.005),neg.log.lik)
 +print(fitted$value*2+2*length(fitted$par))
 +</code>
 +
 +=== 残差をプロットする ===
 +
 +時点jと時点j-1の間の強度は
 +<code>
 +lambda.t(fitted$par[1], fitted$par[2], fitted$par[3], fitted$par[4],j,software.debugging.diff)
 +</code>
 +で求められる。テストデータに付与するには、
 +<code>
 +fitted.int <- function(theta, alpha, beta, beta.2, data) {
 +  lambda <- software.debugging.diff$t
 +  lambda <- 0
 +  J <- dim(software.debugging.diff)[1]
 +  for( j in c(2:J) ) {
 +    lambda[j] <- lambda.t(theta, alpha, beta, beta.2, j, data)
 +  }
 +  Lambda <- cumsum(lambda)
 +  return(list(time=data$t,lambda=lambda, Lambda=Lambda, Nt=cumsum(data$Nt.diff)))
 +}
 +</code>
 +
 +を用いて、
 +
 +<code>
 +fitted.diff <- fitted.int(fitted$par[1], fitted$par[2], fitted$par[3], fitted$par[4],
 +                   software.debugging.diff)
 +</code>
 +
 +とする。グラフに描くには、
 +
 +<code>
 +plot(fitted.diff$t, fitted.diff$Nt, type="l")
 +lines(fitted.diff$t, fitted.diff$Lambda, type="l", lty=2)
 +</code>
 +
 +とする。残差のプロットは
 +
 +<code>
 +plot(fitted.diff$t, fitted.diff$Nt-fitted.diff$Lambda, type="p", pch=20)
 +</code>
 +
 +標準化残差のプロットは
 +
 +<code>
 +plot(fitted.diff$t, (fitted.diff$Nt-fitted.diff$Lambda)/sqrt(fitted.diff$Lambda), type="p", pch=20)
 +</code>
 +
 +で得られる。
 +
 +=== 更に強度関数を変えてみる ===
 +
 +<code>
 +lambda.t <- function(theta, alpha, beta, beta.2, j, data) {
 +  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j]))
 +  data.t <- data$t[1:(j-1)]
 +  beta.s <- data.t
 +  beta.s[(data$t[j])<=300] <- beta
 +  beta.s[(data$t[j])>300] <- beta.2
 +  diff.2 <- (alpha+sum(beta.s*data$Ct.diff[1:(j-1)]*exp(theta*data$t[1:(j-1)])))
 +  return(diff.1*diff.2)
 +}
 +</code>
 +
 +今のところ、これが一番良いモデル。
 +他にないか、頑張ってみて。
 +
 +=== テスト用コード ===
 +<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>
 +