差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
dm:2012 [2012/06/14 07:50] – wataludm:2012 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1
行 4: 行 4:
 ==== 演習の準備 2012.06.14 ==== ==== 演習の準備 2012.06.14 ====
 [[http://stat.inf.uec.ac.jp/library/dm.2011/counting-processes-20100708.pdf|このファイル]](昨年の課題)に収められたコードがすべて動くような環境を整えて下さい。 [[http://stat.inf.uec.ac.jp/library/dm.2011/counting-processes-20100708.pdf|このファイル]](昨年の課題)に収められたコードがすべて動くような環境を整えて下さい。
 +具体的には
  
-まずはRをダウンロードしてインストールします。[[http://cran.r-project.org/|CRAN]]は、Rの開発組織が運営しているウェブサイトで、ここからのダウンロードは、とりあえず信用しても良いです。+  * R言語 (PerlやCと似た文法でMATLABのような機能を持つオブジェクト指向言語、処理系はインタプリタでコンパイラは無し) 
 +  * R言語のsurvivalパッケージ (最近はバンドルされてるはず)
  
-  * WindowsかMac OX Sなら、まず[[http://cran.r-project.org/|CRAN]]にアクセスして下さい。 +の2つです。 
-    - Windowsなら、"Windows"の文字をクリックし、次に現れるページで"base"をクリックすると、"Download R X.Y.Z for Windows"という文字が見つかるので、それをクリックしてダウンロードし、実行する。 + 
-    - Mac OS Xなら、"Mac OS X"の文字をクリックし、次に現れるページで "R-X.Y.Z.pkg (latest version)"という文字が見つかるので、それをクリックしてダウンロードし、実行する。+=== はじめに === 
 +演習を行うのにコーディング量が最も少なくて済む環境として、R言語を用います。 
 +ただし、使って貰うのにR言語をすべて理解する必要はなく、こちらが提供したコードを少し改変する範囲で演習課題が終わるような工夫をしています。 
 +演習を進める上で、コードの改変では済まないけどやってみたいこと、があったら、教員に相談してください。 
 + 
 +Rは、1976年に発表され、1980年頃に外部リリースが始まり、1988年にオブジェクト指向が取り入れられた、AT&Tのベル研究所で開発されたS言語というプログラミング言語の派生言語(文法や関数はほぼ共通だが、実装方針が異なるので、クローンとも言える)です。 
 +S言語の実装のひとつなら、Rと呼ぶべきですが、世間に倣ってR言語と呼びます。 
 +SもRも、線形代数などの基本的な数値計算には、世界中で古くから用いられて現在ではとても枯れている[[http://www.netlib.org/|Netlib]]に含まれているライブラリを用いています。 
 +例えば[[http://www.netlib.org/blas/|BLAS]]は1979年、[[http://www.netlib.org/lapack|LAPACK]]は1992年まで、遡ることができますが、いずれもFortran 77で書かれています。 
 +既存のライブラリが広く普及しているのに、新たに実装することで、既知のバグを新たに引き起こしたり、メンテナンスの労力を分散させるのを避けるためです。 
 + 
 +R言語の今日の時点での最新版は、2012/03/30にリリースされた2.15.0です。この言語に関して、[[http://developer.r-project.org/|開発者向けのウェブサイト]]に、[[http://developer.r-project.org/devel-guidelines.txt|ガイドライン]]が掲載されています。そこを読むと、 
 +<code> 
 +DO fix simple bugs 
 +DO NOT fix bugs that require extensive modification 
 +DO NOT fix exotic bugs that haven't bugged anyone 
 +DO make small enhancements if they are badly needed 
 +DO NOT let one user decide what is badly needed 
 +DO fix configuration for broken platforms 
 +DO NOT break functioning platforms in the process 
 + 
 +DO NOT fix things that are not broken 
 +DO NOT restructure code 
 +DO NOT add experimental code 
 +DO NOT modify the API (unless absolutely sure it is buggy) 
 +DO NOT change defaults without a *very* good reason 
 + 
 +DO clarify documentation 
 +ONLY add or modify examples if needed for clarification 
 +DO NOT reword messages 
 + 
 +DO NOT modify regression tests (except if they were buggy) 
 +Be very careful in adding regression tests and consider using them 
 +  in R-devel only at first. 
 +</code> 
 +と書かれています。メジャーリビジョンを変更するようなドラスティックな設計変更や大規模な改修は行わないこと、実験的なコードは含めないこと、不具合の修正と機能の追加をしつつ定期的にマイナーリビジョンを更新する、などの開発方針が読み取れます。よって、使用するパッケージが要求しているバージョン以降であれば、多少古くても大丈夫です。これまでのところ、全ユーザに対して、特定のバージョン以降を用いるように、という推奨は出されていません。 
 + 
 +まずはR言語をインストールします。 
 +インストールは通常、[[http://cran.r-project.org/|CRAN]]からダウンロードして行います。 
 +Rの国際的な開発組織が運営しているウェブサイトです。ここからのダウンロードは、とりあえず信用しても良いです。 
 +このサイトは、[[http://essrc.hyogo-u.ac.jp/cran/mirrors.html|世界中にミラーサイトを持っており]]、一番近いサイトにアクセスすることが推奨されています。 
 +以下では、この大学から一番ダウンロード速度が速い、[[http://cran.ism.ac.jp/|統計数理研究所のミラーサイト]]にリンクを貼りました。 
 + 
 +  * WindowsかMac OX Sなら、まず[[http://cran.ism.ac.jp/|CRAN(のミラーサイト)]]にアクセスして下さい。 
 +    - Windowsなら、"Windows"の文字をクリックし、[[http://cran.ism.ac.jp/bin/windows/|次に現れるページ]]で"base"をクリックすると、"Download R X.Y.Z for Windows"という文字が見つかるので、それをクリックしてダウンロードし、実行する。 
 +    - Mac OS Xなら、"Mac OS X"の文字をクリックし、[[http://cran.ism.ac.jp/bin/macosx/|次に現れるページ]]で "R-X.Y.Z.pkg (latest version)"という文字が見つかるので、それをクリックしてダウンロードし、実行する。
   * Linuxなら、次の2行のどちらかが、インストールしてくれます。sudoが必要かも。<code>   * Linuxなら、次の2行のどちらかが、インストールしてくれます。sudoが必要かも。<code>
 yum install R R-* yum install R R-*
行 18: 行 65:
  
   * survivalパッケージをインストールするには、次の一行を実行します。<code>   * survivalパッケージをインストールするには、次の一行を実行します。<code>
-install.packages(pkgs=c("survival"), repos='http://cran.md.tsukuba.ac.jp/')+install.packages(pkgs=c("survival"), repos='http://cran.ism.ac.jp/')
 </code> </code>
  
行 47: 行 94:
 で設定したプロキシを使用できます。 で設定したプロキシを使用できます。
  
 +その他の環境については、[[http://appl.stat.inf.uec.ac.jp/doku.php?id=r:how_to:internet_proxy|r;how_to:internet_proxy]]が参考になるかもしれません。
  
 == データファイルの保存方法が原因 == == データファイルの保存方法が原因 ==
行 56: 行 104:
 LaTeXで書いたテキストをdvipdfmxでPDFに変換すると、なぜか「~」の記号のフォントが、置き換えられてしまいます。見た目はほとんど変わりませんが、半角の「~」では無くなっています。改行コードが余計につくこともあります。 LaTeXで書いたテキストをdvipdfmxでPDFに変換すると、なぜか「~」の記号のフォントが、置き換えられてしまいます。見た目はほとんど変わりませんが、半角の「~」では無くなっています。改行コードが余計につくこともあります。
 テキストのコードをそのままRに貼り付けると、プログラムとして認識できない文字が入力されたことになりますので、「~」は半角文字で打ち直すようにして下さい。 テキストのコードをそのままRに貼り付けると、プログラムとして認識できない文字が入力されたことになりますので、「~」は半角文字で打ち直すようにして下さい。
-問題なく動作するコードを下に記します。+あるいはPDFファイルではなく、このページからコピー＆ペーストして下さい。 
 + 
 +課題1について、問題なく動作するコードを下に記します。
  
 <code> <code>
行 92: 行 142:
 </code> </code>
  
-==== 演習課題1 ====+==== 演習課題1 2012.06.14 ====
 |出題|2012.06.14| |出題|2012.06.14|
 |〆切|2012.06.28| |〆切|2012.06.28|
 +
 +レポートの提出について
 +
 +|提出先|mailto:watalu.yamamoto_at_gmail.com (_at_を＠に交換してください)|
 +|提出方法|PDFファイルに変換してメールに添付して送付してください|
 +|提出確認|たまにまとめて返信します|
 +
 +
  
 === データの説明 === === データの説明 ===
行 278: 行 336:
 このデータの再来院に関して、以上の分析から、何が言えるか、まとめよ。 このデータの再来院に関して、以上の分析から、何が言えるか、まとめよ。
  
-==== 演習課題2 ====+==== 演習課題2 2012.07.26 ==== 
 +=== 日程 === 
 +|出題|2012.07.26| 
 +|〆切|2012.08.10| 
 + 
 === 今回の解析の流れ === === 今回の解析の流れ ===
  
行 346: 行 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",
行 352: 行 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"), 
行 364: 行 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>
 === まずは不具合発見数の成長を眺める === === まずは不具合発見数の成長を眺める ===
  
行 496: 行 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本なので、ロバストな推定方法は使わない。
行 604: 行 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>
行 625: 行 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>
行 645: 行 728:
 直近の改修よりも、以前の改修の方が、不具合の発生に大きく影響を与える、という母数の推定値が得られる。 直近の改修よりも、以前の改修の方が、不具合の発生に大きく影響を与える、という母数の推定値が得られる。
 これは、初期に多くの改修を行ったことと、無関係ではないかもしれない。 これは、初期に多くの改修を行ったことと、無関係ではないかもしれない。
 +
 +以下はbetaとbeta.2の切り替え時点を300から500に変更してみる場合の強度関数である。
  
 <code> <code>
行 657: 行 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の間の強度は
行 731: 行 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>