##' A method to run the mcmc in R
##'
##' @title runmcmc
##' @param result An object from the initialize mcmc function
##' @param ob An object from the dataformatting function
##' @param MCMC the number of iterations (Same as the initialize mcmc
##' function used)
##' @param h A control parameter (for avoiding very small values of q
##' ~Dirich, maybe never needed)
##' @param FULL Choose 1 for full Bayesian model (semi-supervised),
##' but Choose 0 for separated estimation (supervised) of relative
##' type frequencies
##' @return A list of objects to be passed onto the plot function
##' @export
##' @author Jukka Ranta
runmcmc <- function(result,
ob,
MCMC,
h = 0,
FULL = 0){
for(mc in 2:MCMC){
## FULL CONDITIONAL DISTRIBUTIONS FOR
## THE RELATIVE ALLELE TYPE FREQUENCIES IN EACH SOURCE POPULATION:
## (qASP, GLN, GLT, GLY, PGM, TKT, UNC: these are defined for each loci)
for(i in 1:ob$inits$ns){ # ns = number of source populations
for(j in 1:ob$inits$nat[1]){
result$ASPpar0[i, j] <- rgamma(1,
ob$data$sourcesASP[i, j] + FULL * sum(result$S[, i] * ob$data$IASP[, j] > 0) + 1/ob$inits$nat[1],
1)
}
for(j in 1:ob$inits$nat[1]){
result$qASP[mc, i, j] <- (result$ASPpar0[i, j] + h) / (sum(result$ASPpar0[i, ]) + h * ob$inits$nat[1])
}
for(j in 1:ob$inits$nat[2]){
result$GLNpar0[i, j] <- rgamma(1,
ob$data$sourcesGLN[i, j] + FULL * sum(result$S[, i] * ob$data$IGLN[, j] > 0) + 1/ob$inits$nat[2],
1)
}
for(j in 1:ob$inits$nat[2]){
result$qGLN[mc, i, j] <- (result$GLNpar0[i, j] + h) / (sum(result$GLNpar0[i, ]) + h * ob$inits$nat[2])
}
for(j in 1:ob$inits$nat[3]){
result$GLTpar0[i, j] <- rgamma(1,
ob$data$sourcesGLT[i, j] + FULL * sum(result$S[, i] * ob$data$IGLT[, j] > 0) + 1/ob$inits$nat[3],
1)
}
for(j in 1:ob$inits$nat[3]){
result$qGLT[mc, i, j] <- (result$GLTpar0[i, j] + h) / (sum(result$GLTpar0[i, ]) + h * ob$inits$nat[3])
}
for(j in 1:ob$inits$nat[4]){
result$GLYpar0[i, j] <- rgamma(1,
ob$data$sourcesGLY[i, j] + FULL * sum(result$S[, i] * ob$data$IGLY[, j] > 0) + 1/ob$inits$nat[4],
1)
}
for(j in 1:ob$inits$nat[4]){
result$qGLY[mc, i, j] <- (result$GLYpar0[i, j] + h) / (sum(result$GLYpar0[i, ]) + h * ob$inits$nat[4])
}
for(j in 1:ob$inits$nat[5]){
result$PGMpar0[i, j] <- rgamma(1,
ob$data$sourcesPGM[i, j] + FULL * sum(result$S[, i] * ob$data$IPGM[, j] > 0) + 1/ob$inits$nat[5],
1)
}
for(j in 1:ob$inits$nat[5]){
result$qPGM[mc, i, j] <- (result$PGMpar0[i, j] + h) / (sum(result$PGMpar0[i, ]) + h * ob$inits$nat[5])
}
for(j in 1:ob$inits$nat[6]){
result$TKTpar0[i, j] <- rgamma(1,
ob$data$sourcesTKT[i, j] + FULL * sum(result$S[, i] * ob$data$ITKT[, j] > 0) + 1/ob$inits$nat[6],
1)
}
for(j in 1:ob$inits$nat[6]){
result$qTKT[mc, i, j] <- (result$TKTpar0[i, j] + h) / (sum(result$TKTpar0[i, ]) + h * ob$inits$nat[6])
}
for(j in 1:ob$inits$nat[7]){
result$UNCpar0[i, j] <- rgamma(1,
ob$data$sourcesUNC[i, j] + FULL * sum(result$S[, i] * ob$data$IUNC[, j] > 0) + 1/ob$inits$nat[7],
1)
}
for(j in 1:ob$inits$nat[7]){
result$qUNC[mc, i, j] <- (result$UNCpar0[i, j] + h) / (sum(result$UNCpar0[i, ]) + h * ob$inits$nat[7])
}
} ## end of ns
## FULL CONDITIONAL DISTRIBUTION FOR EACH Z:
## (Z = group indicators for each isolate)
for(i in 1:result$Nisolates){
for(k in 1:ob$inits$ns){
result$phii0[i, k] <- exp( log(result$phi[mc-1, k]) +
log(result$qASP[mc, k, ob$data$ind[i, 1]]) +
log(result$qGLN[mc, k, ob$data$ind[i, 2]]) +
log(result$qGLT[mc, k, ob$data$ind[i, 3]]) +
log(result$qGLY[mc, k, ob$data$ind[i, 4]]) +
log(result$qPGM[mc, k, ob$data$ind[i, 5]]) +
log(result$qTKT[mc, k, ob$data$ind[i, 6]]) +
log(result$qUNC[mc, k, ob$data$ind[i, 7]])
)
}
for(k in 1:ob$inits$ns){
result$phii[i, k] <- exp(log(result$phii0[i, k])-log(sum(result$phii0[i, ])) ) ## normalize
}
nn <- 1:ob$inits$ns
result$Z[mc, i] <- min(nn[runif(1) < cumsum(result$phii[i, ])]) ## "rcat(1, result$phii[i, ])"
for(k in 1:ob$inits$ns){
result$S[i, k] <- result$Z[mc, i]==k ## membership of ith isolate to group k
}
} ## end Nisolates
## FULL CONDITIONAL DISTRIBUTION FOR phi:
## (phi = proportions of source populations)
for(k in 1:ob$inits$ns){
result$phi0[k] <- rgamma(1,
sum(result$S[, k]) + 1/ob$inits$ns,
1)
}
for(k in 1:ob$inits$ns){
result$phi[mc, k] <- result$phi0[k]/sum(result$phi0[])
}
} ## end of MCMC
result <- list(var_a = result, var_b = ob)
class(result) <- "kilde_rmcmc"
return(result)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.