source("ols_lfn.R") source("simulated_setup.R") # There are really only two distinct likelihood functions: lf1() and # lf1.inlogs. The following test cases are desirable -- # BFGS L-BFGS-B Nelder-Mead SANN # lf1 a,n a,n yes yes # lf1.inlogs n NO yes yes # This function does measurement about the hit rate on reaching # the correct answer, and the time taken (in ms), of alternative MLE # algorithms. # You send in args specific to the method or situation using "..." measurement <- function(method, lf, trans=FALSE, ...) { N <- 100 # How many tests? success <- 0 costs <- system.time( for (i in 1:100) { seeds <- 10*runif(3) # random seeds from 0 to 10 if (trans) { seeds[1] <- log(seeds[1]) } optimised <- optim(seeds, lf, method=method, y=y, X=X, ...) if (trans) { optimised$par[1] <- exp(optimised$par[1]) } if (SSE(optimised$par) < 0.001) {success <- success + 1} } ) return(list(successrate=100*success/N, cost=1000*costs[1]/N)) } performance <- matrix(NA, nrow=9, ncol=2) colnames(performance) <- c("Hit rate", "Cost") rownames(performance) <- rep("Method = ", 9) i <- 1 # Start of 9 runs -- a <- measurement(method="BFGS", lf=ols.lf1) rownames(performance)[i] <- "BFGS, numerical" performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost i <- i + 1 a <- measurement(method="BFGS", lf=ols.lf1, gr=ols.gradient) rownames(performance)[i] <- "BFGS, analytical" performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost i <- i + 1 a <- measurement(method="L-BFGS-B", lf=ols.lf1, lower=c(1e-6,-Inf,-Inf), upper=rep(Inf,3)) rownames(performance)[i] <- "L-BFGS-B, numerical" performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost i <- i + 1 a <- measurement(method="L-BFGS-B", lf=ols.lf1, lower=c(1e-6,-Inf,-Inf), upper=rep(Inf,3), gr=ols.gradient) rownames(performance)[i] <- "L-BFGS-B, analytical" performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost i <- i + 1 a <- measurement(method="BFGS", lf=ols.lf1.inlogs, trans=TRUE) rownames(performance)[i] <- "BFGS, trans., numerical" performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost i <- i + 1 a <- measurement(method="Nelder-Mead", lf=ols.lf1) rownames(performance)[i] <- "Nelder-Mead" performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost i <- i + 1 a <- measurement(method="SANN", lf=ols.lf1) rownames(performance)[i] <- "SANN" performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost i <- i + 1 a <- measurement(method="Nelder-Mead", lf=ols.lf1.inlogs, trans=TRUE) rownames(performance)[i] <- "Nelder-Mead, trans." performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost i <- i + 1 a <- measurement(method="SANN", lf=ols.lf1.inlogs, trans=TRUE) rownames(performance)[i] <- "SANN, trans." performance[i,"Hit rate"] <- a$successrate; performance[i,"Cost"] <- a$cost i <- i + 1 # Print out summary statistics o <- order(performance[,2]) performance[o,]