## 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.