R/Choose_pk.R

In mixAR: Mixture Autoregressive Models

## Do not edit this file manually.
## It has been automatically generated from mixAR.org.

## Give probability of increasing/decreasing pk by one given current pk
## Used within main function

bx_dx <- function(method = c("Ratio", "Poisson", "NULL"), par, pk){
    ## method <- match.arg(method)
    bx <- switch(method,
                 "Poisson" =  dpois(pk + 1, par)/dpois(pk, par),
                 "Ratio" = par / (par + pk),
                 "NULL" = 0.5,
                 stop("Unknown method in 'bx_dx'")
                 )
    dx <- 1 - bx
    list(bx = bx, dx = dx)
}

## "method" calls the method to be used in bx_dx
## "par" is parameter for chosen method
## "par" can be omitted if chosen method is "NULL"

Choose_pk <- function(y, model, fix_shift = FALSE, tau, pmax, method,
                      par = NULL, nsim){
    stopifnot(is(model, "MixARGaussian"))
    
    if(pmax < 0) stop("ar order must be at least 0")
    
    p <- max(model@order)
    
    if(pmax < p) stop("Maximum ar order must be at least p")
    
    if (missing(method))
        method <- "NULL"
    g <- length(model@prob)
    cont <- matrix(0, nrow = nsim, ncol = g)

    for(i in 1:nsim){
        ## Set up model
        mod <- if(i == 1)
                   model
               else
                   new("MixARGaussian", prob = pi, scale = sigma, arcoef = AR, shift = int)

        pk    <- mod@arcoef@p
        p     <- max(pk)
        sim   <- bayes_mixAR(y, mod, .2, 2, tau, 2, 1, fix_shift = fix_shift)
        pi    <- sim$mix_weights
        prec  <- sim$precision
        sigma <- sim$scale
        int   <- sim$shift
        mu    <- sim$mu
        AR    <- sim$ARcoeff
        z     <- sim$LatentZ
        

        comp <- sample(1:g, 1)

        if(pk[comp] == 1)  {      ## Automatically increase by 1
            pkstar <- 2
        }else{
            if(pk[comp] == pmax)
                pkstar <- pmax - 1  ## Automatically decrease by 1
            else{
                bx <-  bx_dx(method = method, par, pk[comp])$bx
                u  <- runif(1)
                if(u < bx) {
                    pkstar <- pk[comp] + 1
                }else {pkstar <- pk[comp] -1}
            }
        }

        ratio <- 0

        ## "if" to handle increase and decrease properly
        if(pkstar > pk[comp]) {
            phistar <- c(AR[[comp]], runif(1,-1.5, 1.5))
            ARcheck <- AR ; ARcheck[[comp]] <- phistar
            Acheck  <- new("MixARGaussian", prob = pi, scale = sigma, arcoef = ARcheck)
            if(isStable(Acheck)){
                ## Required step to increase p to p+1 when new pk > p
                ratio <-
                   exp(- prec[comp]/2 *
                   if(pkstar > p) {
                       sum(-err_k(AR[[comp]], mu[comp], y, z[-1, comp], pkstar, pk[comp]) ^ 2 +
                           err_k(phistar, mu[comp], y, z[-1, comp], pkstar, pkstar) ^ 2)
                   } else {
                       sum(-err_k(AR[[comp]], mu[comp], y, z[ ,comp], p, pk[comp]) ^ 2 +
                           err_k(phistar, mu[comp], y, z[ ,comp], p, pkstar) ^ 2)
                   })
                ratio <- ratio * 3 * bx_dx(method = method, par, pkstar)$dx /
                                     bx_dx(method = method, par, pk[comp])$bx
            }
        }

        if(pkstar < pk[comp]){
            phistar <- AR[[comp]][-pk[comp]]
            ARcheck <- AR
            ARcheck[[comp]] <- phistar
            Acheck  <- new("MixARGaussian", prob = pi, scale = sigma, arcoef = ARcheck)
            if(isStable(Acheck)){
                ratio <- exp(- prec[comp]/2 *
                             sum(- err_k(AR[[comp]], mu[comp], y,
                                                z[,comp], p, pk[comp]) ^ 2 +
                                   err_k(phistar, mu[comp], y,
                                                z[,comp], p, pkstar) ^ 2))
                ratio <- ratio * dnorm(0,0,tau[comp]) *
                    bx_dx(method = method, par, pk[comp] - 1)$bx /
                    bx_dx(method = method, par, pk[comp])$dx
            }
        }

        u <- runif(1)
        if(u < ratio){
            AR[[comp]] <- phistar
            p <- max(Acheck@order)
            pk[comp] <- pkstar
        }
        cont[i,] <- pk
    }

    out <- matrix(0, ncol = g + 1, nrow = nrow(unique(cont)))

    for(i in 1:nrow(unique(cont))){
        out[i, 1:g] <- c(unique(cont)[i, ])
        for(j in 1:nrow(cont)){
            if(all(cont[j, ] == unique(cont)[i,]))
                out[i, g + 1] <- out[i, g + 1] + 1
        }
    }
    out[, g+1] <- out[, g + 1] / nsim
    colnames(out) <- c(paste("Component_", 1:g), "Prop")

    ## out$x is a matrix. Column i corresponds to order of component i.
    ## Last column shows the proportion of times such model was chosen.
    list(x = as.data.frame(out), fix_shift = fix_shift, method = method)
}