R/mixARsim.R

In GeoBosh/mixAR: Mixture Autoregressive Models

## 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)
# }