R/em0.R

In GeoBosh/mixAR: Mixture Autoregressive Models

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

em_tau <- function(stdetk, prob, scale, pdf = dnorm){  # em_tau was previously named 'etk2tau'
    ## assumes etk is standardised and mixComp
    wrk <- prob * (pdf %of% stdetk) / scale
                                                 # todo: why [] below? to keep the dimensions?
    wrk@m[] <- wrk@m / inner(wrk)         # need to work on wrk@m directly since `/'
    wrk                                   #     works differently on mixComp objects
}

em_tau_safe <- function(stdetk, prob, scale, pdf = dnorm){
    nanflag <- any(!is.finite(stdetk))
    if(nanflag)                                 # todo: krapka! appropriate if matching
        stdetk@m[!is.finite(stdetk)] <- Inf     #       p_k, scale_k sre both (close to) zero.
                                                # expects that pdf(Inf) = 0
    wrk <- prob * (pdf %of% stdetk) / scale

    if(nanflag)                            # this makes the above if(nanflag) replacement
        wrk@m[!is.finite(stdetk)] <- 0     # superfluous but that helps to avoid some error
                                           # messages.

    wrk@m[] <- wrk@m / inner(wrk)         # need to work on wrk@m directly since `/'
                                          #   works differently on mixComp objects
    wrk
}

                                          # todo:  see the Mathematica programs for extension.
tau2probhat <- function(tau){          # Note: this does not depend on the noise distribution.
    colSums(tau@m)/sum(tau@m)
}

etk2tau <- function(etk){                               # see my Mathematica function EstepMix
    f <- function(x){
        ax <- abs(x)
        val <- 1 - ax/sum(ax)
        val
    }
                     #    wrk <- apply(etk, 2, f) - tova izglezhda da e v Mathematica.
                     # Ili ne si razchitam Mathematica programata kakto tryabva ili e nesto
                     # gnilo, promenyam kakto mi se struva che sam iskal da bade.
    wrk <- t( apply(etk, 1, f) )
    if(!all(is.finite(wrk)))
        wrk[!is.finite(wrk)] <- 0

    res <- wrk / rowSums(wrk)   # todo: zastita srestu delenie na nula.
    if(!all(is.finite(res)))
        res[!is.finite(res)] <- 0

    # browser()

    res
}

est_templ <- function(model, shift = TRUE, ...){
    if(length(shift) == 1)
        shift <- rep(shift, .nmix(model))
    f <- function(k){
        s <- if(shift[k]) NA
             else         model@shift[k]
        co <- model@arcoef@a[[k]]
        co[] <- NA                    # todo: curently always estimate all coef.
        list(s, co)
    }

    lapply(1:.nmix(model), f)
}