2012-02-10 Rでトービット・モデルによる打ち切りデータの推定 の「2. optim()関数で最尤推定」で ifelse を使っているのは,ダメダメ。
掛け算の方がまだ速い。これで,元のプログラムより 25% 速くなった(元のプログラムが 0.104 sec. で,このプログラムが 0.079 sec. まあ,全然どうってことはないのです(^_^))
> system.time({
+ df <- read.table("kanazawa_accumulation_of_snow.txt",
+ header = TRUE, sep = "\t")
+ frml <- AS ~ TEMP
+ r_ols <- lm(frml, data = df)
+ p <- c(sqrt(deviance(r_ols) / df.residual(r_ols)), coefficients(r_ols))
+ X <- model.matrix(frml, data = df)
+ y <- df$AS
+ print(optim(p, function(p) {
+ sigma <- p[1]
+ xb <- X %*% p[-1]
+ -sum((y <= 0) * log(1 - pnorm(xb / sigma)) +
+ (y > 0) * log(dnorm((y - xb) / sigma) / sigma))
+ }, method = "BFGS", hessian = TRUE))
+ })
$par
(Intercept) TEMP
24.99719 135.28077 -17.73532
$value
[1] 170.5886
$counts
function gradient
105 99
$convergence
[1] 0
$message
NULL
$hessian
(Intercept) TEMP
0.118920177 0.007176688 0.05931874
(Intercept) 0.007176688 0.062764165 0.34602183
TEMP 0.059318737 0.346021832 2.11862555
ユーザ システム 経過
0.079 0.004 0.094
※コメント投稿者のブログIDはブログ作成者のみに通知されます