Nothing
`TMC` <-
function(num.Iter=50, data.Size=100, make.Data = gen.Data, make.Params = gen.Params,
model.List, weight.Vector = rep(1, times = length(model.List)), msc.List,
fit.Model = fit.Models, stepSize = .05, sumstats = list("Median Rank" = median ), huge=FALSE,
var.Frame=data.frame(), par.Sigma=1, data.Sigma=1, barebones=FALSE, allow.Negs = FALSE,
thresholds = c(1,2,3,5,10), test.Size = 0, scale.Frame = TRUE, use.Ranks = TRUE, ...){
ptm <- proc.time()
# Generates a matrix, each row of which contains weights which together form
# a convex combination of the MSCs in msc.List. convCombMat2 is similar,
# but allows for more general linear combinations, including those with negative
# weights, and the weights do not necessarily sum to unity.
convCombMat <- function(numFactors = 3, stepSize = stepSize, tol = 10^-9, huge = FALSE){
if(numFactors == 1) stop("Need more than one MSC")
if((round(1/stepSize))^numFactors > 500000){if(!huge)
{stop("Matrix will be BIG... set huge = TRUE to go on anyway.")}}
tol <- 10^-9
big <- seq(0,1,by=stepSize)
l = length(big)
mat <- matrix(nrow=l^(numFactors-1), ncol = numFactors)
if(numFactors == 2){mat[,1] = big
mat[,2] = 1 - big
return(mat)
}
for (i in 1:(numFactors-1)){
mat[,i] <- rep(big, times = l^(i-1), each = l^(numFactors - 1 - i))
}
mat <- zapsmall(mat)
mat[,numFactors] <- 1 - apply(mat[,1:(numFactors-1)], 1, sum)
zapsmall(mat[mat[,numFactors] >= 0,])
}
convCombMat2 <- function(numFactors = 3, stepSize = stepSize, tol = 10^-9, huge = FALSE){
if(numFactors == 1) stop("Need more than one MSC")
if((round(1/stepSize))^numFactors > 500000){if(!huge)
{stop("Matrix will be BIG... set huge = TRUE to go on anyway.")}}
tol <- 10^-9
big <- seq(-1,1,by=stepSize)
l = length(big)
mat <- matrix(nrow=l^(numFactors), ncol = numFactors)
if(numFactors == 2){mat[,1] = big
mat[,2] = 1 - big
return(mat)
}
for (i in 1:(numFactors)){
mat[,i] <- rep(big, times = l^(i-1), each = l^(numFactors - i))
}
mat <- zapsmall(mat)
mat
}
if(!is.numeric(test.Size) || test.Size != round(test.Size) || test.Size < 0) stop("test.Size must be a nonnegative integer")
if(!is.list(msc.List))stop("msc.List must be a list")
if(!is.list(model.List))stop("model.List must be a list")
if(!is.list(sumstats))stop("sumstats must be a list")
if(any(weight.Vector < 0) || !is.numeric(weight.Vector)) stop("weight.Vector must be a vector of non-negative numbers")
if(!is.function(make.Data))stop("make.Data must be a function")
if(!is.function(make.Params))stop("make.Params must be a function")
if(!is.function(fit.Model))stop("fit.Model must be a function")
if(!is.data.frame(var.Frame))stop("var.Frame must be a data frame")
if(barebones) {
sumstats <- as.list(rep("Internal Function", times = length(thresholds)+1))
names(sumstats) <- c("Average Rank", paste("P(Rank > ", thresholds, ")", sep=''))
}
if(!(deparse(substitute(var.Frame)) %in% search()))attach(var.Frame)
on.exit(try(detach("var.Frame"), silent = TRUE))
if(scale.Frame) var.Frame <- data.frame(scale(var.Frame))
# useful later in plotting functions if all elements of msc.List are named
namevec <- names(msc.List)
if(is.null(namevec))namevec <- rep("", length(msc.List))
namevec[namevec == ""] <- sapply(match.call()$msc.List, deparse)[-1][namevec == ""]
names(msc.List) = namevec
if(("holdout.Mean" %in% names(msc.List) || "holdout.Med" %in% names(msc.List) || "holdout.SS" %in% names(msc.List))
&& test.Size < 1) stop("Cannot use holdout msc functions if test.Size is not positive!")
test.Vector <- numeric()
errs <- rep(0, num.Iter)
weight.Vector <- weight.Vector/sum(weight.Vector)
num.Crits <- length(msc.List)
ifelse(allow.Negs, crits.Mat <- convCombMat2(num.Crits, stepSize=stepSize, huge=huge),
crits.Mat <- convCombMat(num.Crits, stepSize=stepSize, huge=huge))
num.Stats <- ifelse(!barebones, length(sumstats), length(thresholds))
simulated.Models <- list()
simulated.Parameters <- list()
simulated.Data <- list()
simulation.Attempts <- rep(0,num.Iter)
# These will ultimately be parts of the returned msc object.
if(!barebones) ranks.Mat <- matrix(nrow = dim(crits.Mat)[1], ncol = num.Iter)
Avg.Rank <- rep(0, dim(crits.Mat)[1]); Probs <- rep(0,dim(crits.Mat)[1]*length(thresholds))
dim(Probs) <- c(dim(crits.Mat)[1], length(thresholds))
Sum.Stats <- matrix(nrow = dim(crits.Mat)[1], ncol = num.Crits+num.Stats)
for(i in 1:num.Iter){
# in case we need a holdout sample
if(length(var.Frame) != 0 ){
if(test.Size != 0){train.Inds <- 1:(dim(var.Frame)[1] - test.Size)
test.Inds <- (dim(var.Frame)[1] - test.Size + 1):dim(var.Frame)[1]
test.Frame <- data.frame(var.Frame[test.Inds,])
train.Frame <- data.frame(var.Frame[train.Inds,])
}
else {train.Frame <- var.Frame; test.Frame <- NULL}
}
else {test.Frame = NULL; train.Frame = NULL}
success.Sim = FALSE
# Sometimes time-series simulation screws up
while(success.Sim == FALSE){
rand <- runif(1)
# Pick a random model from model.List using weight.Vector
model.Index <- which.max(cumsum(weight.Vector)>rand)
true.Model <- model.List[[model.Index]]
try(pars <- make.Params(true.Model, par.Sigma=par.Sigma, var.Frame = var.Frame, ...),
silent=TRUE)
try({dat <- make.Data(fmla=true.Model, pars=pars,
data.Size=data.Size, var.Frame = var.Frame, data.Sigma=data.Sigma,
train.Frame = train.Frame, test.Size = test.Size, ...);
data.Vector <- dat$data.Vector;
test.Vector <- dat$test.Vector}, silent=TRUE )
simulation.Attempts[i] <- simulation.Attempts[i]+1
if(exists("data.Vector") && is.numeric(data.Vector) && sum(!is.na(data.Vector))>0) success.Sim = TRUE
if(simulation.Attempts[i] > 100)stop("Simulation is screwing up too much")
}
#update records to include in final object
simulated.Models[[i]] <- true.Model
simulated.Parameters[[i]] <- pars
simulated.Data[[i]] <- data.Vector
if(length(train.Frame) != 0) train.Frame$y <- data.Vector
if(length(var.Frame) != 0) var.Frame$y <- c(data.Vector, test.Vector)
# necessary annoyance for calculating Cp later
if("Cp" %in% names(msc.List))
S2 <- (summary(lm(y ~ ., data = train.Frame))$sigma)^2
# Fit all models in model.List to the newly generated data.
fit <- function(x) try(fit.Model(fmla = x, c(data.Vector, test.Vector), var.Frame), silent=TRUE)
fit2 <- function(x) try(fit.Model(fmla = x, data.Vector, train.Frame), silent=TRUE)
fitted.Model.List <- lapply(model.List, fit)
if(test.Size > 0){
train.Model.List <- lapply(model.List, fit2)
for(m in 1:length(model.List))
{
fitted.Model.List[[m]] = list(full = fitted.Model.List[[m]], train = train.Model.List[[m]],
test.Vector = test.Vector, test.Frame = test.Frame)
if(exists("S2")) fitted.Model.List[[m]]$S2 = S2
class(fitted.Model.List[[m]]) <- "fmo"
}
}
else for(m in 1:length(model.List))
{
fitted.Model.List[[m]] = list(full = fitted.Model.List[[m]], train = NULL, test.Vector = NULL,
test.Frame = NULL)
if(exists("S2")) fitted.Model.List[[m]]$S2 = S2
class(fitted.Model.List[[m]]) <- "fmo"
}
err.Fits <- sapply(fitted.Model.List, function(x) "try-error" %in% c(class(x$full), class(x$train)))
# Sometimes fitting the true model to the simulated data gives an error... this is bad. Just throw this iteration out.
if(err.Fits[model.Index]==TRUE){
if(!barebones){
ranks.Mat[,i] <- rep(0, dim(ranks.Mat)[1])
next}
else {errs[i] <- 1; next}}
# update the model list to include only models where the fit actually worked.
# also change the index of the true model to reflect this updated list.
fitted.Model.List <- fitted.Model.List[!err.Fits]
model.Index <- which.max(sapply(model.List[!err.Fits],function(x)identical(x,true.Model)))
# Whew, what a line! Just apply all the msc functions in msc.List
# to each fitted model, and then either take ranks of the resulting vectors
# or standardize them, depending on whether use.Ranks is true.
ifelse(use.Ranks, crit.Vals <- apply(t(sapply(fitted.Model.List, function(x) sapply(msc.List,
function(y)do.call(y,list(x))))), 2, rank), crit.Vals <- apply(t(sapply(fitted.Model.List,
function(x) sapply(msc.List, function(y)do.call(y,list(x))))), 2, scale))
# update ranks.Mat if !barebones, else update the summary stats dynamically and do not retain ranks
# jittering breaks ties, which can make the results a bit discontinuous
temp <- tcrossprod(crits.Mat, crit.Vals)
temp2 <- numeric(dim(temp)[1])
for(j in 1:nrow(temp)) {
if(!barebones) ranks.Mat[j,i] <- sum(temp[j,] <= temp[j,model.Index])
temp2[j] <- sum(temp[j,] <= temp[j,model.Index])
}
Avg.Rank <- ((i-1)/i)*Avg.Rank + (1/i)*temp2
for(j in 1:length(thresholds))
Probs[,j] = ((i-1)/i) * Probs[,j] + (1/i) * (temp2 > thresholds[j])
rm(temp2)
rm(data.Vector, true.Model, pars, fitted.Model.List, model.Index, err.Fits)
}
# Min rank of zero means something screwed up that iteration, so we should toss it out
colnames(Probs) <- paste("P(Rank > ", thresholds, ")", sep='')
colnames(crits.Mat) <- names(msc.List)
if(!barebones) errs <- apply(ranks.Mat, 2, min)==0
# annoying here... if sumstats is only length one, then the result of the "apply" is a vector, rather
# than a matrix, so transposing it gives the wrong dimensions. Hence we need these annoying if conditions.
if(!barebones){
if(length(sumstats)==0) Sum.Stats <- data.frame(cbind(crits.Mat, Avg.Rank, Probs))
if(length(sumstats)==1) Sum.Stats <- data.frame(cbind(crits.Mat, apply(ranks.Mat[,!errs],1,
function(x)sapply(sumstats,function(y)do.call(y,list(x)))), Avg.Rank, Probs))
if(length(sumstats)>1) Sum.Stats <- data.frame(cbind(crits.Mat, t(apply(ranks.Mat[,!errs],1,
function(x)sapply(sumstats,function(y)do.call(y,list(x))))), Avg.Rank, Probs))
colnames(Sum.Stats) <- c(colnames(crits.Mat), names(sumstats), "Average Rank", colnames(Probs))}
else {Sum.Stats <- data.frame(cbind(crits.Mat, Avg.Rank, Probs))
colnames(Sum.Stats) <- c(colnames(crits.Mat), "Average Rank", colnames(Probs))}
if(!barebones){
summary.Functions = as.list(rep("Internal Function", times = 1 + length(thresholds)))
names(summary.Functions) <- c("Average Rank", colnames(Probs))
summary.Functions <- c(sumstats, summary.Functions)
}
if("var.Frame" %in% search()) detach(var.Frame)
var.Frame$y = NULL
if(barebones){
temp <- list(call = match.call(), Sum.Stats = Sum.Stats,
data.Size = data.Size, num.Iter=num.Iter, msc.List=msc.List, summary.Functions = sumstats,
var.Frame = var.Frame, model.List = model.List,
num.Errors = sum(errs), error.Iterations = which(errs==TRUE),
weight.Vector = weight.Vector, make.Data=deparse(substitute(make.Data)),
stepSize=stepSize,make.Params = deparse(substitute(make.Params)),
fit.Model = deparse(substitute(fit.Model)), simulated.Models = simulated.Models,
simulation.Attempts = simulation.Attempts, time.Taken = proc.time() - ptm,
thresholds = thresholds, par.Sigma = par.Sigma, data.Sigma = data.Sigma,
test.Size = test.Size)
class(temp) <- c( "barebones","msc")
return(temp)
}
else{
temp <- list(call = match.call(), Sum.Stats = Sum.Stats, ranks.Mat = ranks.Mat,
data.Size = data.Size, num.Iter=num.Iter, msc.List=msc.List, summary.Functions = summary.Functions,
var.Frame = var.Frame, model.List = model.List,
num.Errors = sum(errs), error.Iterations = which(errs==TRUE),
weight.Vector = weight.Vector, make.Data=deparse(substitute(make.Data)),
stepSize=stepSize,make.Params = deparse(substitute(make.Params)),
fit.Model = deparse(substitute(fit.Model)), simulated.Models = simulated.Models,
simulated.Parameters = simulated.Parameters, simulated.Data = simulated.Data,
simulation.Attempts = simulation.Attempts, time.Taken = proc.time() - ptm,
thresholds = thresholds, par.Sigma = par.Sigma, data.Sigma = data.Sigma,
test.Size = test.Size)
class(temp) <- "msc"
rm(simulated.Data, simulated.Models, simulated.Parameters)
if(!barebones) rm(ranks.Mat)
gc()
return(temp)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.