差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
dm:2012 [2012/06/27 09:52] – [演習課題1 2012.06.14] wataludm:2012 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1
行 104: 行 104:
 LaTeXで書いたテキストをdvipdfmxでPDFに変換すると、なぜか「~」の記号のフォントが、置き換えられてしまいます。見た目はほとんど変わりませんが、半角の「~」では無くなっています。改行コードが余計につくこともあります。 LaTeXで書いたテキストをdvipdfmxでPDFに変換すると、なぜか「~」の記号のフォントが、置き換えられてしまいます。見た目はほとんど変わりませんが、半角の「~」では無くなっています。改行コードが余計につくこともあります。
 テキストのコードをそのままRに貼り付けると、プログラムとして認識できない文字が入力されたことになりますので、「~」は半角文字で打ち直すようにして下さい。 テキストのコードをそのままRに貼り付けると、プログラムとして認識できない文字が入力されたことになりますので、「~」は半角文字で打ち直すようにして下さい。
-問題なく動作するコードを下に記します。+あるいはPDFファイルではなく、このページからコピー＆ペーストして下さい。 
 + 
 +課題1について、問題なく動作するコードを下に記します。
  
 <code> <code>
行 148: 行 150:
 |提出先|mailto:watalu.yamamoto_at_gmail.com (_at_を＠に交換してください)| |提出先|mailto:watalu.yamamoto_at_gmail.com (_at_を＠に交換してください)|
 |提出方法|PDFファイルに変換してメールに添付して送付してください| |提出方法|PDFファイルに変換してメールに添付して送付してください|
 +|提出確認|たまにまとめて返信します|
 +
  
  
行 332: 行 336:
 このデータの再来院に関して、以上の分析から、何が言えるか、まとめよ。 このデータの再来院に関して、以上の分析から、何が言えるか、まとめよ。
  
-==== 演習課題2 ====+==== 演習課題2 2012.07.26 ==== 
 +=== 日程 === 
 +|出題|2012.07.26| 
 +|〆切|2012.08.10| 
 + 
 === 今回の解析の流れ === === 今回の解析の流れ ===
  
行 400: 行 409:
 === 準備 === === 準備 ===
  
-データの読み込み+データファイルは http://stat.inf.uec.ac.jp/library/dm.2011/software.txt に置いてある。そこからダウンロードしても良いし、下記のコードを実行してRに直接読み込ませることもできる。
 <code> <code>
 software.debugging <- 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",
行 406: 行 415:
 </code> </code>
  
-このデータは3つのフィールドを持つ。+先週、紹介した通り、このデータは3つのフィールドから成る。 
 |t|Nt|Ct| |t|Nt|Ct|
 |テスト時間(人時間)|累積不具合発見数|累積改修行数| |テスト時間(人時間)|累積不具合発見数|累積改修行数|
 +|0.0|0|0|
 +|4.8|0|16012|
 +|6.0|0|16012|
 +|以下略|||
 +
 +今回の演習では、尤度関数を自分で書き、それをRのoptim()という最適化の関数を用いて最大化し、パラメータを推定する。
 +その途中で、平滑化のグラフを作成するのに、gamパッケージを用いるので、最初にインストールしておく。
  
-今回は二つのパッケージを用いるので、インストールしておく。 
 <code> <code>
 install.packages(pkgs=c("gam"),  install.packages(pkgs=c("gam"), 
行 418: 行 434:
 </code> </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")
 +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
 +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>
 === まずは不具合発見数の成長を眺める === === まずは不具合発見数の成長を眺める ===
  
行 550: 行 615:
 points(software.debugging.diff$t, software.debugging.diff$dCt, pch=20) points(software.debugging.diff$t, software.debugging.diff$dCt, pch=20)
 </code> </code>
 +
  
 === モデルを当てはめる === === モデルを当てはめる ===
  
-  * このデータの場合、データ観測系列は1本なので、ロバストな推定方法は使わない。 +今回は、用いるモデルがRには用意されていないので、
-  * 先週の説明の通り、飽和しかけているので、曲線を強度関数に用いる (累積データに「飽和する曲線」を当てはめる場合、変曲点後まで観測できているか否かで、当てはまりがとても異なるので注意)+
  
-データの先頭に1行、0を入れておく。+  * モデルの強度関数と、強度関数を含む対数尤度関数を、Rの関数として記述する。 
 +  * 対数尤度関数を、optim()を用いて最大化する
  
-<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本なので、ロバストな推定方法は使わない。
行 658: 行 678:
  
 == 補足 == == 補足 ==
-optimという関数は、目的関数と初期値を与えると最適化してくれる関数、で、幾つかの有名な最適化アルゴリズムを利用できて便利なのですが、探索範囲を制限することができない。+optimという関数は、目的関数と初期値を与えると最適化してくれる関数、で、幾つかの有名な最適化アルゴリズムを利用できて便利だが、探索範囲を制限することができない。 
 + 
 +今回の演習で用いているのはNelder-Meed法という、Simplex法に似たアルゴリズムで、初期値から出発して、探索していく。その際に、大きな歩幅で探索すると、warningが表示されるように計算きない点に到達する。そのため、今回のような「計算可能な探索範囲」が狭い最適化問題では、目的関数を評価できない旨のwarningが沢山、ログに残る。ただし、初期値の近傍の局所最適解に到達している場合には、その解が表示されるので、warningがある旨のメッセージではなく、最終的な結果が出力されているか否かで、最適化が行えたかどうか、判断して欲しい。
  
-今回の演習で用いているのはNelder-Meed法という、Simplex法に似たアルゴリズムで、初期値から出発して、探索していきます。その際に、大きな歩幅で探索すると、warningが表示されるように計算きない点に到達する。そのため、今回のような「計算可能な探索範囲」が狭い最適化問題では、目的関数を評価できない旨のwarningが沢山、ログに残る。ただし、初期値の近傍の局所最適解に到達している場合には、その解が表示されるので、warningがある旨のメッセージではなく、最終的な結果が出力されているか否かで、最適化が行えたかどうか、判断して欲しい。 
  
-(質問を受けたので補足) 
 === 強度関数を変えてみる === === 強度関数を変えてみる ===
 +
 +ここでは、強度関数の変え方の一例を示す。
 +
 +  * 強度関数を変えるにはlambda.t()だけでなく、neg.log.lik()も変える必要が出てくることもある。今回は、パラメータを追加したので、neg.log.lik()の中で x[4] の扱いを追加し、4変数関数としなければならない。
 +  * さらに閾値を変えてみるとして、betaとbeta.2の切り替え点を300から500に変えるなど、してみている。
  
 過去の修正の影響がもっと軽い可能性はないかと、betaを2種類にしてみる。 過去の修正の影響がもっと軽い可能性はないかと、betaを2種類にしてみる。
 +まずは300日以内の影響をbetaとし、300日以降の影響をbeta.2とするように、強度関数を変更する。
  
 <code> <code>
行 679: 行 705:
 lambda.t(0.001704,19.1745,0.003045,0.003045,10,software.debugging.diff) lambda.t(0.001704,19.1745,0.003045,0.003045,10,software.debugging.diff)
 </code> </code>
 +
 +強度関数のパラメータがひとつ増えたので、対数尤度関数の引数xの次元もひとつ増える。
 +そのため、一行追加する。
  
 <code> <code>
行 699: 行 728:
 直近の改修よりも、以前の改修の方が、不具合の発生に大きく影響を与える、という母数の推定値が得られる。 直近の改修よりも、以前の改修の方が、不具合の発生に大きく影響を与える、という母数の推定値が得られる。
 これは、初期に多くの改修を行ったことと、無関係ではないかもしれない。 これは、初期に多くの改修を行ったことと、無関係ではないかもしれない。
 +
 +以下はbetaとbeta.2の切り替え時点を300から500に変更してみる場合の強度関数である。
  
 <code> <code>
行 711: 行 742:
 } }
 </code> </code>
 +
 +このような変更手段だけでも、かなりのバリエーションにはなるので、いろいろ試してみて欲しい。
  
 === モデルを選ぶ === === モデルを選ぶ ===
  
-AICは、モデルの予測精度を評価する指標として導出されており、次のように計算できる。+AICは、モデルの当てはまりの良さを評価する指標のひとつである。 
 + 
 +  * AICは予測精度の近似として導出されている。 
 +  * AICを基準にモデルを選ぶ、とはこの値が最小となるように、モデルを選ぶことである。 
 +  * AICを最小化するモデルが、必ずしも正解とは限らない。 
 + 
 +今回は、次のように計算される量をAICとして用いる。 
 <code> <code>
 fitted <- optim(c(0.001,19,0.005,0.005),neg.log.lik) fitted <- optim(c(0.001,19,0.005,0.005),neg.log.lik)
 print(fitted$value*2+2*length(fitted$par)) print(fitted$value*2+2*length(fitted$par))
 </code> </code>
 +
 +これをなるべく小さくするようなモデルを用いるのが望ましいが、対象についての知識・情報がある場合には、それに基づいて、最小ではないモデルを選択することも少なくない。
  
 === 残差をプロットする === === 残差をプロットする ===
 +
 +残差とは、推定したモデルと推定に用いたデータとの差を表す量である。モデルで説明されなかった、残りの部分、とも言える。
 +
 +  * 残差のばらつきが満足できる小ささであること (予測に使えそうな精度であること)
 +  * 残差が何らかの傾向を持たないこと (他に考慮に入れるべき共変量がないこと)
 +
 +は、モデルを推定した後で確認しなければならない。
  
 時点jと時点j-1の間の強度は 時点jと時点j-1の間の強度は
行 785: 行 834:
 他にないか、頑張ってみて。 他にないか、頑張ってみて。
  
-=== テスト用コード === 
-<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>