差分
このページの2つのバージョン間の差分を表示します。
| 両方とも前のリビジョン前のリビジョン次のリビジョン | 前のリビジョン | ||
| dm:2011 [2011/08/01 13:58] – [準備] watalu | dm:2011 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1 | ||
|---|---|---|---|
| 行 282: | 行 282: | ||
| ==== 演習2 ==== | ==== 演習2 ==== | ||
| + | === 日程 === | ||
| + | |出題|2011.07.21| | ||
| + | |〆切|2011.08.04| | ||
| + | |||
| + | |||
| === 今回の解析の流れ === | === 今回の解析の流れ === | ||
| 行 553: | 行 558: | ||
| 以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。 | 以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。 | ||
| + | === モデルを当てはめる === | ||
| + | |||
| + | * このデータの場合、データ観測系列は1本なので、ロバストな推定方法は使わない。 | ||
| + | * 先週の説明の通り、飽和しかけているので、曲線を強度関数に用いる (累積データに「飽和する曲線」を当てはめる場合、変曲点後まで観測できているか否かで、当てはまりがとても異なるので注意) | ||
| + | |||
| + | データの先頭に1行、0を入れておく。 | ||
| + | |||
| + | < | ||
| + | software.debugging.diff <- rbind(c(0, | ||
| + | </ | ||
| + | |||
| + | 強度関数は，次の通り。 | ||
| + | |||
| + | < | ||
| + | lambda.t <- function(theta, | ||
| + | diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j])) | ||
| + | diff.2 <- (alpha+beta*sum(data$Ct.diff[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, | ||
| + | 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) | ||
| + | } | ||
| + | </ | ||
| + | |||
| + | 対数尤度の最大化は、上の関数の最小化と等しい。 | ||
| + | |||
| + | Rには最小化の関数が幾つかある。nlminb(), | ||
| + | |||
| + | これを最小化するのは、結構、困難であるが、例えば初期値を | ||
| + | |||
| + | < | ||
| + | fitted <- optim(c(0.001, | ||
| + | print(fitted) | ||
| + | </ | ||
| + | |||
| + | のように与えると、それなりの値に収束する。optim()の最初の引数は、最適化の初期値である。これは Cook and Lawless (2007)掲載の推定値を参考に設定した。 | ||
| + | 以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。 | ||
| + | |||
| + | == 補足 == | ||
| + | optimという関数は、目的関数と初期値を与えると最適化してくれる関数、で、幾つかの有名な最適化アルゴリズムを利用できて便利なのですが、探索範囲を制限することができない。 | ||
| + | |||
| + | 今回の演習で用いているのはNelder-Meed法という、Simplex法に似たアルゴリズムで、初期値から出発して、探索していきます。その際に、大きな歩幅で探索すると、warningが表示されるように計算きない点に到達する。そのため、今回のような「計算可能な探索範囲」が狭い最適化問題では、目的関数を評価できない旨のwarningが沢山、ログに残る。ただし、初期値の近傍の局所最適解に到達している場合には、その解が表示されるので、warningがある旨のメッセージではなく、最終的な結果が出力されているか否かで、最適化が行えたかどうか、判断して欲しい。 | ||
| + | (質問を受けたので補足) | ||
| === 強度関数を変えてみる === | === 強度関数を変えてみる === | ||