R/emGaussian.R

In GeoBosh/mixAR: Mixture Autoregressive Models

## Do not edit this file manually.
## It has been automatically generated from *.org sources.

tauetk2sigmahat <- function(tau, etk){                           # this assumes Gaussian noise
    wrk <- colSums(tau@m * etk@m^2) / colSums(tau@m)     # note: the denominator is normally 1
                                                         # todo: proveri gornoto tvardenie!
    if(any(is.nan(wrk)))
        wrk[is.nan(wrk)] <- 0.0000001 # todo: this is a crude fix (KRAPKA)!

    sqrt(wrk)
}

tauCorrelate <- function(y, tau, order){
    p <- max(order)
    ta <- tau@m
    g   <- ncol(ta)
    nmp <- nrow(ta)
    n <- length(y)

    ind <- (n - nmp + 1):n

    Stau <- colSums(ta)

    Stauy <- matrix(0, nrow = g, ncol = 1 + p)
    Stauyy <- array(0, c(g, 1 + p, 1 + p))
    for(i in 0:p){
        tauy <- ta * y[ind - i]
        Stauy[, i + 1] <- colSums(tauy)
        for(j in 0:p){                    # todo: exploit the symmetry w.r.t i and j
            tauyy <- tauy * y[ind - j]
            Stauyy[, i + 1, j + 1] <- colSums(tauyy)
        }
    }
    list(Stau = Stau, Stauy = Stauy, Stauyy = Stauyy)
}

mixSubsolve <- function(k, pk, Stau, Stauy, Stauyy, shift, tol = 1e-7){
    A <- matrix(0, nrow = pk + 1, ncol = pk + 1)
    b <- c(Stauy[k, 0 + 1], Stauyy[k, seq_len(pk) + 1, 0 + 1])

    A[, 1] <- c(Stau[k], Stauy[k, seq_len(pk) + 1])
    for(j in seq_len(pk)){
        A[, j + 1] <- c(Stauy[k, j + 1], Stauyy[k, seq_len(pk) + 1, j + 1])
    }
    res <- if(shift)
               try(     solve(A, b))       # 2011-11-24 - svd to guard against singular A
           else      
               try(c(0,                          # TODO: this works if the shifts are zeroes!
                     ## 2021-08-09 handle dim(A) = c(0,0)
                     if(length(A) > 1) # so, at least 2x2
                         solve(A[-1, -1], b[-1])
                     else
                         numeric(0)
                     ))

    if(inherits(res, "try-error")){
        cat("mixSubsolve: singular system, trying again with svd\n")
        res <- if(shift) pseudoInverse(A, tol) %*% b
               else{
                   ## 2021-08-09 handle dim(A) = c(0,0)
                   c(0, 
                     if(length(A) > 1)
                         c(0, pseudoInverse(A[-1, -1], tol) %*% b[-1])
                     else
                         numeric(0)
                     )
               }
    }
    res
}

tau2arcoef <- function(y, tau, order, est_shift = TRUE){        # estimate arcoef, Gaussian case
    wrk <- tauCorrelate(y, tau, order)                          # attention: arcoef AND shift!

    g <- length(order)

    shift <- numeric(g)
    arcoef <- list()
    for(k in 1:g){
        wrk2 <- mixSubsolve(k, order[k], wrk$Stau, wrk$Stauy, wrk$Stauyy, est_shift)
        shift[k] <- wrk2[1]
        arcoef[[k]] <- wrk2[-1]
    }
    list(shift = shift, arcoef = arcoef)
}

mixMstep <- function(y, tau, order, index, est_shift = TRUE){
    prob <- tau2probhat(tau)

    co <- tau2arcoef(y, tau, order, est_shift = est_shift)
    arcoef <- new("raggedCoef", co$arcoef)
    shift <- co$shift
    scale <- rep(NA_real_, length(order))

    model <- new("MixARGaussian", prob = prob, order = order, shift = shift, scale = scale, 
                 arcoef = arcoef)

    etk <- mix_ek(model, y, index)
    sigmahat <- tauetk2sigmahat(tau, etk)
    model@scale <- sigmahat

    model
}

em_rinit <- function(y, order, partempl){
    etk <- randomMarResiduals(y, order, partempl)

    tau <- etk2tau(etk)
    tau <- new("MixComp", m = tau)

    indx <- (max(order) + 1) : length(y)

    est_shift <-                                              # krapka; todo: obmisli po-dobre
        if(is.logical(partempl))    # TRUE of FALSE in this case
            partempl
        else{
            wrk <- sapply(partempl, function(x) x[[1]])
            if(all(is.na(wrk)))
                TRUE
            else if(any(is.na(wrk)))
                stop("Currently mixMstep estimates either all or none of the 'shift's.")
            else
                FALSE
        }
    mixMstep(y, tau, order, indx, est_shift = est_shift)         # returns MixARGaussian model
}

## 2018-11-24: changing instead of devel_envir, use local()
mixARemFixedPoint <- local({em_global.res <- em_globalAll.res <- list(); function(y, model
                              , est_shift = TRUE
                              , crit = 1e-14
                              , maxniter = 200
                              , minniter = 10
                              , verbose = FALSE   # new: 2011-11-30; 2020-06-13: deprecated
                              ){
    trace <- if(isTRUE(verbose)  &&  interactive()) Inf else 0
        
    y <- as.numeric(y)
    oldmodel <- model

    pm <- max(oldmodel@order)
    n <- length(y)
    indx <- (pm + 1):n

    nmp <- n - pm

    oldvallogf <- cond_loglik(oldmodel, y)
    newvallogf <- NA
    relchange <- crit + 1

    if(trace > 0) cat("Initial vallogf: ", oldvallogf, "\n")

    special_flag <- FALSE
    emtrace <- list(ts = y, init = list(oldmodel, oldvallogf))

    niter <- 0
               # the check for newvallogf is a patch to get the simulations going (2011-11-22)
    while(!is.nan(newvallogf) && (niter <= minniter || niter < maxniter && relchange > crit)){
        niter <- niter + 1
                                                                                      # E-step
        oldetk <- mix_ek(oldmodel, y, indx, scale = TRUE)   # em_tau needs standardised resid.
        newtau <- em_tau(oldetk, oldmodel@prob, oldmodel@scale)

        newmodel <- mixMstep(y, newtau, oldmodel@order, indx, est_shift = est_shift)  # M-step

        oldvallogf <- newvallogf
        newvallogf <- cond_loglik(newmodel, y)

        ## todo: fot testing!
        ##       This records cases when the EM step results in decrease of the likelihood,
        ##       which, in principle shouldn't happen.
        if(is.finite(newvallogf) && is.finite(oldvallogf)  &&
                                                  newvallogf < oldvallogf  && !special_flag ){
            special_flag <- TRUE
            ## 2018-08-30 was:
            ##     if(exists("em_global.res")){
            ##         if(!exists("em_globalAll.res"))
            ##             em_globalAll.res <<- list()
            ##         em_globalAll.res[[length(em_globalAll.res) + 1]] <<- em_global.res
            ##     }
            ## 
            ## Move stored results, if any, from previous run to devel_envir$em_globalAll.res
            if(length(em_global.res) > 0){
                em_globalAll.res[[length(em_globalAll.res) + 1]] <<-
                    em_global.res
                em_global.res <<- list()
            }
            emtrace$niter <- niter
            emtrace[[length(emtrace) + 1]] <- list(oldmodel, oldvallogf)
        }

        if(special_flag )
            emtrace[[length(emtrace) + 1]] <- list(newmodel, newvallogf)

        ## 2020-06-08 was: verbose || niter %% 25 == 0
        if(niter %% 25 == 0  && trace > 0)                                # print some info
            cat("niter: ", niter, "\tvallogf: ", newvallogf, "\n")
        if(is.nan(newvallogf) && trace >= 1)
            cat("!!!! Log-likelihood is NaN, maybe due to singularity.\n")


        relchange <- abs((oldvallogf - newvallogf) / oldvallogf)          # stopping criterion

        oldmodel <- newmodel
    }

    if(trace > 0)
        cat("Final niter: ", niter, "\tvallogf: ", newvallogf, "\n")

    if(special_flag){                              # assign in devel_envir for further testing
        emtrace$niter <- c(emtrace$niter, niter)
        ## 2018-08-30 (see also the note above) was:
        ##    assign("em_global.res", emtrace, envir = .GlobalEnv)
        em_global.res <<- emtrace
    }

    list(model = newmodel, vallogf = newvallogf)
}})