## Do not edit this file manually.
## It has been automatically generated from mixAR.org.
ui <- function(x, t, i){
## 2018-11-13 this is extremely inefficient, replaced by Georgi
##
## out <- numeric(length(i))
## for(k in 1:length(i)){
## out[k] <- if (i[k] == 0)
## 1
## else
## x[t - i[k]]
## }
## out
if(any(t <= i))
stop("all elements of 't - i' must be positive")
## this is a one-liner but may fail if t > length(x):
## res <- ifelse(i == 0, 1, x[t - i])
if(any(i == 0)){
res <- rep.int(1, length(i))
res[i != 0] <- x[t - i[i != 0]]
res
}else
x[t - i]
}
## Version with intercept
mixSARfit <- function(y, model, est_shift = FALSE, tol = 10^-14){
verbose <- interactive() && options("verbose")$verbose
one_or_zero <- if(est_shift) 0 else 1
n <- length(y)
spk <- model@arcoef@ps
sp <- max(spk)
s <- model@arcoef@s
g <- length(model@prob)
pk <- model@order
p <- max(pk)
index <- ((sp * s) + 1):n
oldlik <- cond_loglikS(model, y, index)
if(verbose)
cat("Initial vallogf: ", oldlik, "\n")
nc <- sp + p + one_or_zero
nck_all <- pk + spk + one_or_zero
sigma0 <- model@scale
pi0 <- model@prob
diff <- 1
count <- 0
while(diff > tol){
count <- count+1
## rows contain same time, columns contain same component
e <- matrix(0, nrow = n - sp * s, ncol = g)
for(k in 1:g){
for(t in 1:(n - sp * s)){
e[t, k] <- y[t + sp * s] -
model@shift[k] -
sum(model@arcoef@a[[k]] * y[(t + sp * s - 1):(t + sp * s - pk[k])]) -
sum(model@arcoef@as[[k]] * ui(y, t + sp * s, s * 1:spk[k]))
}
}
dens.e <- dnorm(t( t(e) / sigma0 ))
tau <- t(pi0 / sigma0 * t(dens.e)) / rowSums(t(pi0 / sigma0 * t(dens.e)))
if(g == 2)
pi <- c(sum(tau[ , 1]) / (n - sp * s), 1 - sum(tau[ , 1]) / (n - sp * s))
else if(g > 2){
pi <- colSums(tau[ , 1:(g - 1)]) / (n - sp * s)
# @Davide: there was a space between '<' and '-'
# pi < -c(pi, 1 - sum(pi))
pi <- c(pi, 1 - sum(pi))
}
sigma <- sqrt(colSums(tau * e ^ 2) / colSums(tau))
phimat <- matrix(0, nrow = g, ncol = nc) # matrix to contain coefficients during computation
for(k in 1:g){
nck <- nck_all[k]
yy <- numeric(nck)
## @Davide: is it guaranteed that pk[k] > 0 here?
## @Georgi: I assumed "MixAR" objects by default have minimum pk[k] = 1
## I can rewrite this if needs to handle pk[k] = 0
##
## @Davide: if your code assume this, raise an error - otherwise the program may
## raise an error for unknown reason or, worse, silently produce wrong
## results.
##
## This is yet another reason not to do everything everytime from scratch, since
## the provided functions take care of such things (even if they don't, they can
## be improved to do so.
onepkk <- seq_len(pk[k]) # == 1:pk[k] when pk[k] > 0
sonespk <- s * seq_len(spk[k])
## Left hand side
for(t in 1:(n - sp * s)){
tsps <- t + sp * s
tsps1mk <- tsps - onepkk # == (tsps - 1):(tsps - pk[k]) == tsps - 1:pk[k])
instr <- y[tsps1mk]
if(one_or_zero == 1)
instr <- c(1, instr)
if(spk[k] > 0) ## i.e. if seasonal coefficients
instr <- c(instr, ui(y, tsps, sonespk))
stopifnot(all(tsps1mk == (tsps - 1):(tsps - pk[k]))) # for testing
yy <- yy + tau[t, k] * y[tsps] * instr
}
## Right hand side
gamma <- matrix(0, nrow = nck, ncol = nck)
phi <- numeric(nc)
for(x in 1:nrow(gamma)){
for(t in 1:(n - sp * s)){
tsps <- t + sp * s
tsps1mk <- tsps - onepkk # == (tsps - 1):(tsps - pk[k]) == tsps - 1:pk[k])
instr <- y[tsps1mk]
if(one_or_zero == 1)
instr <- c(1, instr)
if(spk[k] > 0)
instr <- c(instr, ui(y, tsps, sonespk))
if(x == one_or_zero) # i.e. x == 1 and there is intercept
gamma[x, ] <- gamma[x, ] + tau[t, k] * instr
else if(x <= pk[k] + one_or_zero)
gamma[x, ] <- gamma[x, ] +
tau[t, k] * y[tsps - x + one_or_zero] * instr
else if(spk[k] > 0)
gamma[x, ] <- gamma[x, ] +
tau[t, k] * (
if(one_or_zero == 1)
y[tsps - (x - pk[k] - one_or_zero) * s]
else
ui(y, tsps, (x - pk[k] - one_or_zero) * s)
) * instr
}
}
phi[1:nck] <- pseudoInverse(gamma) %*% yy
phimat[k, ] <- phi
}
## Formatting results into "MixAR" object
ar1 <- ar2 <- vector(mode = "list", g) # list()
for(k in 1:g){
ar1[[k]] <- phimat[k, (1 + one_or_zero):(pk[k] + one_or_zero)]
ar2[[k]] <- phimat[k, (pk[k] + 1 + one_or_zero):(pk[k] + spk[k] + one_or_zero)]
}
ar <- new("raggedCoefS", a = ar1, as = ar2, s = s)
newmod <- if(one_or_zero == 1)
new("MixARGaussian", prob = pi, scale = sigma, arcoef = ar,
shift = phimat[, 1])
else
new("MixARGaussian", prob = pi, scale = sigma, arcoef = ar)
newlik <- cond_loglikS(newmod, y, index)
if(count %% 25 == 0 && verbose)
cat("niter: ", count, "\tvallogf: ", newlik, "\n")
if(is.nan(newlik) && verbose)
cat("!!!! Log-likelihood is NaN, maybe due to singularity.\n")
## Update stopping criterion and model
diff <- abs((oldlik - newlik) / oldlik)
oldlik <- newlik
model <- newmod
sigma0 <- model@scale
pi0 <- model@prob
if(count == 200)
break ## max number of iterations
}
if(verbose)
cat("Final niter: ", count, "\tvallogf: ", newlik, "\n")
list("model" = newmod, "vallogf" = newlik)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.