Instant Vectorization - mapply
> setwd("C:/Cauldron/garage/Volatility/Learn/DeliberatePractice/Ahas") |
Learning mapply
For some reason , through out my code, I had conveniently forgotten mappy and was complicating things unnecessary.
Here is my old code for Problem-6.5 from Yudi Pawitan’s book.
> y <- as.numeric(strsplit("8.85 9.40 9.18 8.70 7.53 6.43 5.85 4.73 3.98 3.50 3.10 2.80", " ")[[1]]) > x <- seq(0, 110, 10) > x <- x - mean(x) > X <- cbind(1, x) > lhd <- function(params) { beta0 <- params[1] beta1 <- params[2] sig2 <- params[3] mu <- beta0/(1 + exp(-beta1 * x)) ll <- -n/2 * log(sig2) - 1/(2 * sig2) * sum((y - mu)^2) return(-ll) } > fit <- nlm(lhd, c(0.5, 2, 0.5), hessian = TRUE) > betas <- fit$est[1:2] > H <- fit$hessian > se <- sqrt(diag(solve(H))) > mu <- betas[1]/(1 + exp(-betas[2] * x)) > lhd <- function(params, beta1) { beta0 <- params[1] sig2 <- params[2] mu <- beta0/(1 + exp(-beta1 * x)) ll <- -n/2 * log(sig2) - 1/(2 * sig2) * sum((y - mu)^2) return(-ll) } > betas1 <- seq(-0.1, 0, length = 100) > est <- sapply(betas1, function(b1) { nlm(lhd, c(12, 0.1), beta1 = b1, hessian = TRUE)$est }) > est <- t(est) > ll <- apply(cbind(est, betas1), 1, function(z) { lhd(z[1:2], z[3]) }) > ll <- exp(min(ll) - ll) > residuals <- mu - x |
Here is my new code for Problem-6.5 from Yudi Pawitan’s book.
> betas1 <- seq(-0.1, 0, length = 100) > est <- sapply(betas1, function(b1) { nlm(lhd, c(12, 0.1), beta1 = b1, hessian = TRUE)$est }) > est <- t(est) > ll <- mapply(lhd, est, betas1) > ll <- exp(min(ll) - ll) > residuals <- mu - x |
The difference is one line of statement
> ll <- apply(cbind(est, betas1), 1, function(z) { lhd(z[1:2], z[3]) }) > ll <- mapply(lhd, est, betas1) |
Instead of cbinding the columns and passing them to a function, I can use mappy and instantly vectorize