## Do not edit this file manually.
## It has been automatically generated from mixAR.org.
## Function to sample latentZ and probabilities
sampZpi <- function(y, pk, prob, mu, AR, sigma, nsim, d){
g <- length(pk)
if(missing(d)) d <- rep(1, g)
else if (length(d) == 1) d <- rep(d, g)
if(length(d) != g){
stop("One hyperparameter for each component is needed.\n
length(d) must be 1 or g")
}
pid <- matrix(ncol = g, nrow = nsim)
n <- length(y)
p <- max(pk)
for(j in 1:nsim){
z <- pz <- matrix(0, n-p, g)
for(t in 1:(n-p)){
for(k in 1:g){
pz[t, k] <- prob[k] *
dnorm(y[t+p], mu[k] +
(y[(t + p - 1) : (t+p-pk[k])]-mu[k])%*%AR[[k]],
sigma[k])
}
}
fu <- function(x) sample.int(g, 1L, replace = TRUE, prob = x/sum(x))
obs <- apply(pz, 1, fu)
ind <- matrix(c(1:(n-p), obs), ncol=2)
z[ind] <- 1
## Update proportions
nk <- colSums(z)
pi <- rdirichlet(1, d + nk)
pid[j, ] <- pi
}
nam <- nai <- nank <- 0
for(i in 1:g){
nam[i] <- paste("prob ", i)
nai[i] <- paste("latent Z ", i)
nank[i] <- paste("nk ", i)
}
colnames(pid) <- nam
colnames(z) <- nai
names(nk) <- nank
list(mix_weights = pid, latentZ = z, nk = nk)
}
## Function to sample mu and shift
sampMuShift <- function(y, pk, prec, nk, shift, z, AR, nsim){
Ry <- max(y) - min(y)
zeta <- min(y) + Ry / 2
ka <- 1 / Ry
g <- length(pk)
mu <- numeric(g)
bk <- sapply(AR, function(x) 1 - sum(x))
mu[which(shift != 0)] <- shift[which(shift != 0)] / (bk[which(shift != 0)])
mud <- intd <- matrix(ncol = g, nrow = nsim)
nam <- nai <- 0
for(i in 1:g){
nam[i] <- paste("mu ", i)
nai[i] <- paste("shift ", i)
}
colnames(mud) <- nam
colnames(intd) <- nai
int <- numeric(g)
n <- length(y)
p <- max(pk)
dens_mu <- rep(1, nsim)
## The following 'e' is because "e[t,k]" here is assumed to be
## y[t] - sum(y[t-1 : t-pk] * AR[[k]])
## i.e. e is equal to err(....)[t,k] + shift[k] (for elements different from 0)
e <- err(AR, mu, y, z, pk)
fu <- function(x){k <- which(x != 0); x[k] <- x[k] +shift[k]; x }
e <- t(apply(e, 1, fu))
for(j in 1:nsim){
for(k in 1:g){
index <- which(e[,k] !=0)
mu[k] <- rnorm(1, (prec[k] * nk[k] * mean(e[index, k]) *
bk[k] + ka * zeta) /
(prec[k] * nk[k] * bk[k]^2 + ka),
1 / sqrt((prec[k] * nk[k] * bk[k]^2 + ka)))
int[k] <- mu[k] * (bk[k])
}
intd[j, ] <- int
mud[j, ] <- mu
}
list(shift = intd, mu = mud)
}
## Function to sample from sigma/tau
sampSigmaTau <- function(y, pk, prec, nk, AR, mu, z, a, c, nsim){
Ry <- max(y) - min(y)
b <- 100 * a / (c * Ry^2)
g <- length(pk)
n <- length(y)
p <- max(pk)
lambda <- numeric(nsim)
sigmad <- precd <- matrix(ncol = g, nrow = nsim)
for(j in 1:nsim){
lambda[j] <- rgamma(1, a + g * c, b + sum(prec))
for(k in 1:g){
prec[k] <- rgamma(1, c + nk[k] / 2,
lambda + 0.5 * sum(err(AR, mu, y, z, pk)[ ,k]^2))
}
sigma <- 1 / sqrt(prec)
sigmad[j, ] <- sigma
precd[j, ] <- prec
}
nam <- nai <- 0
for(i in 1:g){
nam[i] <- paste("sigma ", i)
nai[i] <- paste("precision ", i)
}
colnames(sigmad) <- nam
colnames(precd) <- nai
list(scale = sigmad, precision = precd, lambda = lambda)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.