R/samp_functions.R

In GeoBosh/mixAR: Mixture Autoregressive Models

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