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