## Do not edit this file manually.
## It has been automatically generated from mixAR.org.
## Two functions to calculate component-specific error terms
err <- function(AR, mu, y, z, pk){
n <- length(y)
p <- max(pk)
g <- length(pk)
e <- matrix(0, n-p, g)
for (t in 1 : (n - p)){
for(k in 1:g){
e[t, k] <-
(y[t+p] - mu[k] -
(y[(t + p - 1) : (t + p - pk[k])] - mu[k]) %*% AR[[k]]) * z[t,k]
}
}
e
}
err_k <- function(AR, mu, y, z, p, pk){
n <- length(y)
e <- numeric(n - p)
for (t in 1 : (n - p)){
e[t] <- (y[t + p] - mu - (y[(t + p - 1) : (t + p - pk)] - mu) %*% AR) * z[t]
}
e
}
## Main function
bayes_mixAR <- function(y, model, fix_shift = FALSE, a = 0.2, c = 2, tau, nsim, burnin){
d <- rep(1,length(model@prob))
n <- length(y)
prob <- model@prob
sigma <- model@scale
prec <- 1 / sigma ^ 2
g <- length(model@prob)
pk <- model@arcoef@p
p <- max(pk)
shift <-
mu <- model@shift
AR <- model@arcoef@a
bk <- sapply(AR, function(x) 1 - sum(x))
mu[which(shift != 0)] <- shift[which(shift != 0)] / (bk[which(shift != 0)])
if(length(tau) == 1) tau <- rep(tau, g)
if(length(tau) != 1 & length(tau) != g){
stop("tau should have length equal to 1 or g")
}
## Hyperparameters
Ry <- max(y) - min(y)
zeta <- min(y) + Ry / 2
ka <- 1 / Ry
b <- 100 * a / (c * Ry^2)
## Matrices to store output
pid <-
sigmad <-
precd <-
mud <- matrix(nrow = nsim, ncol = g)
lnames <- NULL
ard <- list()
for(k in 1:g){
ard[[k]] <- matrix(0, nrow = nsim, ncol = pk[k])
lnames <- c(lnames, paste("Component_", k, sep = ""))
}
names(ard) <- lnames
intd <- matrix(nrow = nsim, ncol = g)
count <- rep(0, g)
for(j in 1:nsim){
## Draw z_t's
samp1 <- sampZpi(y, pk, prob, mu, AR, sigma, 1, d)
z <- samp1$latentZ
nk <- samp1$nk
## Update proportions
prob <- samp1$mix_weights
pid[j, ] <- prob
## Update mu
if(fix_shift == FALSE){
samp2 <- sampMuShift(y, pk, prec, nk, shift, z, AR, 1)
mu <- samp2$mu
shift <- samp2$shift
}
intd[j, ] <- shift
mud[j, ] <- mu
## Update precision
samp3 <- sampSigmaTau(y, pk, prec, nk, AR, mu, z, a, c, 1)
lambda <- samp3$lambda
prec <- samp3$precision
sigma <- samp3$scale
sigmad[j, ] <- sigma
precd[j, ] <- prec
## Update phi
ratio <- rep(0, g)
phi_k <- list()
for(k in 1:g){
phi_k[[k]] <- rnorm(pk[k], AR[[k]], tau[k])
ratio[k] <-
min(1, exp(- prec[k]/2 *
sum(err_k(phi_k[[k]], mu[k], y, z[ ,k], p, pk[k]) ^ 2 -
err_k(AR[[k]], mu[k], y, z[ ,k], p, pk[k]) ^ 2
)))
}
cc <- rep(0, g) ## counter vector
u <- runif(g) ## to accept or reject a set of new AR values
for(k in 1:g){
if(u[k] > ratio[k]){ ## if rejected
phi_k[[k]] <- AR[[k]] ## phi_k is set to current values
}else {
if(j > burnin) cc[k] <- cc[k] + 1
}
}## does not count acceptance rate over burnin period
A <- new("MixARGaussian", prob = as.numeric(prob), scale = as.numeric(sigma),
arcoef = phi_k, shift = as.numeric(shift))
if(isStable(A)){ ## check for stability
for(k in 1:g){
ard[[k]][j, ] <- phi_k[[k]] ## add row to output
}
AR <- phi_k ## update current values
if(j > burnin) {count <- count + cc}
} else{
for(k in 1:g){ ## if model is not stable, reject
ard[[k]][j, ] <- as.matrix(AR[[k]])
prob <- pid[j, ] <- pid[j-1, ]
intd[j, ] <- shift <- intd[j-1, ]
mud[j, ] <- mu <- mud[j-1, ]
sigmad[j, ] <- sigma <- sigmad[j-1, ]
precd[j, ] <- prec <- precd[j-1, ]
}
}
}
ard <- lapply(ard, function(x) as.matrix(x[-c(1:burnin), ]))
list(mix_weights = pid[-c(1:burnin), ],
scale = sigmad[-c(1:burnin), ],
precision = precd[-c(1:burnin), ],
shift = intd[-c(1:burnin), ],
mu = mud[-c(1:burnin), ],
ARcoeff = ard, ## now returns a matrix for each K, also when pk=1
acc_rate = count / (nsim - burnin),
n_samp = nsim - burnin,
LatentZ = z,
ncomp = g,
fix_shift = fix_shift)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.