## Do not edit this file manually.
## It has been automatically generated from *.org sources.
mixgenMstep <- function(y, tau, model, index, fix = NULL # 2012-11-02 est_shift=TRUE
, comp_sigma = FALSE
, method = "BBsolve", maxit = 100, trace = FALSE, lessverbose = TRUE,
... ){
verbose <- lessverbose && interactive()
n <- length(y)
g <- .nmix(model)
order <- model@order
sigma <- model@scale
p <- order
pm <- max(order)
# 2012-11-02 tova beshe vav fit_mixAR
# ako ostavya taka, ste raboti bez izmenenie
shift_ar_params <- unlist(c(model@shift, model@arcoef@a))
## 2021-03-31 new
## TODO: this could be merged better with the code for the other cases
if(identical(fix, "white.noise")) {
## TODO: the calling function should check that the last component has order pmax !!!
white_noise <- TRUE
fix <- NULL # for now; could allow fixed params later!
constr_params <- numeric(g + sum(model@order))
constr_indx <- c(g, g + sum(model@order[1:(g-1)]) + 1:model@order[g])
unconstr_indx <- seq_along(constr_params)[-constr_indx]
}else
white_noise <- FALSE
dont_g <- .nmix(model)
dont_ar <- sum(model@order)
dontfix_shift <- if(identical(fix,"shift")) rep(FALSE, dont_g) else rep(TRUE, dont_g)
dontfix_scale <- rep(TRUE, dont_g)
dontfix_arcoef <- rep(TRUE, dont_ar)
if(is.list(fix)){ # should contain named logical vectors with TRUE for fixed elements
# dontfix_prob <- if(!is.null(fix$prob)) !fix$prob else rep(TRUE, dont_g)
if(!is.null(fix$shift)) dontfix_shift <- !fix$shift
if(!is.null(fix$scale)) dontfix_scale <- !fix$scale
if(!is.null(fix$arcoef)) dontfix_arcoef <- !fix$arcoef
}
dontfix_shift_ar <- c(dontfix_shift, dontfix_arcoef)
est_shift <- all(dontfix_shift_ar) # compatibility
est_ar_all <- all(dontfix_arcoef) && all(!dontfix_shift)
armaxit <- maxit # todo: make an argument for this?
dist <- get_edist(model) # dist <- model@dist
Fscore <- lapply( dist, function(x) x$Fscore ) # Fscore <- noise_dist(model, "Fscore")
fpdf <- lapply( dist, function(x) x$pdf ) # fpdf <- noise_dist(model, "pdf" )
# 2012-10-26 TODO: slagam vremenno na TRUE, za da testvam novite raboti
# todo: otchitay i novite raboti!
# 2012-10-31 estdist_flag <- TRUE # any( sapply(dist, function(x) x$any_param()) )
# 2021-03-31 was: estdist_flag <- with(environment(dist[[1]]$pdf), any(param_flags))
# it seems that I have changed the way this is retrieved at the time
#browser()
## 20211-06-26 TODO: this seems to assume that any_param() for the 1st component gives an
## overall result. Shouldn't it be any(sapply(dist, function(x) x$any_param)) ?
##
## For now leaving as is but this will not work correctly if the first component has no
## params to estimate but some of the others do.
estdist_flag <- dist[[1]]$any_param()
prob <- tau2probhat(tau) # TODO: new mixing weights, but model@prob is updated after the other parameters,
# i.e., prob is not used in their updates but rather model@prob.
# TODO: why? (seems ok!)
lcond <-
if(white_noise){ # 2021-03-31 new - modified from the 'else' part below
function(par){
# shift_ar_params[dontfix_shift_ar] <- par
# model@shift <- shift_ar_params[1:g]
# model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])
cur_probs <- model@prob # TODO: Dali tova ne tryava da e 'prob'?
# may ne, 'prob' e veroyatno initial value
shift_ar_params[unconstr_indx] <- par
model@shift <- c(shift_ar_params[1:(g-1)],
- sum(cur_probs[-g] * shift_ar_params[1:(g-1)]) / cur_probs[g])
# tmp_m <- matrix(shift_ar_params[-(1:g)], nrow = g - 1, ncol = pm)
tmp_m <- matrix(0, nrow = g - 1, ncol = pm)
tmp_v <- shift_ar_params[-(1:g)] # only ar param
for(i in 1:(g-1)){
r_ind <- 1:model@order[i]
tmp_m[i, r_ind] <- tmp_v[r_ind]
tmp_v <- tmp_v[-r_ind]
}
constr_ar_param <- numeric(pm)
for(k in 1:(g-1)) {
for(i in 1:pm){
constr_ar_param[i] <-
constr_ar_param[i] + cur_probs[k] * tmp_m[k, i]
}
}
shift_ar_params[constr_indx[-1]] <- - constr_ar_param / cur_probs[g]
model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])
stdetk <- mix_ek(model, y, index, scale=TRUE) # standardised resid.
sc <- tau * (log %of% ((fpdf %of% stdetk)/ sigma))
# tuk nyama nuzhda da vzemam podmnozhestvo ponezhe se
# vrasta samo edna stoynost!
sum(sc@m) + sum( (tau * log(model@prob))@m )
}
} else if(est_shift){
function(par){
model@shift <- par[1:g]
model@arcoef <- rag_modify(model@arcoef, par[-(1:g)])
stdetk <- mix_ek(model, y, index, scale=TRUE) # standardised resid.
sc <- tau * (log %of% ((fpdf %of% stdetk)/ sigma))
res <- sum(sc@m) + sum( (tau * log(model@prob))@m )
res
}
}else if(est_ar_all){
function(par){
model@arcoef <- rag_modify(model@arcoef, par)
stdetk <- mix_ek(model, y, index, scale=TRUE) # standardised resid.
log_etk <- log %of% ((fpdf %of% stdetk))
if(any(log_etk@m==-Inf)){
bad <- log_etk@m== -Inf
log_etk@m[bad] <- - 750 # todo: needs better thought!
# Arbitrary constant!
} # (but see R transcript above)
sc <- tau * log_etk
# if(any(is.infinite(sc))){
# bad <- sc==-Inf
# sc[bad] <- ifelse(tau[bad]<1e-10, 0, - 750 )
# }
res <- sum(sc@m) + sum( (tau * log(model@prob/sigma))@m )
if(any(is.infinite(res))){
## 2020-06-12 was: cat("infinite values in res!\n")
## browser()
if(verbose)
message("infinite values in res!")
}
res[res==-Inf] <- - 10000 # see the R snipped above
res <- sum(res)
if(res == -Inf){
## 2020-06-12 was: cat("res is: ", res, "\n")
## browser()
if(verbose)
message("res is: ", res, "\n")
}
res
}
}else{
function(par){
shift_ar_params[dontfix_shift_ar] <- par
model@shift <- shift_ar_params[1:g]
model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])
stdetk <- mix_ek(model, y, index, scale=TRUE) # standardised resid.
sc <- tau * (log %of% ((fpdf %of% stdetk)/ sigma))
# tuk nyama nuzhda da vzemam podmnozhestvo ponezhe se
# vrasta samo edna stoynost!
sum(sc@m) + sum( (tau * log(model@prob))@m )
}
}
eqns <- matrix(NA, nrow = g, ncol = 1 + pm)
feqns <-
if(est_shift){ # estimate all shift and ar parameters (ie all of them are not fixed)
function(par){
model@shift <- par[1:g]
model@arcoef <- rag_modify(model@arcoef, par[-(1:g)])
stdetk <- mix_ek(model, y, index, scale=TRUE) # standardised resid.
sc <- (tau * (Fscore %of% stdetk)) / sigma
# TODO: tova tyrabva da se vidi vnimatelno!
# 2012-10-31 smenyam znaka; todo: no dali e pravilno?
eqns[,1] <- - colSums(sc@m) # separate formula for intercepts
for(i in seq_len(pm)){ # 1:pm
eqns[,i+1] <- - inner( y[(pm+1 -i):(n-i)], sc)
}
res <- eqns[,1]
for(i in 1:g){
res <- c(res, eqns[i,1+seq_len(p[i])])
}
# cat("feqns\n")
res # -res ne mozhe tuk da smenyam znaka, ponezhe smenya i na interceptite!
# 2012-10-31 todo: No zasto da ne go smenya za tyach? !!! check!!!
}
}else if(est_ar_all){ # all shift parameters fixed and all ar parameters not fixed.
function(par){
model@arcoef <- rag_modify(model@arcoef, par)
stdetk <- mix_ek(model, y, index, scale=TRUE) # standardised resid.
# cat("feqns 2\n")
# todo: need methods for unitary "-"
# currently the sign is changed just before returning the value.
# sc <- - (tau * (Fscore[[1]] %of% stdetk)) / sigma
sc <- (tau * (Fscore %of% stdetk)) / sigma
# cat("feqns 3\n")
for(i in seq_len(pm)){ # 1:pm
eqns[,i] <- - inner( y[(pm+1 -i):(n-i)], sc)
}
# cat("feqns 4\n")
res <- numeric(0)
for(i in 1:g){
res <- c(res, eqns[i,seq_len(p[i])]) # [i,1:p[i]]
}
# cat("feqns\n")
res
}
}else{ # estimate subset of the shift-ar parameters
function(par){
shift_ar_params[dontfix_shift_ar] <- par
model@shift <- shift_ar_params[1:g]
model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])
stdetk <- mix_ek(model, y, index, scale=TRUE) # standardised resid.
sc <- (tau * (Fscore %of% stdetk)) / sigma
# TODO: tova tyrabva da se vidi vnimatelno!
# 2012-10-31 smenyam znaka; todo: no dali e pravilno?
eqns[,1] <- - colSums(sc@m) # separate formula for intercepts
for(i in seq_len(pm)){ # 1:pm
eqns[,i+1] <- - inner( y[(pm+1 -i):(n-i)], sc)
}
res <- eqns[,1]
for(i in 1:g){
res <- c(res, eqns[i,1+seq_len(p[i])])
}
res[dontfix_shift_ar] # return the required subset of the equations.
}
}
# 2012-11-02 smen yam s dolnoto zaradi promenite
# if(est_shift)
# param <- c( model@shift, ragged2vec(model@arcoef) )
# else
# param <- ragged2vec(model@arcoef)
if(white_noise){
tmp_v <- ragged2vec(model@arcoef)
param <- c( model@shift[-g], tmp_v[1:(length(tmp_v) - model@order[g])] )
}else if(est_shift)
param <- c( model@shift, ragged2vec(model@arcoef) )
else if(est_ar_all)
param <- ragged2vec(model@arcoef)
else{ # subset of shift-ar parameters
param <- c( model@shift, ragged2vec(model@arcoef) )
param <- param[dontfix_shift_ar]
}
# print(param) # cat(feqns(param), "\n", sep = " ")
if(verbose){
cat("Estimating")
cat(if(lessverbose) "." else " AR parameters...")
}
oldparam <- param
param <- switch(method,
# BBsolve = {print("Ouch!"); BBsolve(par=param, fn = feqns, ...) },
BBsolve = {BBsolve(par=param, fn = feqns, quiet = TRUE,
control = list(maxit=armaxit, trace = trace), ...) },
# requires monotonicity (set M=1 in "control"?)
BBoptim = BBoptim(par=param, fn = lcond,
control=list(maximize=TRUE, trace = trace), ...),
BBoptimgr = BBoptim(par=param, fn = lcond, gr = feqns,
control=list(maximize=TRUE, trace = trace), ...),
BBspg = spg(par=param, fn = lcond, gr = feqns,
control=list(maximize=TRUE, trace = trace), ...),
stop("no method specified")
)
if(verbose)
cat(if(lessverbose) param$convergence else c("convergence =", param$convergence, "fn.reduction =", param$fn.reduction, "\n"))
if(param$convergence > 0 && !white_noise){ # TODO: tova e krapka!
oldlik <- cond_loglik(model, y)
tmpmodel <- model
if(est_shift){
tmpmodel@shift <- param$par[1:g]
tmpmodel@arcoef <- rag_modify(model@arcoef, param$par[-(1:g)])
}else if(est_ar_all){
tmpmodel@arcoef <- rag_modify(model@arcoef, param$par)
}else{
shift_ar_params[dontfix_shift_ar] <- param$par
tmpmodel@shift <- shift_ar_params[1:g]
tmpmodel@arcoef <- rag_modify(tmpmodel@arcoef, shift_ar_params[-(1:g)])
}
tmplik <- cond_loglik(tmpmodel, y)
# todo: more work needed, maybe enforce nondecrease?
if(tmplik < oldlik){ # not only lack of convergence but also worse likelihood,
# try again enforcing monotonicity (M=1)
if(verbose)
cat("tmplik:", tmplik, "is smaller than the old lik,", oldlik, "\n")
tmpparam <- BBsolve(par=oldparam, fn = feqns, quiet = TRUE,
control = list(maxit=max(armaxit,100),
trace = trace, M = 1), ...)
if(verbose)
cat("\n\tconvergence =", tmpparam$convergence, "fn.reduction =",
tmpparam$fn.reduction, "\n")
param <- tmpparam
}
}
# if(est_shift){
# model@shift <- param$par[1:g]
# model@arcoef <- rag_modify(model@arcoef, param$par[-(1:g)])
# }else{
# model@arcoef <- rag_modify(model@arcoef, param$par)
# }
if(white_noise){ # 2021-03-31 new
shift_ar_params[unconstr_indx] <- param$par
model@shift <- shift_ar_params[1:g]
model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])
}else if(est_shift){
model@shift <- param$par[1:g]
model@arcoef <- rag_modify(model@arcoef, param$par[-(1:g)])
}else if(est_ar_all){
model@arcoef <- rag_modify(model@arcoef, param$par)
}else{
shift_ar_params[dontfix_shift_ar] <- param$par
model@shift <- shift_ar_params[1:g]
model@arcoef <- rag_modify(model@arcoef, shift_ar_params[-(1:g)])
}
curlik <- cond_loglik(model, y)
tmpmodel2 <- model
if(white_noise){
## update mixing weights
ui <- rbind(rep(-1, g-1), diag(g-1))
ci <- c(-1, numeric(g-1))
f_wn_pi <- function(par) {
tmpmodel2@prob[1:(g-1)] <- par
tmpmodel2@prob[g] <- 1 - sum(par)
tmpmodel2@shift[g] <- - sum(par * tmpmodel2@shift[1:(g-1)]) / (1 - sum(par))
constr_ar_param <- - par %*% tmpmodel2@arcoef[1:(g-1), ] / (1 - sum(par))
tmpmodel2@arcoef[g, ] <- as.numeric(constr_ar_param)
- cond_loglik(tmpmodel2, y) # 'minus' since constrOptim minimises
}
new_probs <- constrOptim(tmpmodel2@prob[1:(g-1)], f_wn_pi, grad = NULL, ui = ui, ci = ci)
tmpmodel2@prob[1:(g-1)] <- new_probs$par
tmpmodel2@prob[g] <- 1 - sum(new_probs$par)
tmpmodel2@shift[g] <- - sum(new_probs$par * tmpmodel2@shift[1:(g-1)]) / (1 - sum(new_probs$par))
tmpmodel2@arcoef[g, ] <-
as.numeric( - new_probs$par %*% tmpmodel2@arcoef[1:(g-1), ] / (1 - sum(new_probs$par)) )
}else{
tmpmodel2@prob <- prob
}
problik <- cond_loglik(tmpmodel2, y)
if(!white_noise && problik < curlik){
if(verbose) cat("\tBoshwarning: using the estimated prob reduces the likelihood.\n")
}
# 2021-04-04 was: model@prob <- prob # update the mixture probabilitties
model <- tmpmodel2
if(verbose) cat(if(lessverbose) "." else " scale parameters...")
etk <- mix_ek(model, y, index)
if(any(dontfix_scale)){ # update the scale parameters
sigmahat <- em_est_sigma(tau, etk, Fscore, model@scale,
dontfix = dontfix_scale,
compwise = comp_sigma)
# was: tauetk2sigmahat(tau,etk)
model@scale <- sigmahat
}
if(estdist_flag){
if(verbose)
cat(if(lessverbose) "." else " noise parameter(s)...")
nu <- noise_params(model)
# print("nu is: ")
# print(nu)
parscore <- lapply( dist, function(x) x$Parscore )
wrklogpdf <- lapply( dist, function(x) x$logpdf )
# cat("Entering em_est_dist\n")
nu <- em_est_dist(tau, etk, parscore, sigma, nu, wrklogpdf)
model <- set_noise_params(model, nu)
}
if(verbose)
cat("\n")
model
}
# 2011-11-03 name was em_sigma
# 2012-10-23 prerabotvam osnovno
em_est_sigma <- function(tau, etk, Fscore, sigma, dontfix = rep(TRUE, length(sigma)),
compwise = FALSE){
# cat("input: ", sigma, "\n")
# cat("colSums tau: ", colSums(tau@m), "\n")
tauetk <- tau * etk # no need to compute it each time inside f
mtauetk <- tauetk@m
metk <- etk@m
tausums <- colSums(tau@m)
# todo: check if Fscore is guaranteed to be a list of g components
getFscore <- function(k){ # min() here is in case Fscore is of length one
if(is.function(Fscore)) Fscore else Fscore[[min(length(Fscore),k)]]
}
# todo: transformations for 'par' to allow, e.g., common parameters.
# no za da mozhe da se izpolzva tova tryabva i obst mechanizam.
if(all(dontfix)){
f <- function(par){ # multi-variate
wrk <- tauetk * (Fscore %of% (etk/par))
wrk <- - colSums(wrk@m) / tausums
wrk - par # equation is wrk - par = 0
}
fk <- function(par){ # component-wise # par here is scalar
wrk <- mtauetk[,k_cur] * sapply(metk[,k_cur]/par, getFscore(k_cur))
wrk <- - sum(wrk) / tausums[k_cur]
wrk - par # equation is wrk - par = 0
}
}else{
si <- sigma
f <- function(par){ # multi-variate
si[dontfix] <- par
wrk <- tauetk * (Fscore %of% (etk/si))
wrk <- - colSums(wrk@m) / tausums
wrk[dontfix] - par # equation is wrk - par = 0
}
}
## (2021-04-03) TODO: temporary patch !!! Find more permanent solution!
sum_f_sq <- function(par){
sum(f(par)^2)
}
# the equations for sigma can be solved independently
# for each component. need to change `f()' above
# though... a lazy solution would be (wrk-par)[k]
if(compwise){
newsigma <- sigma
for(k_cur in seq_along(sigma)){
sigmak_cur <- BBsolve(par=sigma[k_cur], fn = fk, quiet = TRUE,
control=list(trace = FALSE))
newsigma[k_cur] <- sigmak_cur$par
}
# print(cbind( oldsigma=sigma, newsigma=newsigma, diffsigma=newsigma-sigma))
sigma <- newsigma
}else if(all(dontfix)){
## (2021-04-03) TODO: temporary patch !!! Find more permanent solution!
cat("sigma = ", sigma, "\n")
if(any(sigma < 0.01))
print("inside!\n")
sigma[sigma < 0.01] <- 0.01
## (2021-04-03) TODO: temporary patch !!! Find more permanent solution!
# newsigma <- BBsolve(par=sigma, fn = f, quiet = TRUE,
# control=list(trace = FALSE))
newsigma <- try( BBsolve(par=sigma, fn = f, quiet = TRUE,
control=list(trace = FALSE)) )
if(inherits(newsigma, "try-error") || any( newsigma$par < 0.01 ))
newsigma <- optim(par=sigma, fn = sum_f_sq, method = "L-BFGS-B", lower = 0.01)
newsigma <- newsigma$par
# print(cbind( oldsigma=sigma, newsigma=newsigma, diffsigma=newsigma-sigma))
sigma <- newsigma
}else{
## (2021-04-03) TODO: temporary patch !!! Find more permanent solution!
if(any(sigma[dontfix] < 0.01))
sigma[dontfix][sigma[dontfix] < 0.01] <- 0.01
newsigma <- BBsolve(par=sigma[dontfix], fn = f, quiet = TRUE,
control=list(trace = FALSE))
newsigma <- newsigma$par
# print(cbind( oldsigma=sigma, newsigma=newsigma, diffsigma=newsigma-sigma))
sigma[dontfix] <- newsigma
}
# todo: krapka, uravnenieto se reshava v sigma i -sigma
if(any(sigma<0)){ # kogato Fscore e nechetna function.
sigma <- abs(sigma)
## 2010-06-12 was: print("Changed sign of some sigma[k] to make them positive")
}
# cat("newsigma: ", sigma, "\n")
# if(any(is.nan(wrk))) wrk[is.nan(wrk)] <- 0.0000001
sigma
}
em_est_dist <- function(tau, etk, parscore, sigma, nu, logpdf){
# cat("input: ", nu, "\n")
# cat("colSums tau: ", colSums(tau@m), "\n")
lp <- logpdf
xeval <- etk/sigma
f <- function(nu){
wrk <- tau * (parscore %of% xeval)
wrk <- colSums(wrk@m)
wrk # equation is wrk = 0
}
# origlpnu <- get("nu", environment(lp)) - this was when lp =logpdf[[1]]
parenv <- parent.env(environment(lp[[1]]))
if(exists("creator_fun", parent.env(parenv), inherits=FALSE)) # todo: krapka!
parenv <- parent.env(parenv)
else if(exists("creator_fun", parent.env(parent.env(parenv)), inherits=FALSE))
parenv <- parent.env(parent.env(parenv)) # TODO: krapka! 2012-10-30 !!!
# tova stana s vavezhdaneto na ed_skeleton
origlpnu <- get("nu", parenv) # cat("origlpnu is: ", origlpnu, "\n")
if(exists("set_non_fixed", parenv, inherits=FALSE)){ # todo: krapka 2012-10-26
assign("nu", nu, parenv) # todo: krapka, 'nu' currently comes as the
nu <- evalq(get_non_fixed(), parenv) # whole nu, not as non-fixed only.
f0 <- function(nu){
# evalq does not do here since we need nu to be evaluated before the call.
# evalq(set_non_fixed(nu), parenv) # assign("nu", nu, parenv)
do.call("set_non_fixed", list(nu), envir = parenv)
wrk <- tau * (lp %of% xeval)
if(any(is.nan(wrk@m))) # todo: krapka, no: NaN appear from 0*(-Inf)
wrk@m[is.nan(wrk@m)] <- 0 # i.e. elem of tau = 0 and corresp value = -Inf
# todo: tova sigurno ne e edinstvenata vazmozhnost.
wrk <- colSums(wrk@m)
sum(wrk)
}
}else
f0 <- function(nu){
assign("nu", nu, parenv) # assign("nu", nu, environment(lp))
wrk <- tau * (lp %of% xeval)
if(any(is.nan(wrk@m))) # todo: krapka, no: NaN appear from 0*(-Inf)
wrk@m[is.nan(wrk@m)] <- 0 # i.e. elem of tau = 0 and corresp value = -Inf
# todo: tova sigurno ne e edinstvenata vazmozhnost.
wrk <- colSums(wrk@m)
sum(wrk)
}
# the equations for nu can be solved independently for each component.
# need to change `f()' above though... a lazy solution would be wrk[k]
# todo: nu must be > 2 !!!
# newnu <- BBsolve(par=nu, fn = f)
newnu <- BBoptim(par=nu, fn = f0, lower = 2.2,
control = list(maximize=TRUE, trace = FALSE), quiet = TRUE)
assign("nu", origlpnu, parenv) # restore the original parameter(s) of lp
newnu <- newnu$par
# print(cbind( oldnu=nu, newnu=newnu, diffnu=newnu-nu))
if(exists("set_non_fixed", parenv, inherits=FALSE)){ # todo: krapka 2012-10-26
# vzh i v nachaloto na tazi funktsiya
do.call("set_non_fixed", list(newnu), envir = parenv)
nu <- evalq(get_nu(), parenv)
}else
nu <- newnu
nu
}
mixAR_cond_probs <- function(model, y, indx = NULL){
if(is.null(indx))
indx <- (max(model@order) + 1) : length(y)
stdetk <- mix_ek(model, y, indx, scale=TRUE) ## em_tau needs standardised resid.
# em_tau_safe below takes care of this but is currently
# commented out. Probably need better solution.
# if(any(!is.finite(stdetk))) # todo: krapka! appropriate if matching
# stdetk@m[!is.finite(stdetk)] <- Inf # p_k, scale_k are both (close to) zero.
fpdf <- noise_dist(model, "pdf") # note: this may change if noise dist. has parameters
# so, in general needs to be computed every time
# em_tau_safe(stdetk, model@prob, model@scale, fpdf)
em_tau(stdetk, model@prob, model@scale, fpdf) # todo: dali za slozha i argument indx?
}
mixARgenemFixedPoint <- function(y, model # crit # opts
# ne, zasega chrez '...', fix = "" # 2012-11-02 new arg.
# 2012-11-01 comment out, they are for '...'
# , est_shift = TRUE
# , comp_sigma = FALSE
, crit = 1e-14 # maybe set to a larger value?
, maxniter = 200
, minniter = 10
, verbose = FALSE
, ... ){
verbose <- verbose && interactive()
y <- as.numeric(y)
oldmodel <- model
p <- oldmodel@order
pm <- max(oldmodel@order)
n <- length(y)
indx <- (pm+1):n
nmp <- n - pm
oldvallogf <- cond_loglik(oldmodel, y)
newvallogf <- NA
relchange <- crit + 1
all_relchange <- numeric(0) # for testing
all_vallogf <- oldvallogf
if(verbose)
cat("niter: ", 0, "\t", "vallogf: ", oldvallogf, "\n")
niter <- 0
while( niter <= minniter || niter <= maxniter && relchange > crit ){
niter <- niter + 1
# 2012-10-23 zamenyam nyakolko komandi s edin call na mixAR_cond_probs
newtau <- mixAR_cond_probs(oldmodel, y, indx)
# if(any(colSums(newtau@m) < 2)){ # todo: arbitrary threshold here!
# print(newtau)
# stop("One or more components are too improbable")
# }
newmodel <- mixgenMstep(y, newtau, oldmodel, indx, ...)
oldvallogf <- newvallogf
newvallogf <- cond_loglik(newmodel, y)
relchange <- abs((oldvallogf - newvallogf)/oldvallogf)
all_vallogf <- c(all_vallogf, newvallogf)
all_relchange <- c(all_relchange, relchange)
if(niter %% 10 == 1){ # todo: replace 10 with an argument.
if(verbose){
cat("niter: ", niter, "\t", "vallogf: ", newvallogf, "\n")
show(newmodel)
}
}
oldmodel <- newmodel
}
list(model=newmodel, vallogf=newvallogf, niter = niter
, all_vallogf = all_vallogf
, all_relchange = all_relchange)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.