Nothing
errorModelSelection <- function(project=NULL, pen.coef=NULL, nb.model=1, f.min=0.001) {
if (!is.null(project))
mlx.loadProject(project)
go <- mlx.getObservationInformation()
go$name <- names(go$mapping)
go$prediction <- mlx.getContinuousObservationModel()$prediction
go$type <- go$type[1:length(go$name)]
ic <- which(go$type=="continuous")
y.name <- go$name[ic]
y.obs <- go[y.name]
y.pred <- go$prediction[y.name]
n.out <- length(ic)
pred <- getSimulatedPredictions()
# while( (sum(apply(pred[[1]],2,is.nan))+sum(apply(pred[[1]],2,is.na)) + sum(apply(pred[[1]],2,is.infinite))) > 0) {
# browser()
# runConditionalDistributionSampling()
# pred <- getSimulatedPredictions()
# }
# obs.model <- mlx.getContinuousObservationModel()
# obs.names <- mlx.getData()$observationNames
# fit.info <- mlx.getFit()
# for (noj in seq_along(obs.names)) {
# nojk <- grep(obs.names[noj],fit.info$data)
# if (length(nojk)>0)
# obs.names[noj] <- fit.info$model[nojk]
# }
# i.contObs <- which(mlx.getData()$observationTypes=="continuous")
# i.contModel <- which(names(obs.model$prediction) %in% obs.names[i.contObs])
# n.out <- length(i.contModel)
#d <- mlx.getObservationInformation()
if (nb.model==1) {
res.errorModel <- NULL
for (i.out in (1:n.out)) {
# name.predi <- obs.model$prediction[[i.contModel[i.out]]]
# name.obsi <- names(obs.model$prediction[i.contModel[i.out]])
# y.obsi <- d[[name.obsi]][[name.obsi]]
# y.predi <- pred[[name.predi]][[name.predi]]
y.obsi <- y.obs[[y.name[i.out]]][[y.name[i.out]]]
y.predi <- pred[[y.pred[i.out]]][[y.pred[i.out]]]
test <- T
kt <- 0
while(test | kt<10) {
kt <- kt+1
resi <- tryCatch(computeBIC(y.obs=y.obsi, y.pred=y.predi, pen.coef=pen.coef[i.out], nb.model=nb.model, f.min=f.min),
error = function(e) e)
if (inherits(resi, "error")) {
print("error !!!")
mlx.saveProject("foo.mlxtran")
mlx.loadProject("foo.mlxtran")
pred <- getSimulatedPredictions()
y.predi <- pred[[y.pred[i.out]]][[y.pred[i.out]]]
} else {
test <- F
}
}
res.errorModel <- c(res.errorModel, as.character(resi[['error.model']]))
names(res.errorModel)[i.out] <- y.name[i.out]
}
} else {
res.errorModel <- list()
for (i.out in (1:n.out)) {
# name.predi <- obs.model$prediction[[i.contModel[i.out]]]
# name.obsi <- names(obs.model$prediction[i.contModel[i.out]])
# y.obsi <- d[[name.obsi]][[name.obsi]]
# y.predi <- pred[[name.predi]][[name.predi]]
y.obsi <- y.obs[[y.name[i.out]]][[y.name[i.out]]]
y.predi <- pred[[y.pred[i.out]]][[y.pred[i.out]]]
resi <- tryCatch(computeBIC(y.obs=y.obsi, y.pred=y.predi, pen.coef=pen.coef[i.out], nb.model=nb.model, f.min=f.min),
error = function(e) e)
if (inherits(resi, "error")) {
print("error !!!")
mlx.saveProject("foo.mlxtran")
mlx.loadProject("foo.mlxtran")
pred <- getSimulatedPredictions()
y.predi <- pred[[y.pred[i.out]]][[y.pred[i.out]]]
resi <- computeBIC(y.obs=y.obsi, y.pred=y.predi, pen.coef=pen.coef[i.out], nb.model=nb.model, f.min=f.min)
}
res.errorModel[[i.out]] <- resi
names(res.errorModel)[i.out] <- y.name[i.out]
}
}
return(res.errorModel)
}
e.min1 <- function(x,y.pred,y.obs) {
sigma2 <- (x[1]+x[2]*y.pred)^2
e <- sum((y.obs-y.pred)^2/sigma2) + sum(log(sigma2))
return(e)
}
e.min2 <- function(x,y.pred,y.obs) {
sigma2 <- x[1]^2+(x[2]*y.pred)^2
e <- sum((y.obs-y.pred)^2/sigma2) + sum(log(sigma2))
return(e)
}
computeBIC <- function(y.obs, y.pred, pen.coef, nb.model, f.min) {
y.pos <- (min(y.obs) >= 0)
nrep <- length(y.pred)/length(y.obs)
y.obs <- rep(y.obs, nrep)
if (y.pos) {
iy <- (y.pred>0)
y.obs <- y.obs[iy]
y.pred <- y.pred[iy]
}
yb <- mean(abs(y.obs))
# if (sum(is.na(y.pred))+sum(is.nan(y.pred))>0)
# browser()
a.cons <- sqrt(mean((y.obs-y.pred)^2, na.rm=T))
x.min <- nlm(e.min2,c(a.cons/10,0.2),y.pred,y.obs)
a.comb2 <- abs(x.min$estimate[1])
b.comb2 <- abs(x.min$estimate[2])
x.min <- nlm(e.min1,c(a.comb2,b.comb2),y.pred,y.obs)
a.comb1 <- x.min$estimate[1]
b.comb1 <- x.min$estimate[2]
x.min <- nlm(e.min2,c(a.comb1,b.comb1),y.pred,y.obs)
a.comb2 <- abs(x.min$estimate[1])
b.comb2 <- abs(x.min$estimate[2])
if (a.comb2/yb < b.comb2*f.min)
a.comb2 <- 0
else if (b.comb2 < a.comb2/yb*f.min)
b.comb2 <- 0
if (a.comb1/yb < b.comb1*f.min)
a.comb1 <- 0
else if (b.comb1 < a.comb1/yb*f.min)
b.comb1 <- 0
if (y.pos == TRUE ) {
b.prop <- sqrt(mean(((y.obs/y.pred-1)^2)))
a.expo <- sqrt(mean((log(y.obs)-log(y.pred))^2))
error.model=c("constant", "proportional", "combined1","combined2","exponential")
# sigma2 <- cbind(a.cons^2, (b.prop*y.pred)^2, (a.comb1+b.comb1*y.pred)^2,
# a.comb2^2+(b.comb2*y.pred)^2, a.expo^2)
sigma2 <- cbind(a.cons^2, (b.comb2*y.pred)^2, (a.comb1+b.comb1*y.pred)^2,
a.comb2^2+(b.comb2*y.pred)^2, a.expo^2)
df <- c(1, 1, 2, 2, 1)
} else {
error.model=c("constant","combined1","combined2")
sigma2 <- cbind(a.cons^2, (a.comb1+b.comb1*y.pred)^2, a.comb2^2+(b.comb2*y.pred)^2)
df <- c(1, 2, 2)
}
ll <- pen <- bic <- NULL
for (k in 1:length(error.model)) {
if (error.model[k] == "exponential")
ll <- c(ll, - 0.5*sum( (log(y.obs)-log(y.pred))^2/(a.expo^2) + log(2*pi*(a.expo^2)) + 2*log(y.obs) )/nrep)
else
ll <- c(ll, - 0.5*sum( (y.obs-y.pred)^2/sigma2[,k] + log(2*pi*sigma2[,k]) )/nrep)
}
bic <- -2*ll + pen.coef*df
E <- data.frame(error.model=error.model,
ll=ll, df=df, criterion= bic)
nb.model <- min(nb.model, length(bic))
E <- E[order(bic)[1:nb.model],]
row.names(E) <- 1:nrow(E)
return(E)
}
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.