差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
dm:2011 [2011/06/15 13:16] – [治療群毎の強度関数の推定] wataludm:2011 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1
行 141: 行 141:
  
 <code> <code>
-NPfit <- coxph(Surv(start,stop,event)~1, data=bladder.data, subset=(rx==1),method="breslow") +par(mfrow=c(2,1)) 
-KM <- survfit(NPfit, conf.int=.95, type="aalen") +NPfit.1 <- coxph(Surv(start,stop,event)~1, data=bladder.data, subset=(rx==1),method="breslow") 
-plot(KM) +KM.1 <- survfit(NPfit.1, conf.int=.95, type="aalen") 
-NPfit <- coxph(Surv(start,stop,event)~rx, data=bladder.data, subset=(rx==0),method="breslow") +plot(KM.1) 
-KM <- survfit(NPfit, conf.int=.95, type="aalen") +NPfit.2 <- coxph(Surv(start,stop,event)~rx, data=bladder.data, subset=(rx==0),method="breslow") 
-lines(KM, lty=3)+KM.2 <- survfit(NPfit.2, conf.int=.95, type="aalen") 
 +lines(KM.2, lty=3) 
 +plot(KM.2) 
 +lines(KM.1, lty=3) 
 +</code>
  
-NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv))) +95%信頼区間を重ね書きできないので、Treatment群の生存関数と信頼区間を描いてPlacebo群の生存関数を上書きしたものと、 
-plot(NA.MF, conf.int=TRUE, type="l", lty=1) +Placebo群の生存関数と信頼区間を描いてTreatment群の生存関数を上書きしたもの、の2枚を用意した。
-</code>+
  
-次に、上記と同等だが、それぞれの強度関数を推定してみる。+次に、上記と同等だが、それぞれの累積強度関数を推定してみる。
 ここではNelson-Aalen推定量という推定方法を用いているが、詳細は省略する。 ここではNelson-Aalen推定量という推定方法を用いているが、詳細は省略する。
  
 <code> <code>
 +par(mfrow=c(1,1))
 NPfit <- coxph(Surv(start,stop,event)~1, data=bladder.data, subset=(rx==1),method="breslow") NPfit <- coxph(Surv(start,stop,event)~1, data=bladder.data, subset=(rx==1),method="breslow")
 KM <- survfit(NPfit, conf.int=.95, type="aalen") KM <- survfit(NPfit, conf.int=.95, type="aalen")
行 244: 行 248:
 <code> <code>
 NPfit <- coxph(Surv(start, stop, event)~factor(rx)+number+size NPfit <- coxph(Surv(start, stop, event)~factor(rx)+number+size
-               +factor(rx)1:number+factor(rx)1:size+number:size+factor(rx)1:number:size+cluster(id), +               +factor(rx):number+factor(rx):size+number:size+factor(rx):number:size+cluster(id), 
                data=bladder.data, method="breslow")                data=bladder.data, method="breslow")
 summary(NPfit) summary(NPfit)
行 276: 行 280:
  
 このデータの再来院に関して、以上の分析から、何が言えるか、まとめよ。 このデータの再来院に関して、以上の分析から、何が言えるか、まとめよ。
 +
 +==== 演習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時間、の二種類の与え方は例示している。
 +
 +これを更に拡張するなど、できたら試みて欲しいとは思っている。
 +
 +
 +=== 準備 ===
 +
 +データの読み込み
 +<code>
 +software.debugging <- read.table("http://stat.inf.uec.ac.jp/library/dm.2011/software.txt",
 +                           sep=",", header=T)
 +</code>
 +
 +このデータは3つのフィールドを持つ。
 +|t|Nt|Ct|
 +|テスト時間(人時間)|累積不具合発見数|累積改修行数|
 +
 +今回は二つのパッケージを用いるので、インストールしておく。
 +<code>
 +install.packages(pkgs=c("gam"), 
 +                 repos='http://cran.md.tsukuba.ac.jp/',
 +                 dependencies=TRUE)
 +library(gam)
 +</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>
 +