==== 演習課題2 2012.07.26 ====
=== 日程など ===
|出題|2012.07.26|
|〆切|2012.08.10|

  * 提出先は watalu.yamamoto_at_gmail まで電子メールに添付して送付すること。
  * レポートの作成はWordでもGoogle DocsでもLibreOfficeでもLaTeXでもWikiでも、何を用いても良いが、送付するのはPDFファイルかWordファイルにして欲しい。

=== 今回の解析の流れ ===

  - 解析対象は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時間、の二種類の与え方は例示している。

これを更に拡張するなど、できたら試みて欲しいとは思っている。


=== 準備 ===

データファイルは http://stat.inf.uec.ac.jp/library/dm.2011/software.txt に置いてある。そこからダウンロードしても良いし、下記のコードを実行してRに直接読み込ませることもできる。
<code>
software.debugging <- read.table("http://stat.inf.uec.ac.jp/library/dm.2011/software.txt",
                           sep=",", header=T)
</code>

先週、紹介した通り、このデータは3つのフィールドから成る。

|t|Nt|Ct|
|テスト時間(人時間)|累積不具合発見数|累積改修行数|
|0.0|0|0|
|4.8|0|16012|
|6.0|0|16012|
|以下略|||

今回の演習では、尤度関数を自分で書き、それをRのoptim()という最適化の関数を用いて最大化し、パラメータを推定する。
その途中で、平滑化のグラフを作成するのに、gamパッケージを用いるので、最初にインストールしておく。

<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")
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")
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>

=== まずは不具合発見数の成長を眺める ===

累積不具合発見数のグラフは次の一行で描ける。

<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>


=== モデルを当てはめる ===

今回は、用いるモデルがRには用意されていないので、

  * モデルの強度関数と、強度関数を含む対数尤度関数を、Rの関数として記述する。
  * 対数尤度関数を、optim()を用いて最大化する

という手順で、パラメータの推定値を得る。

  * このデータの場合、データ観測系列は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がある旨のメッセージではなく、最終的な結果が出力されているか否かで、最適化が行えたかどうか、判断して欲しい。


=== 強度関数を変えてみる ===

ここでは、強度関数の変え方の一例を示す。

  * 強度関数を変えるにはlambda.t()だけでなく、neg.log.lik()も変える必要が出てくることもある。今回は、パラメータを追加したので、neg.log.lik()の中で x[4] の扱いを追加し、4変数関数としなければならない。
  * さらに閾値を変えてみるとして、betaとbeta.2の切り替え点を300から500に変えるなど、してみている。

過去の修正の影響がもっと軽い可能性はないかと、betaを2種類にしてみる。
まずは300日以内の影響をbetaとし、300日以降の影響を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>

強度関数のパラメータがひとつ増えたので、対数尤度関数の引数xの次元もひとつ増える。
そのため、一行追加する。

<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>

直近の改修よりも、以前の改修の方が、不具合の発生に大きく影響を与える、という母数の推定値が得られる。
これは、初期に多くの改修を行ったことと、無関係ではないかもしれない。

以下はbetaとbeta.2の切り替え時点を300から500に変更してみる場合の強度関数である。

<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は、モデルの当てはまりの良さを評価する指標のひとつである。

  * AICは予測精度の近似として導出されている。
  * AICを基準にモデルを選ぶ、とはこの値が最小となるように、モデルを選ぶことである。
  * AICを最小化するモデルが、必ずしも正解とは限らない。

今回は、次のように計算される量を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>

今のところ、これが一番良いモデル。
他にないか、頑張ってみて。

=== レポートで検討して欲しいこと ===

Q1. 作成した幾つかのグラフに基づいてデータについて考察せよ。

Q2. 幾つかのモデルを当てはめ、AICを計算し、残差も眺めながら検討せよ。

Q3. 最終的にひとつ、モデルを提案せよ。