## Do not edit this file manually.
## It has been automatically generated from mixARsim.org.
simuExperiment <- function(model, simu, est,
N = 100, use_true = FALSE, raw = FALSE, init_name = "init",
keep = identity, summary_fun = .fsummary,
...){ # todo: arguments subject to change!
# maybe ... should be for 'sim' or 'est'.
# todo: zasega predavam '...' na .fsummary
fsimu <- match.fun(simu$fun)
fsimu_args <- simu$args # todo: kakvo ste stane ako nyakoy e slozhil '...' tuk?
fest <- match.fun(est$fun)
festargs <- est$args
usepartrue <- use_true
fkeep <- keep
fsummary <- summary_fun
if(usepartrue){
festargstrue <- festargs
festargstrue[[init_name]] <- model # todo: krapka! was: partrue - but what was the
} # intent? there is no `partrue' variable.
wrk <- vector(mode = "list", N)
currpartrue <- list() # in case usepartrue is FALSE when this would be left undefined
# but does this matter, there is an `if' checking this below?
for(i in 1:N){
data <- do.call(fsimu, fsimu_args) # generate data
if(usepartrue) # use the true parameters as starting values.
currpartrue <- do.call(fest, festargstrue)
currpar <- do.call(fest, c(list(data), festargs))
wrk[[i]] <- if(usepartrue) fkeep(currpar, currpartrue)
else fkeep(currpar)
}
res <- list(Summary = fsummary(wrk, ...))
if(raw)
res$Raw <- wrk
res
}
# based on simG2a
mixARExperiment <- function(model, imodel= NULL, simargs = NULL, estargs = NULL, fix, ...){
sa <- list(n = 200, init = rep(0, max(model@arcoef@p) ))
if(!is.null(simargs))
sa[names(simargs)] <- simargs
sa <- c(list(model), sa)
## 2020-06-12 was: cat("sa$init is: ", sa$init, "\n")
ea <- list(crit = 10^(-4))
if(!missing(fix))
ea$fix <- fix
if(!is.null(estargs))
ea[names(estargs)] <- estargs
if(is.null(imodel))
ea <- c(list(model), ea)
else
ea <- c(list(imodel), ea)
keep <- function(x) x$model
fsum <- .fsumA # 2018-11-24 was: mixARsim:::.fsumA
sim <- list(fun = "mixAR_sim", args = sa) # argsWL_II$sim
est <- list(fun = "fit_mixAR", args = ea) # argsWL_II$est
simuExperiment(TRUE, simu = sim
, est = est
, keep = keep
, summary_fun = fsum
, ...)
}
.fsummary <- function(obj, relist = TRUE, merge = FALSE){
skel <- as.relistable(obj[[1]])
wrk <- sapply(obj, function(x) unlist(x))
flag <- !is.null(dim(wrk)) # res <- rowMeans(wrk)
res <- list( mean = if(flag) apply(wrk, 1, mean) else mean(wrk)
, sd = if(flag) apply(wrk, 1, sd) else sd(wrk)
)
if(relist)
res <- lapply(res, function(x) relist(x, skel))
if(merge)
res$merged <- sapply(res, identity)
res$N <- length(obj)
res
}
.fsumA <- function(models, relabel = TRUE){ # todo: add argument 'fstat' - a list of
# functions to apply
pat <- models[[1]]
params <- sapply(models, parameters)
# 2018-11-24 changed set_parameters() below to "parameters<-()".
# The change is mechnical, to avoid errors from R check.
# TODO: make this better.
# resOld <- list( mean = set_parameters(pat, apply(params, 1, mean))
# , sd = set_parameters(pat, apply(params, 1, sd ))
# )
res <- list( mean = {parameters(pat) <- apply(params, 1, mean); pat}
, sd = {parameters(pat) <- apply(params, 1, sd ); pat}
)
# stopifnot(identical(resOld$mean, res$mean))
# stopifnot(identical(resOld$sd, res$sd ))
if(relabel){
usmodels <- unswitch(models, pat)
usparams <- sapply(usmodels, parameters)
# resOld <- c(resOld, list(
# meanUnswitched = set_parameters(pat, apply(usparams, 1, mean))
# , sdUnswitched = set_parameters(pat, apply(usparams, 1, sd ))
# ))
if(identical(usparams, list()))
res <- c(res, list(meanUnswitched = NA, sdUnswitched = NA))
else
res <- c(res, list( meanUnswitched = {parameters(pat) <- apply(usparams, 1, mean); pat}
, sdUnswitched = {parameters(pat) <- apply(usparams, 1, sd ); pat}
))
# stopifnot(identical(resOld$meanUnswitched, res$meanUnswitched))
# stopifnot(identical(resOld$sdUnswitched, res$sdUnswitched ))
}
# stopifnot(identical(resOld, res))
res
}
permn_cols <- function(m){
permn(ncol(m), function(j) list(j, m[,j]))
}
unswitch <- function(models, true_model, Nref = 100, simargs = NULL){
sa <- list(n = Nref, init = numeric(10))
if(!is.null(simargs))
sa[names(simargs)] <- simargs
sa <- c(list(true_model), sa)
refts <- do.call("mixAR_sim", sa)
reftau <- mixAR_cond_probs(true_model, refts) # todo: check if "representative"? (and
reftau <- reftau@m # posibly simulate another ts)
refranks <- t(apply(reftau, 1, rank))
# lapply(models, unswitch_1, reftau, refts)
allperm <- matrix(0, nrow = length(models), ncol = ncol(reftau))
for(i in seq_along(models)){
f <- function(x){
tf <- x == refranks
sum( apply(tf, 1, all) ) # simply counts matching rows; todo: improve!
}
tau <- mixAR_cond_probs(models[[i]], refts)
tau <- tau@m
permtau <- permn_cols(tau) # todo: smyatay permutatsiyata vednazh! !!!
ranks <- lapply(permtau, function(x) t(apply(x[[2]], 1, rank))) # x[[2]] is the matrix
wrk <- sapply(ranks, f)
## 2020-06-12 was: print(wrk)
# if(sum(wrk)< nrow(reftau)) # todo: if sigma_k for some component is (close
# browser() # to) 0, bad things happen (all taus are NaN).
#
# 2011-12-06 - a new patch in mixAR_cond_probs,
# todo: check if it reduces the number of
# problem cases (maybe the number quoted as
# "Dropping XXX unsuccesful ..."
best <- if(sum(wrk) == nrow(reftau)) # todo: krapka!
which.max(wrk)
else # if this kicks in, most of wrk and tau is NaN.
NA
allperm[i, ] <- if(is.na(best)) NA
else permtau[[c(best,1)]] # save the "best" permutaion
}
# browser()
good <- !is.na(allperm[,1])
## 2020-06-12 was: if(sum(good) < length(models))
## cat("Dropping", length(models) - sum(good), "unsuccessful runs.\n")
models <- models[good]
allperm <- allperm[good,] # a matrix containing the best permutation of models[[i]]
# in its i-th row (not true).
fp <- function(i){
mixAR_permute(models[[i]], allperm[i,])
}
res <- lapply(seq_along(models), fp) # todo: message if length(models)==0 ?
res
}
test_unswitch <- function(models, true_model, allperm, ...){ # todo: introduce weights (for the call to
# sample.int)
if(missing(allperm))
allperm <- permn( length(models[[1]]@prob) )
np <- length(allperm)
f <- function(model){
curp <- allperm[[ sample.int(np,1) ]] # permute randomly;
mixAR_switch(model, curp)
}
wrk <- lapply(models, f)
wrkall <- list(original = models,
u_original = unswitch(models, true_model, ...),
messed = wrk,
u_messed = unswitch(wrk, true_model, ...)
)
# 2018-11-24 was: mixARsim:::.fsumA(x, relabel = FALSE))
summaries <- lapply(wrkall, function(x) .fsumA(x, relabel = FALSE))
list(mean = sapply(summaries, function(x) parameters(x$mean, namesflag = TRUE)),
sd = sapply(summaries, function(x) parameters(x$sd, namesflag = TRUE)),
summaries = summaries,
Raw = wrkall
)
}
# par(mfrow = c(1,3))
# for (i in 1:3){
# title <- substitute(gamma[j], list(j=i))
# plot(1:10, main = title)
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.