hidden.r

使用R语言的马尔科夫链蒙特卡洛模拟（MCMC）源代码程序。

    K <- ncol(X)
    if (is.na(psi.start)){
      Psi <- 0.5 * diag(diag(var(X)))
    }
    else if (is.double(psi.start) &&
             (length(psi.start==1) || length(psi.start==K))){
      Psi <- diag(K) * psi.start
    }
    else {
      cat("psi.start neither NA, double. nor appropriately sized matrix.\n")
      stop("Please respecify and call ", calling.function, " again.\n")
    }
    if (nrow(Psi) != K || ncol(Psi) != K){
      cat("Psi starting value not K by K matrix.\n")
      stop("Please respecify and call ", calling.function, " again.\n")    
    }

    return(Psi)
  }


## form the ind. normal prior for a factor loading matrix
"form.factload.norm.prior" <-
  function(l0, L0, K, factors, X.names){

    ## prior means
    if (is.matrix(l0)){ # matrix input for l0
      if (nrow(l0)==K && ncol(l0)==factors){
        Lambda.prior.mean <- l0
        rownames(Lambda.prior.mean) <- X.names
      }
      else {
        cat("l0 not of correct size for model specification.\n")
        stop("Please respecify and call ", calling.function(), " again.\n")
      }
    }
    else if (is.list(l0)){ # list input for l0
      Lambda.prior.mean <- matrix(0, K, factors)
      rownames(Lambda.prior.mean) <- X.names
      l0.names <- names(l0)
      for (i in 1:length(l0.names)){
        name.i <- l0.names[i]
        l0.i <- l0[[i]]
        col.index <- l0.i[[1]]
        replace.element <- l0.i[[2]]
        if (is.numeric(replace.element)){
          Lambda.prior.mean[rownames(Lambda.prior.mean)==name.i,
                            col.index] <- replace.element
        }   
      }
    }
    else if (length(l0)==1 && is.numeric(l0)){ # scalar input for l0
      Lambda.prior.mean <- matrix(l0, K, factors)
      rownames(Lambda.prior.mean) <- X.names
    }
    else {
      cat("l0 neither matrix, list, nor scalar.\n")
      stop("Please respecify and call ", calling.function(), " again.\n")
    }
    
    ## prior precisions
    if (is.matrix(L0)){ # matrix input for L0
      if (nrow(L0)==K && ncol(L0)==factors){
        Lambda.prior.prec <- L0
        rownames(Lambda.prior.prec) <- X.names
      }
      else {
        cat("L0 not of correct size for model specification.\n")
        stop("Please respecify and call ", calling.function(), " again.\n")
      }
    }
    else if (is.list(L0)){ # list input for L0
      Lambda.prior.prec <- matrix(0, K, factors)  
      rownames(Lambda.prior.prec) <- X.names
      L0.names <- names(L0)
      for (i in 1:length(L0.names)){
        name.i <- L0.names[i]
        L0.i <- L0[[i]]
        col.index <- L0.i[[1]]
        replace.element <- L0.i[[2]]
        if (is.numeric(replace.element)){
          Lambda.prior.prec[rownames(Lambda.prior.prec)==name.i,
                            col.index] <- replace.element
        }   
      }
    }
    else if (length(L0)==1 && is.numeric(L0)){ # scalar input for L0
      Lambda.prior.prec <- matrix(L0, K, factors)
      rownames(Lambda.prior.prec) <- X.names
    }
    else {
      cat("L0 neither matrix, list, nor scalar.\n")
      stop("Please respecify and call ", calling.function(), " again.\n")
    }
    if (min(Lambda.prior.prec) < 0) {
      cat("L0 contains negative elements.\n")
      stop("Please respecify and call ", calling.function(), " again.\n")
    } 


    return( list(Lambda.prior.mean, Lambda.prior.prec))    
  }

## form ind. inv. gamma prior for a diagonal var. cov. matrix
"form.ig.diagmat.prior" <-
  function(a0, b0, K){

    ## setup prior for diag(Psi)
    if (length(a0)==1 && is.double(a0))
      a0 <- matrix(a0, K, 1)
    else if (length(a0) == K && is.double(a0))
      a0 <- matrix(a0, K, 1)
    else {
      cat("a0 not properly specified.\n")
      stop("Please respecify and call ", calling.function, " again.\n")
    }
    if (length(b0)==1 && is.double(b0))
      b0 <- matrix(b0, K, 1)
    else if (length(b0) == K && is.double(b0))
      b0 <- matrix(b0, K, 1)
    else {
      cat("b0 not properly specified.\n")
      stop("Please respecify and call ", calling.function(), " again.\n")
    }
    
    ## prior for Psi error checking
    if(min(a0) <= 0) {
      cat("IG(a0/2,b0/2) prior parameter a0 less than or equal to zero.\n")
      stop("Please respecify and call ", calling.function, " again.\n")
    }
    if(min(b0) <= 0) {
      cat("IG(a0/2,b0/2) prior parameter b0 less than or equal to zero.\n")
      stop("Please respecify and call ", calling.function(), " again.\n")      
    }  

    return(list(a0, b0) )
  }

# pull together the posterior density sample
"form.mcmc.object" <-
  function(posterior.object, names, title, ...) {
    holder <- matrix(posterior.object$sampledata,
                     posterior.object$samplerow,
                     posterior.object$samplecol,
                     byrow=FALSE)
      
    output <- mcmc(data=holder, start=(posterior.object$burnin+1),
                   end=(posterior.object$burnin+posterior.object$mcmc),
                   thin=posterior.object$thin)
    varnames(output) <- as.list(names)
    attr(output,"title") <- title
    
    attribs <- list(...)
    K <- length(attribs)
    attrib.names <- names(attribs)

    if (K>0){
      for (i in 1:K){
        attr(output, attrib.names[i]) <- attribs[[i]]
      }
    }
    
    return(output)  
  }

# form multivariate Normal prior
"form.mvn.prior" <-
   function(b0, B0, K) {
  
     # prior mean
     if(is.null(dim(b0))) {
       b0 <- b0 * matrix(1,K,1)  
     } 
     if((dim(b0)[1] != K) || (dim(b0)[2] != 1)) {
       cat("Error: N(b0,B0^-1) prior b0 not conformable.\n")
       stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
     }
     
     # prior precision
     if(is.null(dim(B0))) {
       if (length(B0) > K){
         stop("B0 was passed as a vector longer than K.\nB0 must be either a scalar or a matrix.\nPlease respecify and call ", calling.function(), " again.\n", call.=FALSE)
       }
       B0 <- B0 * diag(K)    
     }
     if((dim(B0)[1] != K) || (dim(B0)[2] != K)) {
       cat("Error: N(b0,B0^-1) prior B0 not conformable.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)
     }

     ## check B0 for symmetry
     symproblem <- FALSE
     for (i in 1:K){
       for (j in i:K){
         if (B0[i,j] != B0[j,i]){
           symproblem <- TRUE
         }
       }
     }
     if (symproblem){
       cat("B0 is not symmetric.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)       
     }
     
     return(list(b0,B0))
   }

# parse the passed seeds
# 1] if a scalar is passed, it is used by Mersennse twister
# 2] if a list of length two is passed, a parallel-friendly stream is
#    created using L'Ecuyer
"form.seeds" <-
   function(seed) {
      if(length(seed)==1) {
         if(is.na(seed)) seed <- 12345
         seed <- as.integer(seed)
         if(seed < 0) {
            cat("Error: Mersenne seed negative.\n")
            stop("Please respecify and call ", calling.function(), " again.",
              call.=FALSE)                       
         }
         seeds <- list(0, rep(seed,6), 0)
      }
      if(length(seed)==2) {
         if(!is.list(seed)) {
            cat("Error: List must be passed to use L'Ecuyer.\n")
            stop("Please respecify and call ", calling.function(), " again.",
              call.=FALSE)          
         }
         lec.seed <- seed[[1]]
         lec.substream <- as.integer(seed[[2]])
         if(is.na(lec.seed[1])) lec.seed <- rep(12345, 6)
         if(length(lec.seed) != 6) {
            cat("Error: L'Ecuyer seed not of length six.\n")
            stop("Please respecify and call ", calling.function(), " again.",
              call.=FALSE)          
         }
         if(!all(lec.seed >= 0))  {
             cat("Error: At least one L'Ecuyer seed negative.\n")
            stop("Please respecify and call ", calling.function(), " again.",
              call.=FALSE)          
         }
         if( max(lec.seed[1:3]) >= 4294967087){
           cat("Error: At least one of first three L'Ecuyer seeds\n")
           cat("  greater than or equal to 4294967087\n")
           stop("Please respecify and call ", calling.function(), " again.",
                call.=FALSE)          
         }
         if( all(lec.seed[1:3]) == 0 ){
           cat("Error: first three L'Ecuyer seeds == 0\n")
           stop("Please respecify and call ", calling.function(), " again.",
                call.=FALSE)          
         }
         if( max(lec.seed[4:6]) >= 4294944443){
           cat("Error: At least one of last three L'Ecuyer seeds\n")
           cat("  greater than or equal to 4294944443\n")
           stop("Please respecify and call ", calling.function(), " again.",
                call.=FALSE)          
         }         
         if( all(lec.seed[4:6]) == 0 ){
           cat("Error: last three L'Ecuyer seeds == 0\n")
           stop("Please respecify and call ", calling.function(), " again.",
                call.=FALSE)          
         }
         if(lec.substream < 1) {
            cat("Error: L'Ecuyer substream number not positive.\n")
            stop("Please respecify and call ", calling.function(), " again.",
                 call.=FALSE)               
         }
         seeds <- list(1, lec.seed, lec.substream) 
      }
      if(length(seed)>2) {
            cat("Error: Seed passed as length greater than two.\n")
            stop("Please respecify and call ", calling.function(), " again.",
              call.=FALSE)        
      }
      return(seeds)
   }

# form Wishart prior
"form.wishart.prior" <-
    function(v, S, K) {
    
    # check to see if degrees of freedom produces proper prior
    if(v < K) {
      cat("Error: Wishart(v,S) prior v less than or equal to K.\n")
      stop("Please respecify and call ", calling.function(), " again.\n")
    } 
    
    # form the prior scale matrix
    if(is.null(dim(S))) {
      S <- S * diag(K)
    }
    if((dim(S)[1] != K) | (dim(S)[2] != K)) {
      cat("Error: Wishart(v,S) prior S not comformable [K times K].\n")
      stop("Please respecify and call ", calling.function(), " again.\n")
    }
           
    return(list(v,S))
}

# parse formula and return a list that contains the model response
# matrix as element one, and the model matrix as element two
"parse.formula" <- 
   function(formula, data, intercept=TRUE, justX=FALSE) {

    # extract Y, X, and variable names for model formula and frame
    mt <- terms(formula, data=data)
    if(missing(data)) data <- sys.frame(sys.parent())
    mf <- match.call(expand.dots = FALSE)
    mf$intercept <- mf$justX <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    if (!intercept){
      attributes(mt)$intercept <- 0
    }

    # null model support
    X <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts)
    X <- as.matrix(X)         # X matrix
    xvars <- dimnames(X)[[2]] # X variable names
    xobs  <- dimnames(X)[[1]] # X observation names
    if (justX){
      Y <- NULL
    }
    else {
      Y <- as.matrix(model.response(mf, "numeric")) # Y matrix
    }
    return(list(Y, X, xvars, xobs))
   }

# setup tuning constant for scalar parameter
"scalar.tune" <- function(mcmc.tune){
  if (max(is.na(mcmc.tune))){
    cat("Error: Scalar tuning parameter cannot contain NAs.\n")
    stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)    
  }
  if (length(mcmc.tune) != 1){
    cat("Error: Scalar tuning parameter does not have length = 1.\n")
    stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
  }
  if (mcmc.tune <= 0) {
    cat("Error: Scalar tuning parameter not positive.\n")
    stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
  } 
  return(mcmc.tune)
}
  
# put together starting values for sigma2
"sigma2.start" <-
   function(sigma2.start, formula, data) {
   
     if(is.na(sigma2.start)){ # use MLE
       lm.out <- lm(formula, data=data)
       sigma2.start <- var(residuals(lm.out))
     }   
     else if(sigma2.start <= 0) {
       cat("Error: Starting value for sigma2 negative.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)
     }
     else if (length(sigma2.start) != 1){
       cat("Error: Starting value for sigma2 not a scalar.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)
     }
     else if (!is.numeric(sigma2.start)){
       cat("Error: Starting value for sigma2 neither numeric nor NA.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)
     }
     return(sigma2.start)   
   }




## setup diagonal tuning matrix for vector parameters
"vector.tune" <- function(mcmc.tune, K){
  if (max(is.na(mcmc.tune))){
    cat("Error: Vector tuning parameter cannot contain NAs.\n")
    stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)    
  }
  if (length(mcmc.tune) == 1){
    mcmc.tune <- rep(mcmc.tune, K)
  }
  if (length(mcmc.tune) != K){
    cat("Error: length(vector tuning parameter) != length(theta) or 1.\n")
    stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
  }
  if (sum(mcmc.tune <= 0) != 0) {
    cat("Error: Vector tuning parameter cannot contain negative values.\n")
    stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
  }
  if (length(mcmc.tune)==1){
    return(matrix(mcmc.tune, 1, 1))
  }
  else{
    return(diag(as.double(mcmc.tune)))
  }
}