.packageName <- "MCMCpack"
# sample from the posterior of Quinn's dynamic ecological inference model
# in R using linked C++ code in Scythe
#
# KQ 10/25/2002

"MCMCdynamicEI" <-
  function(r0, r1, c0, c1, burnin=5000, mcmc=50000,
           thin=1, verbose=FALSE, seed=NA,
           W=0, a0=0.825, b0=0.0105, a1=0.825,
           b1=0.0105, ...){
    
    # Error checking
    if (length(r0) != length(r1)){
      cat("length(r0) != length(r1).\n")
      stop("Please check data and try MCMCdynamicEI() again.\n")
    }

    if (length(r0) != length(c0)){
      cat("length(r0) != length(c0).\n")
      stop("Please check data and try MCMCdynamicEI() again.\n")
    }

    if (length(r0) != length(c1)){
      cat("length(r0) != length(c1).\n")
      stop("Please check data and try MCMCdynamicEI() again.\n")
    }

    if (length(r1) != length(c0)){
      cat("length(r1) != length(c0).\n")
      stop("Please check data and try MCMCdynamicEI() again.\n")
    }

    if (length(r1) != length(c1)){
      cat("length(r1) != length(c1).\n")
      stop("Please check data and try MCMCdynamicEI() again.\n")
    }

    if (length(c0) != length(c1)){
      cat("length(c0) != length(c1).\n")
      stop("Please check data and try MCMCdynamicEI() again.\n")
    }

    if (min((r0+r1) == (c0+c1))==0){
      cat("Rows and columns do not sum to same thing.\n")
      stop("Please check data and try MCMCdynamicEI() again.\n")
    }

    check.mcmc.parameters(burnin, mcmc, thin)
   
    # seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]

    if (a0 <= 0 ){
      cat("Parameter a0 <= 0.\n")
      stop("Please respecify and try MCMCdynamicEI() again.\n")
    }
    
    if (b0 <= 0 ){
      cat("Parameter b0 <= 0.\n")
      stop("Please respecify and try MCMCdynamicEI() again.\n")
    }
    
    if (a1 <= 0 ){
      cat("Parameter a1 <= 0.\n")
      stop("Please respecify and try MCMCdynamicEI() again.\n")
    }
    
    if (b1 <= 0 ){
      cat("Parameter b1 <= 0.\n")
      stop("Please respecify and try MCMCdynamicEI() again.\n")
    }
    
    ntables = length(r0)
    
    if (W==0){ # construct weight matrix for a simple random walk assuming
               # tables are temporally ordered and 1 time unit apart
      W <- matrix(0, ntables, ntables)
      for (i in 2:(ntables)){
        W[i,i-1] <- 1
        W[i-1,i] <- 1
      }
    }

    # setup matrix to hold output from sampling
    sample <- matrix(0, mcmc/thin, ntables*2+2)

    # call C++ code to do the sampling
    C.sample <- .C("dynamicEI",
                   samdata = as.double(sample),
                   samrow = as.integer(nrow(sample)),
                   samcol = as.integer(ncol(sample)),
                   r0 = as.double(r0),
                   r1 = as.double(r1),
                   c0 = as.double(c0),
                   c1 = as.double(c1),
                   ntables = as.integer(ntables),
                   burnin = as.integer(burnin),
                   mcmc = as.integer(mcmc),
                   thin = as.integer(thin),
                   W = as.double(W),
                   a0 = as.double(a0),
                   b0 = as.double(b0),
                   a1 = as.double(a1),
                   b1 = as.double(b1),
                   verbose = as.integer(verbose),
                   lecuyer = as.integer(lecuyer),
                   seedarray = as.integer(seed.array),
                   lecuyerstream = as.integer(lecuyer.stream),
                   PACKAGE="MCMCpack"
                   )
    
    sample <- matrix(C.sample$samdata, C.sample$samrow, C.sample$samcol,
                     byrow=TRUE)
    output <- mcmc(data=sample, start=1, end=mcmc, thin=thin)
    p0names <- paste("p0table", 1:ntables, sep="")
    p1names <- paste("p1table", 1:ntables, sep="")
    varnames(output) <- c(p0names, p1names, "sigma^2_0", "sigma^2_1")
    
    attr(output, "title") <- "MCMCpack Quinn's Dynamic EI Model Posterior Density Sample" 
    
    
    return(output)
    
  }
##########################################################################
## sample from the posterior distribution of a factor analysis model
## model in R using linked C++ code in Scythe.
##
## The model is:
##
## x_i = \Lambda \phi_i + \epsilon_i,   \epsilon_i \sim N(0, \Psi)
##
## where \Psi is diagonal and the priors are:
##
## \lambda_{ij} \sim N(l_{ij}, L^{-1}_{ij})
## \phi_i \sim N(0,I)
## \psi^2_{jj} \sim IG(a0_j/2, b0_j/2)
##
##
## Andrew D. Martin
## Washington University
##
## Kevin M. Quinn
## Harvard University
##
## May 7, 2003
## revised to accomodate new spec 7/5/2004 KQ
##
##########################################################################

"MCMCfactanal" <-
  function(x, factors, lambda.constraints=list(),
           data=parent.environment(), burnin = 1000, mcmc = 20000,
           thin=1, verbose = FALSE, seed = NA,
           lambda.start = NA, psi.start = NA,
           l0=0, L0=0, a0=0.001, b0=0.001,
           store.scores = FALSE, std.var=TRUE, ... ) {

    ## check for an offset
    check.offset(list(...))
    
    ## get data matrix and associated constants
    if (is.matrix(x)){
      X <- x
      xvars <- dimnames(X)[[2]]
      xobs <- dimnames(X)[[1]]
      N <- nrow(X)
      K <- ncol(X)     
    }
    else {
      holder <- parse.formula(formula=x, data=data,
                              intercept=FALSE, justX=TRUE)
      X <- holder[[2]]
      xvars <- holder[[3]]
      xobs <- holder[[4]]
      N <- nrow(X)
      K <- ncol(X)
    }

    ## standardize X
    if (std.var){
      for (i in 1:K){
        X[,i] <- (X[,i]-mean(X[,i]))/sd(X[,i])
      }
    }
    else{
      for (i in 1:K){
        X[,i] <- X[,i]-mean(X[,i])
      }
    }

    ## take care of the case where X has no row names
    if (is.null(xobs)){
      xobs <- 1:N
    }

    ## check mcmc parameters
    check.mcmc.parameters(burnin, mcmc, thin)

    # seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]

    ## setup constraints on Lambda
    holder <- build.factor.constraints(lambda.constraints, X, K, factors)
    Lambda.eq.constraints <- holder[[1]]
    Lambda.ineq.constraints <- holder[[2]]
    X.names <- holder[[3]]
  
    ## setup prior on Lambda
    holder <- form.factload.norm.prior(l0, L0, K, factors, X.names)
    Lambda.prior.mean <- holder[[1]]
    Lambda.prior.prec <- holder[[2]]

    ## setup and check prior on Psi
    holder <- form.ig.diagmat.prior(a0, b0, K)
    a0 <- holder[[1]]
    b0 <- holder[[2]]
    
    ## starting values for Lambda
    Lambda <- factload.start(lambda.start, K, factors,
                             Lambda.eq.constraints, Lambda.ineq.constraints)
    
    ## starting values for Psi
    Psi <- factuniqueness.start(psi.start, X)
    

    ## define holder for posterior density sample
    if(store.scores == FALSE) {
      sample <- matrix(data=0, mcmc/thin, K*factors+K)
    }
    else {
      sample <- matrix(data=0, mcmc/thin, K*factors+K+factors*N)
    }
   
    ## call C++ code to do the sampling
    auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCfactanal",
                     sample=sample, X=X, burnin=as.integer(burnin),
                     mcmc=as.integer(mcmc), thin=as.integer(thin), 
                     lecuyer=as.integer(lecuyer), 
                     seedarray=as.integer(seed.array),
                     lecuyerstream=as.integer(lecuyer.stream),
                     verbose=as.integer(verbose),
                     Lambda=Lambda, Psi=Psi, Lameq=Lambda.eq.constraints,
                     Lamineq=Lambda.ineq.constraints,
                     Lampmean=Lambda.prior.mean, Lampprec=Lambda.prior.prec,
                     a0=a0, b0=b0, storescores=as.integer(store.scores))

    Lambda.names <- paste(paste("Lambda",
                                rep(X.names,
                                    each=factors), sep=""),
                          rep(1:factors,K), sep="_")
    Psi.names <- paste("Psi", X.names, sep="")
    par.names <- c(Lambda.names, Psi.names)
    if (store.scores==TRUE){
      phi.names <- paste(paste("phi",
                               rep(xobs, each=factors), sep="_"),
                         rep(1:factors,factors), sep="_")
      par.names <- c(par.names, phi.names)
    }
    title <- "MCMCpack Factor Analysis Posterior Density Sample"
    output <- form.mcmc.object(posterior, par.names, title)
    ## get rid of columns for constrained parameters
    output.df <- as.data.frame(as.matrix(output))
    output.var <- diag(var(output.df))
    output.df <- output.df[,output.var != 0]
    
    output <- mcmc(as.matrix(output.df), start=1, end=mcmc, thin=thin)

    ## add constraint info so this isn't lost
    attr(output, "constraints") <- lambda.constraints
    attr(output, "n.manifest") <- K
    attr(output, "n.factors") <- factors
    return(output)
  }
# sample from the posterior distribution of Wakefield's baseline model
# for ecological inference in R using linked C++ code in Scythe
#
# KQ 10/22/2002

"MCMChierEI" <-
  function(r0, r1, c0, c1, burnin=5000, mcmc=50000, thin=1,
           verbose=FALSE, seed=NA,
           m0=0, M0=2.287656,
           m1=0, M1=2.287656,
           a0=0.825, b0=0.0105,
           a1=0.825, b1=0.0105, ...){
    
    # Error checking
    if (length(r0) != length(r1)){
      cat("length(r0) != length(r1).\n")
      stop("Please check data and try MCMChierEI() again.\n")
    }

    if (length(r0) != length(c0)){
      cat("length(r0) != length(c0).\n")
      stop("Please check data and try MCMChierEI() again.\n")
    }
    
    if (length(r0) != length(c1)){
      cat("length(r0) != length(c1).\n")
      stop("Please check data and try MCMChierEI() again.\n")
    }
    
    if (length(r1) != length(c0)){
      cat("length(r1) != length(c0).\n")
      stop("Please check data and try MCMChierEI() again.\n")
    }
    
    if (length(r1) != length(c1)){
      cat("length(r1) != length(c1).\n")
      stop("Please check data and try MCMChierEI() again.\n")
    }
    
    if (length(c0) != length(c1)){
      cat("length(c0) != length(c1).\n")
      stop("Please check data and try MCMChierEI() again.\n")
    }
    
    if (min((r0+r1) == (c0+c1))==0){
      cat("Rows and columns do not sum to same thing.\n")
      stop("Please check data and try MCMChierEI() again.\n")
    }

    check.mcmc.parameters(burnin, mcmc, thin)

    # seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]
   
    
    if (M0 <= 0 ){
      cat("Parameter M0 <= 0.\n")
      stop("Please respecify and try MCMChierEI() again.\n")
    }
    
    if (M1 <= 0 ){
      cat("Parameter M1 <= 0.\n")
      stop("Please respecify and try MCMChierEI() again.\n")
    }
    
    if (a0 <= 0 ){
      cat("Parameter a0 <= 0.\n")
      stop("Please respecify and try MCMChierEI() again.\n")
    }
    
    if (a1 <= 0 ){
      cat("Parameter a1 <= 0.\n")
      stop("Please respecify and try MCMChierEI() again.\n")
    }
    
    if (b0 <= 0 ){
      cat("Parameter b0 <= 0.\n")
      stop("Please respecify and try MCMChierEI() again.\n")
    }
    
    if (b1 <= 0 ){
      cat("Parameter b1 <= 0.\n")
      stop("Please respecify and try MCMChierEI() again.\n")
    }
   
    # setup matrix to hold output from sampling
    ntables = length(r0)
    sample <- matrix(0, mcmc/thin, ntables*2+4)

    # call C++ code to do the sampling
    C.sample <- .C("hierEI",
                   samdata = as.double(sample),
                   samrow = as.integer(nrow(sample)),
                   samcol = as.integer(ncol(sample)),
                   r0 = as.double(r0),
                   r1 = as.double(r1),
                   c0 = as.double(c0),
                   c1 = as.double(c1),
                   ntables = as.integer(ntables),
                   burnin = as.integer(burnin),
                   mcmc = as.integer(mcmc),
                   thin = as.integer(thin),
                   mu0priormean = as.double(m0),
                   mu0priorvar = as.double(M0),
                   mu1priormean = as.double(m1),
                   mu1priorvar = as.double(M1),
                   a0 = as.double(a0),
                   b0 = as.double(b0),
                   a1 = as.double(a1),
                   b1 = as.double(b1),
                   verbose = as.integer(verbose),
                   lecuyer = as.integer(lecuyer),
                   seedarray = as.integer(seed.array),
                   lecuyerstream = as.integer(lecuyer.stream),
                   PACKAGE="MCMCpack"
                   )

    sample <- matrix(C.sample$samdata, C.sample$samrow, C.sample$samcol,
                     byrow=TRUE)
    
    output <- mcmc(data=sample, start=1, end=mcmc, thin=thin)
    p0names <- paste("p0table", 1:ntables, sep="")
    p1names <- paste("p1table", 1:ntables, sep="")
    varnames(output) <- c(p0names, p1names, "mu0", "mu1", "sigma^2.0",
                          "sigma^2.1")
    
    attr(output, "title") <- "MCMCpack Wakefield's Hierarchical EI Model Posterior Density Sample"
        
    return(output)
    
  }
## sample from the posterior distribution of a one-dimensional item
## response theory model in R using linked C++ code in Scythe.
##
## ADM and KQ 1/23/2003
## updated extensively ADM & KQ 7/28/2004
 
"MCMCirt1d" <-
  function(datamatrix, theta.constraints=list(), burnin = 1000,
           mcmc = 20000, thin=1, verbose = FALSE, seed = NA,
           theta.start = NA, alpha.start = NA, beta.start = NA, t0 = 0,
           T0 = 1, ab0=0, AB0=.25, store.item = FALSE, ... ) {
    
    ## checks
    check.offset(list(...))
    check.mcmc.parameters(burnin, mcmc, thin)

    ## check vote matrix and convert to work with C++ code 
    datamatrix <- as.matrix(datamatrix)   
    K <- ncol(datamatrix)   # cases, bills, items, etc
    J <- nrow(datamatrix)   # justices, legislators, subjects, etc
    if(sum(datamatrix==1 | datamatrix==0 | is.na(datamatrix)) != (J*K)) {
      cat("Error: Data matrix contains elements other than 0, 1 or NA.\n")
      stop("Please check data and try MCMCirt1d() again.\n",
              call.=FALSE)
    }
    datamatrix[is.na(datamatrix)] <- 9   
    item.names <- colnames(datamatrix)
    subject.names <- rownames(datamatrix)
    
    ## setup constraints on theta
    for (i in 1:length(theta.constraints)){
      theta.constraints[[i]] <- list(as.integer(1), theta.constraints[[i]][1])
    }
    holder <- build.factor.constraints(theta.constraints, t(datamatrix), J, 1)
    theta.eq.constraints <- holder[[1]]
    theta.ineq.constraints <- holder[[2]]
    subject.names <- holder[[3]]
    ## names
    item.names <- colnames(datamatrix)
    if (is.null(item.names)){
      item.names <- paste("item", 1:K, sep="")
    }

    ## prior for theta 
    holder <- form.mvn.prior(t0, T0, 1)
    t0 <- holder[[1]]
    T0 <- holder[[2]]

    ## prior for (alpha, beta)
    holder <- form.mvn.prior(ab0, AB0, 2)
    ab0 <- holder[[1]]
    AB0 <- holder[[2]]
    
    ## starting values for theta error checking
    theta.start <- factor.score.start.check(theta.start, datamatrix,
                                            t0, T0,
                                            theta.eq.constraints,
                                            theta.ineq.constraints, 1)
    
    ## starting values for (alpha, beta)
    ab.starts <- matrix(NA, K, 2)
    for (i in 1:K){
      local.y <-  datamatrix[,i]
      local.y[local.y==9] <- NA
      if (var(na.omit(local.y))==0){
        ab.starts[i,] <- c(0,10)
      }
      else {
        ab.starts[i,] <- coef(suppressWarnings(glm(local.y~theta.start,
                                                   family=binomial(probit),
                                                   control=glm.control(
                                                     maxit=8, epsilon=1e-3)
                                                   )))
      }
    }
    ab.starts[,1] <- -1 * ab.starts[,1] # make this into a difficulty param
 
    ## starting values for alpha and beta error checking
    if (is.na(alpha.start)) {
      alpha.start <- ab.starts[,1]
    }
    else if(is.null(dim(alpha.start))) {
      alpha.start <- alpha.start * matrix(1,K,1)  
    }
    else if((dim(alpha.start)[1] != K) || (dim(alpha.start)[2] != 1)) {
      cat("Error: Starting value for alpha not conformable.\n")
      stop("Please respecify and call MCMCirt1d() again.\n",
           call.=FALSE)
    }      
    if (is.na(beta.start)) {
      beta.start <- ab.starts[,2]
    }
    else if(is.null(dim(beta.start))) {
      beta.start <- beta.start * matrix(1,K,1)  
    }
    else if((dim(beta.start)[1] != K) || (dim(beta.start)[2] != 1)) {
      cat("Error: Starting value for beta not conformable.\n")
      stop("Please respecify and call MCMCirt1d() again.\n",
              call.=FALSE)
    }    

    ## define holder for posterior density sample
    if(store.item == FALSE) {
      sample <- matrix(data=0, mcmc/thin, J)
    }
    else {
      sample <- matrix(data=0, mcmc/thin, J + 2 * K)
    }

    ## seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]
    
    # call C++ code to draw sample
    posterior <- .C("MCMCirt1d",
                    sampledata = as.double(sample),
                    samplerow = as.integer(nrow(sample)),
                    samplecol = as.integer(ncol(sample)),
                    Xdata = as.integer(datamatrix),
                    Xrow = as.integer(nrow(datamatrix)),
                    Xcol = as.integer(ncol(datamatrix)),    
                    burnin = as.integer(burnin),
                    mcmc = as.integer(mcmc),
                    thin = as.integer(thin),
                    lecuyer = as.integer(lecuyer),
                    seedarray = as.integer(seed.array),
                    lecuyerstream = as.integer(lecuyer.stream),
                    verbose = as.integer(verbose),
                    thetastartdata = as.double(theta.start),
                    thetastartrow = as.integer(nrow(theta.start)),
                    thetastartcol = as.integer(ncol(theta.start)),
                    astartdata = as.double(alpha.start),
                    astartrow = as.integer(length(alpha.start)),
                    astartcol = as.integer(1),
                    bstartdata = as.double(beta.start),
                    bstartrow = as.integer(length(beta.start)),
                    bstartcol = as.integer(1),
                    t0 = as.double(t0),
                    T0 = as.double(T0),
                    ab0data = as.double(ab0),
                    ab0row = as.integer(nrow(ab0)),
                    ab0col = as.integer(ncol(ab0)),
                    AB0data = as.double(AB0),
                    AB0row = as.integer(nrow(AB0)),
                    AB0col = as.integer(ncol(AB0)),
                    thetaeqdata = as.double(theta.eq.constraints),
                    thetaeqrow = as.integer(nrow(theta.eq.constraints)),
                    thetaeqcol = as.integer(ncol(theta.eq.constraints)),
                    thetaineqdata = as.double(theta.ineq.constraints),
                    thetaineqrow = as.integer(nrow(theta.ineq.constraints)),
                    thetaineqcol = as.integer(ncol(theta.ineq.constraints)),
                    store = as.integer(store.item),
                    PACKAGE="MCMCpack"
                  )
    
    theta.names <- paste("theta.", subject.names, sep = "")
    alpha.beta.names <- paste(rep(c("alpha.","beta."), K),
                              rep(item.names, each = 2),
                              sep = "")
   
    # put together matrix and build MCMC object to return
    sample <- matrix(posterior$sampledata, posterior$samplerow,
                     posterior$samplecol,
                     byrow=TRUE)
    output <- mcmc(data=sample, start=1, end=mcmc, thin=thin)
    
    if(store.item == FALSE) {
      names <- theta.names
    }
    else {
      names <- c(theta.names, alpha.beta.names)
    }
    varnames(output) <- names
    attr(output,"title") <-
      "MCMCirt1d Posterior Density Sample"
    return(output)
    
  }
##########################################################################
## sample from a K-dimensional two-parameter item response model with
## probit link. This is just a wrapper function that calls
## MCMCordfactanal.
##
## Andrew D. Martin
## Washington University
##
## Kevin M. Quinn
## Harvard University
##
## June 8, 2003
##
##########################################################################

"MCMCirtKd" <-
  function(datamatrix, dimensions, item.constraints=list(),
           burnin = 1000, mcmc = 10000,
           thin=1, verbose = FALSE, seed = NA,
           alphabeta.start = NA, b0=0, B0=0,
           store.item=FALSE, store.ability=TRUE,
           drop.constantvars=TRUE, ... ) {

    datamatrix <- t(as.matrix(datamatrix))   
    
    post <- MCMCordfactanal(x=datamatrix, factors=dimensions,
                            lambda.constraints=item.constraints,
                            burnin=burnin, mcmc=mcmc, thin=thin,
                            tune=NA, verbose=verbose, seed=seed,
                            lambda.start=alphabeta.start,
                            l0=b0, L0=B0, store.lambda=store.item,
                            store.scores=store.ability,
                            drop.constantvars=drop.constantvars,
                            model="MCMCirtKd")
    return(post)
  }

## sample from the posterior distribution of a logistic regression
## model in R using linked C++ code in Scythe
##
## KQ 1/23/2003
##
## Modified to meet new developer specification 7/15/2004 KQ
## Modified for new Scythe and rngs 7/25/2004 KQ
## note: B0 is now a precision

"MCMClogit" <-
  function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
           thin=1, tune=1.1, verbose = FALSE, seed = NA, beta.start = NA,
           b0 = 0, B0 = 0, ...) {

    ## checks
    check.offset(list(...))
    check.mcmc.parameters(burnin, mcmc, thin)
    
    ## seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]

    ## form response and model matrices
    holder <- parse.formula(formula, data)
    Y <- holder[[1]]
    X <- holder[[2]]
    xnames <- holder[[3]]    
    K <- ncol(X)  # number of covariates
        
    ## starting values and priors
    beta.start <- coef.start(beta.start, K, formula, family=binomial, data)
    mvn.prior <- form.mvn.prior(b0, B0, K)
    b0 <- mvn.prior[[1]]
    B0 <- mvn.prior[[2]]

    ## form the tuning parameter
    tune <- vector.tune(tune, K)
    V <- vcov(glm(formula=formula, data=data, family=binomial))
      
    ## y \in {0, 1} error checking
    if (sum(Y!=0 & Y!=1) > 0) {
      cat("Elements of Y equal to something other than 0 or 1.\n")
      stop("Check data and call MCMClogit() again. \n") 
    }
   
   
    ## define holder for posterior density sample
    sample <- matrix(data=0, mcmc/thin, dim(X)[2] )

    ## call C++ code to draw sample
    auto.Scythe.call(output.object="posterior", cc.fun.name="MCMClogit",
                     sample=sample, Y=Y, X=X, burnin=as.integer(burnin),
                     mcmc=as.integer(mcmc), thin=as.integer(thin),
                     tune=tune, lecuyer=as.integer(lecuyer),
                     seedarray=as.integer(seed.array),
                     lecuyerstream=as.integer(lecuyer.stream),
                     verbose=as.integer(verbose), betastart=beta.start,
                     b0=b0, B0=B0, V=V) 
    
    ## put together matrix and build MCMC object to return
    output <- form.mcmc.object(posterior, names=xnames,
                               title="MCMClogit Posterior Density Sample")
    return(output)    
  }

##########################################################################


## samples from a user-written posterior code in R using a
## random walk Metropolis algorithm
##
## KQ 6/24/2004
##

"MCMCmetrop1R" <- function(fun, theta.init,
                           burnin=500, mcmc=20000, thin=1,
                           tune=1, verbose=TRUE, seed=NA, logfun=TRUE,
                           optim.trace=0, optim.REPORT=10,
                           optim.maxit=500, ...){
  
  ## error checking here
  check.offset(list(...))
  check.mcmc.parameters(burnin, mcmc, thin)
    
  ## form the tuning vector
  tune <- vector.tune(tune, length(theta.init))
  
  ## form seed
  seeds <- form.seeds(seed) 
  lecuyer <- seeds[[1]]
  seed.array <- seeds[[2]]
  lecuyer.stream <- seeds[[3]]
  
  
  ## setup the environment so that fun can see the things passed as ...
  ## it should be the case that users can specify arguments with the same
  ## names as variables defined in MCMCmetrop1R without causing problems.
  my.env <- new.env()
  environment(fun) <- my.env
  dots <- list(...)
  dotnames <- names(dots)
  ndots <- length(dots)
  if (ndots >= 1){
    for (i in 1:ndots){
      assign(x=dotnames[[i]], value=dots[[i]], inherits=FALSE, envir=my.env)
    }
  }
  
  ## find approx mode and Hessian using optim()
  opt.out <- optim(theta.init, fun,
                   control=list(fnscale=-1, trace=optim.trace,
                     REPORT=optim.REPORT, maxit=optim.maxit),
                   method="BFGS", hessian=TRUE)
  if(opt.out$convergence!=0){
    warning("Mode and Hessian were not found with call to optim().\nSampling proceeded anyway. \n") 
  }
  
  propvar <- tune %*% solve(-1*opt.out$hessian) %*% tune
  
  ## Call the C++ function to do the MCMC sampling 
  sample <- .Call("MCMCmetrop1R_cc", fun, as.double(theta.init),
                  my.env, as.integer(burnin), as.integer(mcmc),
                  as.integer(thin),
                  as.logical(verbose),
                  lecuyer=as.integer(lecuyer), 
                  seedarray=as.integer(seed.array),
                  lecuyerstream=as.integer(lecuyer.stream),
                  as.logical(logfun),
                  as.matrix(propvar),
                  PACKAGE="MCMCpack")

  ## turn sample into an mcmc object
  sample <- mcmc(data=sample, start=1, end = mcmc, thin=thin)
  return(sample)
}
 
##########################################################################
## sample from the posterior distribution of a factor analysis model
## model in R using linked C++ code in Scythe.
##
## The model is:
##
## x*_i = \Lambda \phi_i + \epsilon_i,   \epsilon_i \sim N(0, \Psi)
##
## \lambda_{ij} \sim N(l0_{ij}, L0^{-1}_{ij})
## \phi_i \sim N(0,I)
##
## and x*_i is the latent variable formed from the observed ordinal
## variable in the usual (Albert and Chib, 1993) way and is equal to
## x_i when x_i is continuous. When x_j is ordinal \Psi_jj is assumed
## to be 1.
##
## Andrew D. Martin
## Washington University
##
## Kevin M. Quinn
## Harvard University
##
## 12/2/2003
## Revised to accommodate new spec 7/20/2004
##
##########################################################################

"MCMCmixfactanal" <-
  function(x, factors, lambda.constraints=list(),
           data=parent.environment(), burnin = 1000, mcmc = 20000,
           thin=1, tune=NA, verbose = FALSE, seed = NA,
           lambda.start = NA, psi.start=NA,
           l0=0, L0=0, a0=0.001, b0=0.001,
           store.lambda=TRUE, store.scores=FALSE,
           std.mean=TRUE, std.var=TRUE, ... ) {
    
    call <- match.call()
    mt <- terms(x, data=data)
    if (attr(mt, "response") > 0) 
      stop("Response not allowed in formula in MCMCmixfactanal().\n")
    if(missing(data)) data <- sys.frame(sys.parent())
    mf <- match.call(expand.dots = FALSE)
    mf$factors <- mf$lambda.constraints <- mf$burnin <- mf$mcmc <- NULL
    mf$thin <- mf$tune <- mf$verbose <- mf$seed <- NULL
    mf$lambda.start <- mf$l0 <- mf$L0 <- mf$a0 <- mf$b0 <- NULL
    mf$store.lambda <- mf$store.scores <- mf$std.var <- mf$... <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    attributes(mt)$intercept <- 0
    Xterm.length <- length(attr(mt, "variables"))
    X <- subset(mf,
                select=as.character(attr(mt, "variables"))[2:Xterm.length])

    N <- nrow(X)	      # number of observations      
    K <- ncol(X)              # number of manifest variables
    ncat <- matrix(NA, K, 1)  # vector of number of categ. in each man. var. 
    for (i in 1:K){ 
      if (is.numeric(X[,i])){
        ncat[i] <- -999
        X[is.na(X[,i]), i] <- -999
      }
      else if (is.ordered(X[,i])){
        ncat[i] <- nlevels(X[,i])      
        X[,i] <- as.integer(X[,i])
        X[is.na(X[,i]), i] <- -999
      }
      else {
        stop("Manifest variable ", dimnames(X)[[2]][i],
             " neither ordered factor nor numeric variable.\n")
      }
    }
    
    X <- as.matrix(X)
    xvars <- dimnames(X)[[2]] # X variable names
    xobs <- dimnames(X)[[1]]  # observation names
    
    if (is.null(xobs)){
      xobs <- 1:N
    }

    # standardize X
    if (std.mean){
      for (i in 1:K){
        if (ncat[i] == -999){
          X[,i] <- X[,i]-mean(X[,i])
        }
      }
    }
    if (std.var){
      for (i in 1:K){
        if (ncat[i] == -999){
          X[,i] <- (X[,i])/sd(X[,i])
        }
      }
    }
    
    n.ord.ge3 <- 0
    for (i in 1:K)
      if (ncat[i] >= 3) n.ord.ge3 <- n.ord.ge3 + 1
   
    check.mcmc.parameters(burnin, mcmc, thin)
    
    ## setup constraints on Lambda
    holder <- build.factor.constraints(lambda.constraints, X, K, factors+1)
    Lambda.eq.constraints <- holder[[1]]
    Lambda.ineq.constraints <- holder[[2]]
    X.names <- holder[[3]]
    
    ## if subtracting out the mean of continuous X then constrain
    ## the mean parameter to 0
    for (i in 1:K){
      if (ncat[i] < 2 && std.mean==TRUE){
        if ((Lambda.eq.constraints[i,1] == -999 ||
             Lambda.eq.constraints[i,1] == 0.0) &&
            Lambda.ineq.constraints[i,1] == 0.0){
          Lambda.eq.constraints[i,1] <- 0.0
        }
        else {
          cat("Constraints on Lambda are logically\ninconsistent with std.mean==TRUE.\n")
          stop("Please respecify and call MCMCmixfactanal() again\n")          
        }
      }
    }    


    ## setup and check prior on Psi
    holder <- form.ig.diagmat.prior(a0, b0, K)
    a0 <- holder[[1]]
    b0 <- holder[[2]]

    ## setup prior on Lambda
    holder <- form.factload.norm.prior(l0, L0, K, factors+1, X.names)
    Lambda.prior.mean <- holder[[1]]
    Lambda.prior.prec <- holder[[2]]

    # seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]
    
    # Starting values for Lambda
    Lambda <- matrix(0, K, factors+1)
    if (is.na(lambda.start)){# sets Lambda to equality constraints & 0s
      for (i in 1:K){
        for (j in 1:(factors+1)){
          if (Lambda.eq.constraints[i,j]==-999){
            if(Lambda.ineq.constraints[i,j]==0){
              if (j==1){
                if (ncat[i] < 2){
                  Lambda[i,j] <- mean(X[,i]!=-999)
                }
                if (ncat[i] == 2){
                  probit.out <- glm(as.factor(X[X[,i]!=-999,i])~1,
                                    family=binomial(link=probit))
                  probit.beta <- coef(probit.out)
                  Lambda[i,j] <- probit.beta[1]
                }
                if (ncat[i] > 2){
                  polr.out <- polr(ordered(X[X[,i]!=-999,i])~1)
                  Lambda[i,j] <- -polr.out$zeta[1]*.588
                }
              }
            }
            if(Lambda.ineq.constraints[i,j]>0){
              Lambda[i,j] <- 1.0
            }
            if(Lambda.ineq.constraints[i,j]<0){
              Lambda[i,j] <- -1.0
            }          
          }
          else Lambda[i,j] <- Lambda.eq.constraints[i,j]
        }
      }    
    }
    else if (is.matrix(lambda.start)){
      if (nrow(lambda.start)==K && ncol(lambda.start)==(factors+1))
        Lambda  <- lambda.start
      else {
        cat("Starting values not of correct size for model specification.\n")
        stop("Please respecify and call ", echo.name, "() again\n")
      }
    }
    else if (length(lambda.start)==1 && is.numeric(lambda.start)){
      Lambda <- matrix(lambda.start, K, factors+1)
      for (i in 1:K){
        for (j in 1:(factors+1)){
          if (Lambda.eq.constraints[i,j] != -999)
            Lambda[i,j] <- Lambda.eq.constraints[i,j]
        }
      }    
    }
    else {
      cat("Starting values neither NA, matrix, nor scalar.\n")
      stop("Please respecify and call ", echo.name, "() again\n")
    }
    
    # check MH tuning parameter
    if (is.na(tune)){
      tune <- matrix(NA, K, 1)
      for (i in 1:K){
        tune[i] <- abs(0.05/ncat[i])
      }
    }
    else if (is.double(tune)){
      tune <- matrix(abs(tune/ncat), K, 1)
    }
  
    # starting values for gamma (note: not changeable by user)
    if (max(ncat) <= 2){
      gamma <- matrix(0, 3, K)
    }
    else {
      gamma <- matrix(0, max(ncat)+1, K)
    }
    for (i in 1:K){
      if (ncat[i]<=2){
        gamma[1,i] <- -300
        gamma[2,i] <- 0
        gamma[3,i] <- 300
      }
      if(ncat[i] > 2) {
        polr.out <- polr(ordered(X[X[,i]!=-999,i])~1)
        gamma[1,i] <- -300
        gamma[2,i] <- 0
        gamma[3:ncat[i],i] <- (polr.out$zeta[2:(ncat[i]-1)] -
                               polr.out$zeta[1])*.588
        gamma[ncat[i]+1,i] <- 300
      }
    }

    ## starting values for Psi
    Psi <- factuniqueness.start(psi.start, X)
    for (i in 1:K){
      if (ncat[i] >= 2){
        Psi[i,i] <- 1.0
      }
    }
          
    # define holder for posterior density sample
    if (store.scores == FALSE && store.lambda == FALSE){
      sample <- matrix(data=0, mcmc/thin, length(gamma)+K)
    }
    else if (store.scores == TRUE && store.lambda == FALSE){
      sample <- matrix(data=0, mcmc/thin, (factors+1)*N + length(gamma)+K)
    }
    else if(store.scores == FALSE && store.lambda == TRUE) {
      sample <- matrix(data=0, mcmc/thin, K*(factors+1)+length(gamma)+K)
    }
    else { # store.scores==TRUE && store.lambda==TRUE
      sample <- matrix(data=0, mcmc/thin, K*(factors+1)+(factors+1)*N +
                       length(gamma)+K)
    }
    
    # Call the C++ code to do the real work
    posterior <- .C("mixfactanalpost",
                    samdata = as.double(sample),
                    samrow = as.integer(nrow(sample)),
                    samcol = as.integer(ncol(sample)),
                    X = as.double(X),
                    Xrow = as.integer(nrow(X)),
                    Xcol = as.integer(ncol(X)),
                    burnin = as.integer(burnin),
                    mcmc = as.integer(mcmc),
                    thin = as.integer(thin),
                    tune = as.double(tune),
                    lecuyer = as.integer(lecuyer),
                    seedarray = as.integer(seed.array),
                    lecuyerstream = as.integer(lecuyer.stream),
                    verbose = as.integer(verbose),
                    Lambda = as.double(Lambda),
                    Lambdarow = as.integer(nrow(Lambda)),
                    Lambdacol = as.integer(ncol(Lambda)),
                    gamma = as.double(gamma),
                    gammarow = as.integer(nrow(gamma)),
                    gammacol = as.integer(ncol(gamma)),
                    Psi = as.double(Psi),
                    Psirow = as.integer(nrow(Psi)),
                    Psicol = as.integer(ncol(Psi)),
                    ncat = as.integer(ncat),
                    ncatrow = as.integer(nrow(ncat)),
                    ncatcol = as.integer(ncol(ncat)),
                    Lameq = as.double(Lambda.eq.constraints),
                    Lameqrow = as.integer(nrow(Lambda.eq.constraints)),
                    Lameqcol = as.integer(ncol(Lambda.ineq.constraints)),
                    Lamineq = as.double(Lambda.ineq.constraints),
                    Lamineqrow = as.integer(nrow(Lambda.ineq.constraints)),
                    Lamineqcol = as.integer(ncol(Lambda.ineq.constraints)),
                    Lampmean = as.double(Lambda.prior.mean),
                    Lampmeanrow = as.integer(nrow(Lambda.prior.mean)),
                    Lampmeancol = as.integer(ncol(Lambda.prior.prec)),
                    Lampprec = as.double(Lambda.prior.prec),
                    Lampprecrow = as.integer(nrow(Lambda.prior.prec)),
                    Lamppreccol = as.integer(ncol(Lambda.prior.prec)),
                    a0 = as.double(a0),
                    a0row = as.integer(nrow(a0)),
                    a0col = as.integer(ncol(a0)),
                    b0 = as.double(b0),
                    b0row = as.integer(nrow(b0)),
                    b0col = as.integer(ncol(b0)),
                    storelambda = as.integer(store.lambda),
                    storescores = as.integer(store.scores),
                    accepts = as.integer(0),
                    PACKAGE="MCMCpack"
                    )

    cat(" overall acceptance rate = ",
        posterior$accepts / ((posterior$burnin+posterior$mcmc)*n.ord.ge3), "\n")

    
    # put together matrix and build MCMC object to return
    sample <- matrix(posterior$samdata, posterior$samrow, posterior$samcol,
                     byrow=TRUE)
    output <- mcmc(data=sample,start=1, end=mcmc, thin=thin)
    
    par.names <- NULL
    if (store.lambda==TRUE){
      Lambda.names <- paste(paste("Lambda",
                                  rep(X.names,
                                      each=(factors+1)), sep=""),
                            rep(1:(factors+1),K), sep=".")
      par.names <- c(par.names, Lambda.names)
    }
    
    gamma.names <- paste(paste("gamma",
                               rep(0:(nrow(gamma)-1),
                                   each=K), sep=""),
                         rep(X.names,  nrow(gamma)), sep=".")
    par.names <- c(par.names, gamma.names)
    
    if (store.scores==TRUE){
      phi.names <- paste(paste("phi",
                               rep(xobs, each=(factors+1)), sep="."),
                         rep(1:(factors+1),(factors+1)), sep=".")
      par.names <- c(par.names, phi.names)
    }

    Psi.names <- paste("Psi", X.names, sep=".")
    par.names <- c(par.names, Psi.names)
    
    varnames(output) <- par.names

    # get rid of columns for constrained parameters
    output.df <- as.data.frame(as.matrix(output))
    output.var <- diag(var(output.df))
    output.df <- output.df[,output.var != 0]
    output <- mcmc(as.matrix(output.df), start=1, end=mcmc, thin=thin)
    
    # add constraint info so this isn't lost
    attr(output, "constraints") <- lambda.constraints
    attr(output, "n.manifest") <- K
    attr(output, "n.factors") <- factors
    attr(output, "accept.rate") <- posterior$accepts /
      ((posterior$burnin+posterior$mcmc)*n.ord.ge3)
      attr(output,"title") <-
        "MCMCpack Mixed Data Factor Analysis Posterior Density Sample"
    
    return(output)
    
  }

## sample from the posterior distribution of an ordered probit model
## via the data augmentation approach of Cowles (1996)
##
## KQ 1/25/2003
## Modified to meet new developer specification 7/26/2004 KQ
## Modified for new Scythe and rngs 7/26/2004 KQ


"MCMCoprobit" <-
  function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
           thin = 1, tune = NA, verbose = FALSE, seed = NA, beta.start = NA,
           b0 = 0, B0 = 0, ...) {


    ## checks
    check.offset(list(...))
    check.mcmc.parameters(burnin, mcmc, thin)
    
    ## seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]

    
    ## extract X, Y, and variable names from the model formula and frame
    call <- match.call()
    mt <- terms(formula, data=data)
    if(missing(data)) data <- sys.frame(sys.parent())
    mf <- match.call(expand.dots = FALSE)
    mf$burnin <- mf$mcmc <- mf$b0 <- mf$B0 <- NULL
    mf$thin <- mf$... <- mf$tune <- mf$verbose <- mf$seed <- NULL
    mf$beta.start <- NULL
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, sys.frame(sys.parent()))
    vars <- as.character(attr(mt, "variables"))[-1] # y varname and x varnames
    
    ## null model support
    X <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts)# else NULL
    X.names <- dimnames(X)[[2]]
    Y    <- model.response(mf, "numeric")
    Y    <- factor(Y, ordered=TRUE)
    ncat <- nlevels(Y)             # number of categories of y
    cat  <- levels(Y)              # values of categories of y
    N <- nrow(X)	                  # number of observations
    K <- ncol(X)	                  # number of covariates
    if (length(Y) != N){
      cat("X and Y do not have same number of rows.\n")
      stop("Please respecify and call MCMCoprobit() again.\n")      
    }
    
    ## convert data to matrices to be passed
    Y <- as.matrix(as.integer(Y))
    X <- as.matrix(X)
    
    ## check tuning parameter
    if (is.na(tune)){
      tune <- 0.05/ncat
    }
    
    xint <- match("(Intercept)", colnames(X), nomatch=0)
    if (xint > 0){
      new.X <- X[, -xint, drop=FALSE]
    }
    else warning("An intercept is needed and assumed in MCMCoprobit()\n.")
    if (ncol(new.X) == 0){
      polr.out <- polr(ordered(Y)~1)
    }
    else {
      polr.out <- polr(ordered(Y)~new.X)
    }
    
    ## starting values for beta error checking
    if (is.na(beta.start)){
      beta.start <- matrix(0, K, 1)
      beta.start[1] <- -.588 * polr.out$zeta[1]
      if( ncol(new.X) > 0){
        beta.start[2:K] <- .588 * coef(polr.out)
      }
    }
    else if(is.null(dim(beta.start))) {
      beta.start <- beta.start * matrix(1,K,1)  
    }
    else if((dim(beta.start)[1] != K) || (dim(beta.start)[2] != 1)) {
      cat("Starting value for beta not conformable.\n")
      stop("Please respecify and call MCMCoprobit() again.\n")
    }
    
    ## prior for beta error checking
    if(is.null(dim(b0))) {
      b0 <- b0 * matrix(1,K,1)  
    }
    if((dim(b0)[1] != K) || (dim(b0)[2] != 1)) {
      cat("N(b0,B0) prior b0 not conformable.\n")
      stop("Please respecify and call MCMCoprobit() again.\n") 
    }   
    if(is.null(dim(B0))) {
      B0 <- B0 * diag(K)    
    }
    if((dim(B0)[1] != K) || (dim(B0)[2] != K)) {
      cat("N(b0,B0) prior B0 not conformable.\n")
      stop("Please respecify and call MCMCoprobit() again.\n")
    }
    
    ## form gamma starting values (note: not changeable)
    gamma <- matrix(NA,ncat+1,1)
    gamma[1] <- -300
    gamma[2] <- 0
    gamma[3:ncat] <- (polr.out$zeta[2:(ncat-1)] - polr.out$zeta[1])*.588
    gamma[ncat+1] <- 300
    
    ## posterior density sample
    sample <- matrix(data=0, mcmc/thin, K + ncat + 1)

    ## call C++ code to draw sample
    auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCoprobit",
                     sample=sample, Y=as.integer(Y), X=X, 
                     burnin=as.integer(burnin),
                     mcmc=as.integer(mcmc), thin=as.integer(thin),
                     tune=as.double(tune), lecuyer=as.integer(lecuyer),
                     seedarray=as.integer(seed.array),
                     lecuyerstream=as.integer(lecuyer.stream),
                     verbose=as.integer(verbose), beta=beta.start,
                     gamma=gamma, b0=b0, B0=B0) 
    
    ## put together matrix and build MCMC object to return
    sample <- matrix(posterior$sampledata, posterior$samplerow,
                     posterior$samplecol, byrow=TRUE)
    sample <- sample[,c(1:K, (K+3):(K+ncat))]
    output <- mcmc(data=sample, start=1, end=mcmc, thin=thin)
    xnames <- c(X.names, paste("gamma", 2:(ncat-1), sep=""))
    varnames(output) <- xnames
    attr(output, "title") <- "MCMCoprobit Posterior Density Sample"
        
    return(output)
  }
###########################################################################
## sample from the posterior distribution of a factor analysis model
## model in R using linked C++ code in Scythe.
##
## The model is:
##
## x*_i = \Lambda \phi_i + \epsilon_i,   \epsilon_i \sim N(0, I)
##
## \lambda_{ij} \sim N(l0_{ij}, L0^{-1}_{ij})
## \phi_i \sim N(0,I)
##
## and x*_i is the latent variable formed from the observed ordinal
## variable in the usual (Albert and Chib, 1993) way.
##
## Andrew D. Martin
## Washington University
##
## Kevin M. Quinn
## Harvard University
##
## May 12, 2003
## Revised to accommodate new spec 7/13/2004
## 
##########################################################################

"MCMCordfactanal" <-
  function(x, factors, lambda.constraints=list(),
           data=parent.environment(), burnin = 1000, mcmc = 20000,
           thin=1, tune=NA, verbose = FALSE, seed = NA,
           lambda.start = NA, l0=0, L0=0,
           store.lambda=TRUE, store.scores=FALSE,
           drop.constantvars=TRUE, ... ) {

    ## check for MCMCirtKd special case, this is used to tell the R
    ## and C++ code what to echo (1 if regular, 2 if MCMCirtKd)
    ## the test is based on the existence of model="MCMCirtKd"
    ## passed through ...
    args <- list(...)

    if (length(args$model) > 0){ # if model arg is passed
      if (args$model=="MCMCirtKd"){
        case.switch <- 2
        echo.name <- "MCMCirtKd"
      }
      ## could allow for other possibities here but not clear what they
      ## would be 
    }
    else { # if model arg not passed then assume MCMCordfactanal
      case.switch <- 1
      echo.name <- "MCMCordfactanal"
    }
    
 
    # extract X and variable names from the model formula and frame       
    if (is.matrix(x)){
      if (drop.constantvars==TRUE){
        x.col.var <- apply(x, 2, var, na.rm=TRUE)
        x <- x[,x.col.var!=0]
        x.row.var <- apply(x, 1, var, na.rm=TRUE)
        x <- x[x.row.var!=0,]
      }
      X <- as.data.frame(x)
      xvars <- dimnames(X)[[2]]
      xobs <- dimnames(X)[[1]]
      N <- nrow(X)    # number of observations
      K <- ncol(X)    # number of manifest variables
      ncat <- matrix(NA, K, 1) # vector of number of categ. in each man. var. 
      for (i in 1:K){
        X[,i] <- factor(X[,i], ordered=TRUE)
        ncat[i] <- nlevels(X[,i])
        X[,i] <- as.integer(X[,i])
        X[is.na(X[,i]), i] <- -999
      }
      X <- as.matrix(X)
    }
    else {
      call <- match.call()
      mt <- terms(x, data=data)
      if (attr(mt, "response") > 0) 
        stop("Response not allowed in formula in ", echo.name, "().\n")
      if(missing(data)) data <- sys.frame(sys.parent())
      mf <- match.call(expand.dots = FALSE)
      mf$factors <- mf$lambda.constraints <- mf$burnin <- mf$mcmc <- NULL
      mf$thin <- mf$tune <- mf$verbose <- mf$seed <- NULL
      mf$lambda.start <- mf$l0 <- mf$L0 <- NULL
      mf$store.lambda <- mf$store.scores <- mf$drop.constantvars <- NULL
      mf$... <- NULL
      mf$drop.unused.levels <- TRUE
      mf[[1]] <- as.name("model.frame")
      mf <- eval(mf, sys.frame(sys.parent()))
      attributes(mt)$intercept <- 0
      Xterm.length <- length(attr(mt, "variables"))
      X <- subset(mf,
                  select=as.character(attr(mt, "variables"))[2:Xterm.length])
      if (drop.constantvars==TRUE){
        x.col.var <- apply(X, 2, var, na.rm=TRUE)
        X <- X[,x.col.var!=0]
      }
      N <- nrow(X)	        # number of observations      
      K <- ncol(X)              # number of manifest variables
      ncat <- matrix(NA, K, 1)  # vector of number of categ. in each man. var. 
      for (i in 1:K){
        X[,i] <- factor(X[,i], ordered=TRUE)
        ncat[i] <- nlevels(X[,i])      
        X[,i] <- as.integer(X[,i])
        X[is.na(X[,i]), i] <- -999
      }
      X <- as.matrix(X)
      xvars <- dimnames(X)[[2]] # X variable names
      xobs <- dimnames(X)[[1]]  # observation names
    }

    ## take care of the case where X has no row names
    if (is.null(xobs)){
      xobs <- 1:N
    }
    
    check.offset(list(...))
    check.mcmc.parameters(burnin, mcmc, thin)
    
    ## setup constraints on Lambda
    holder <- build.factor.constraints(lambda.constraints, X, K, factors+1)
    Lambda.eq.constraints <- holder[[1]]
    Lambda.ineq.constraints <- holder[[2]]
    X.names <- holder[[3]]
  
    ## setup prior on Lambda
    holder <- form.factload.norm.prior(l0, L0, K, factors+1, X.names)
    Lambda.prior.mean <- holder[[1]]
    Lambda.prior.prec <- holder[[2]]

    # seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]

    ## Starting values for Lambda
    Lambda <- matrix(0, K, factors+1)
    if (is.na(lambda.start)){# sets Lambda to equality constraints & 0s
      for (i in 1:K){
        for (j in 1:(factors+1)){
          if (Lambda.eq.constraints[i,j]==-999){
            if(Lambda.ineq.constraints[i,j]==0){
              if (j==1){
                if (ncat[i] == 2){
                  probit.out <- glm(as.factor(X[X[,i]!=-999,i])~1,
                                    family=binomial(link=probit))
                  probit.beta <- coef(probit.out)
                  Lambda[i,j] <- probit.beta[1]
                }
                if (ncat[i] > 2){
                  polr.out <- polr(ordered(X[X[,i]!=-999,i])~1)
                  Lambda[i,j] <- -polr.out$zeta[1]*.588
                }
              }
            }
            if(Lambda.ineq.constraints[i,j]>0){
              Lambda[i,j] <- 1.0
            }
            if(Lambda.ineq.constraints[i,j]<0){
              Lambda[i,j] <- -1.0
            }          
          }
          else Lambda[i,j] <- Lambda.eq.constraints[i,j]
        }
      }    
    }
    else if (is.matrix(lambda.start)){
      if (nrow(lambda.start)==K && ncol(lambda.start)==(factors+1))
        Lambda  <- lambda.start
      else {
        cat("Starting values not of correct size for model specification.\n")
        stop("Please respecify and call ", echo.name, "() again\n")
      }
    }
    else if (length(lambda.start)==1 && is.numeric(lambda.start)){
      Lambda <- matrix(lambda.start, K, factors+1)
      for (i in 1:K){
        for (j in 1:(factors+1)){
          if (Lambda.eq.constraints[i,j] != -999)
            Lambda[i,j] <- Lambda.eq.constraints[i,j]
        }
      }    
    }
    else {
      cat("Starting values neither NA, matrix, nor scalar.\n")
      stop("Please respecify and call ", echo.name, "() again\n")
    }

    ## check MH tuning parameter
    if (is.na(tune)){
      tune <- matrix(NA, K, 1)
      for (i in 1:K){
        tune[i] <- 0.05/ncat[i]
      }
    }
    else if (is.double(tune)){
      tune <- matrix(tune/ncat, K, 1)
    }
    if(min(tune) < 0) {
      cat("Tuning parameter is negative.\n")
      stop("Please respecify and call ", echo.name, "() again\n")
    }

    ## starting values for gamma (note: not changeable by user)
    gamma <- matrix(0, max(ncat)+1, K)
    for (i in 1:K){
      if (ncat[i]<=2){
        gamma[1,i] <- -300
        gamma[2,i] <- 0
        gamma[3,i] <- 300
      }
      if(ncat[i] > 2) {
        polr.out <- polr(ordered(X[X[,i]!=-999,i])~1)
        gamma[1,i] <- -300
        gamma[2,i] <- 0
        gamma[3:ncat[i],i] <- (polr.out$zeta[2:(ncat[i]-1)] -
                               polr.out$zeta[1])*.588
        gamma[ncat[i]+1,i] <- 300
      }
    }

    ## define holder for posterior density sample
    if (store.scores == FALSE && store.lambda == FALSE){
      sample <- matrix(data=0, mcmc/thin, length(gamma))
    }
    else if (store.scores == TRUE && store.lambda == FALSE){
      sample <- matrix(data=0, mcmc/thin, (factors+1)*N + length(gamma))
    }
    else if(store.scores == FALSE && store.lambda == TRUE) {
      sample <- matrix(data=0, mcmc/thin, K*(factors+1)+length(gamma))
    }
    else { # store.scores==TRUE && store.lambda==TRUE
      sample <- matrix(data=0, mcmc/thin, K*(factors+1)+(factors+1)*N +
                       length(gamma))
    }

    
    ## Call the C++ code to do the real work
    posterior <- .C("ordfactanalpost",
                    samdata = as.double(sample),
                    samrow = as.integer(nrow(sample)),
                    samcol = as.integer(ncol(sample)),
                    X = as.integer(X),
                    Xrow = as.integer(nrow(X)),
                    Xcol = as.integer(ncol(X)),
                    burnin = as.integer(burnin),
                    mcmc = as.integer(mcmc),
                    thin = as.integer(thin),
                    tune = as.double(tune),
                    lecuyer = as.integer(lecuyer),
                    seedarray = as.integer(seed.array),
                    lecuyerstream = as.integer(lecuyer.stream),
                    verbose = as.integer(verbose),
                    Lambda = as.double(Lambda),
                    Lambdarow = as.integer(nrow(Lambda)),
                    Lambdacol = as.integer(ncol(Lambda)),
                    gamma = as.double(gamma),
                    gammarow = as.integer(nrow(gamma)),
                    gammacol = as.integer(ncol(gamma)),
                    ncat = as.integer(ncat),
                    ncatrow = as.integer(nrow(ncat)),
                    ncatcol = as.integer(ncol(ncat)),
                    Lameq = as.double(Lambda.eq.constraints),
                    Lameqrow = as.integer(nrow(Lambda.eq.constraints)),
                    Lameqcol = as.integer(ncol(Lambda.ineq.constraints)),
                    Lamineq = as.double(Lambda.ineq.constraints),
                    Lamineqrow = as.integer(nrow(Lambda.ineq.constraints)),
                    Lamineqcol = as.integer(ncol(Lambda.ineq.constraints)),
                    Lampmean = as.double(Lambda.prior.mean),
                    Lampmeanrow = as.integer(nrow(Lambda.prior.mean)),
                    Lampmeancol = as.integer(ncol(Lambda.prior.prec)),
                    Lampprec = as.double(Lambda.prior.prec),
                    Lampprecrow = as.integer(nrow(Lambda.prior.prec)),
                    Lamppreccol = as.integer(ncol(Lambda.prior.prec)),
                    storelambda = as.integer(store.lambda),
                    storescores = as.integer(store.scores),
                    accepts = as.integer(0),
                    outswitch = as.integer(case.switch),
                    PACKAGE="MCMCpack"
                    )
    if(case.switch==1) {
      cat(" overall acceptance rate = ",
          posterior$accepts / ((posterior$burnin+posterior$mcmc)*K), "\n")
    }
    
    # put together matrix and build MCMC object to return
    sample <- matrix(posterior$samdata, posterior$samrow, posterior$samcol,
                     byrow=TRUE)
    output <- mcmc(data=sample,start=1, end=mcmc, thin=thin)
    
    par.names <- NULL
    if (store.lambda==TRUE){
      if(case.switch==1) {
      Lambda.names <- paste(paste("Lambda",
                                  rep(X.names,
                                      each=(factors+1)), sep=""),
                            rep(1:(factors+1),K), sep=".")
      }
      if(case.switch==2) {
        alpha.hold <- paste("alpha", X.names, sep=".")
        beta.hold <- paste("beta", X.names, sep = ".")
        beta.hold <- rep(beta.hold, factors, each=factors)
        beta.hold <- paste(beta.hold, 1:factors, sep=".")
                
        Lambda.names <- t(cbind(matrix(alpha.hold, K, 1), 
                                matrix(beta.hold,K,factors,byrow=TRUE)))  
        dim(Lambda.names) <- NULL
      }
      par.names <- c(par.names, Lambda.names)
    }
    
    gamma.names <- paste(paste("gamma",
                               rep(0:(nrow(gamma)-1),
                                   each=K), sep=""),
                         rep(X.names,  nrow(gamma)), sep=".")
    par.names <- c(par.names, gamma.names)
    
    if (store.scores==TRUE){
      if(case.switch==1) {
      phi.names <- paste(paste("phi",
                               rep(xobs, each=(factors+1)), sep="."),
                         rep(1:(factors+1),(factors+1)), sep=".")
      par.names <- c(par.names, phi.names)
      }
      if(case.switch==2) {
      phi.names <- paste(paste("theta",
                               rep(xobs, each=(factors+1)), sep="."),
                         rep(0:factors,(factors+1)), sep=".")
      par.names <- c(par.names, phi.names)      
      
      }
    }
    
    varnames(output) <- par.names

    # get rid of columns for constrained parameters
    output.df <- as.data.frame(as.matrix(output))
    output.var <- diag(var(output.df))
    output.df <- output.df[,output.var != 0]
    output <- mcmc(as.matrix(output.df), start=1, end=mcmc, thin=thin)
    
    # add constraint info so this isn't lost
    attr(output, "constraints") <- lambda.constraints
    attr(output, "n.manifest") <- K
    attr(output, "n.factors") <- factors
    attr(output, "accept.rate") <- posterior$accepts /
      ((posterior$burnin+posterior$mcmc)*K)
    if(case.switch==1) {
      attr(output,"title") <-
        "MCMCpack Ordinal Data Factor Analysis Posterior Density Sample"
    }
    if(case.switch==2) {
      attr(output,"title") <-
        "MCMCpack K-Dimensional Item Response Theory Model Posterior Density Sample"
    }
    return(output)
    
  }

# sample from the posterior distribution of general linear panel
# model in R using linked C++ code in Scythe
#
# NOTE: this does not take a data argument because only
#       matrices are passed.  This should probably be fixed.  Also,
#       re-implementing using something likes Bates' syntax
#       would be nice.  The helper functions could also be used 
#       all over the place here. This is another good project for a grad
#       student.
#
#
# ADM and KQ 8/1/2002
# updated with Ben Goodrich's feedback and new spec ADM 7/28/2004


"MCMCpanel" <-
  function(obs, Y, X, W, burnin = 1000, mcmc = 10000, thin = 5, 
           verbose = FALSE, seed = NA, sigma2.start = NA,
           D.start = NA, b0 = 0, B0 = 1, eta0, R0, nu0 = 0.001,
           delta0 = 0.001, ...) {

    # DESCRIPTION:
    #
    #   MCMCpanel fits a general linear panel model using Algorithm 2 of
    #   Chib and Carlin (1999).  The program calls a compiled C++ shared
    #   library to perform the actual sampling. The model takes the
    #   following form:
    #
    #      y_i = X_i \beta + W_i b_i + \varepsilon_i
    #
    #      b_i \sim N_q(0,D)
    #
    #      \varepsilon_i \sim N_k(0,\sigma^2 I_n)
    #
    #   With conjugate priors:
    #
    #      \beta \sim N_p(\beta_0, \B_0^-1)
    #
    #      D^-1 \sim Wishart(\nu_0^{-1} R_0, \nu_0)
    #
    #      \sigma^-2 \sim Gamma(\nu_00/2, \delta_00/2) 
    #
    #   The model is defined in terms of k (the number of responses
    #   per subject, assumed to be constant across subjects), p (the
    #   number of columns in the design matrix of covariates), and 
    #   q (the number of columns in the design matrix), and n (the
    #   number of subjects).  The components of the model are the
    #   following:
    #
    #      y_i (k \times 1) vector of responses for subject i
    #
    #      X_i (k \times p) matrix of covariates for subject i      
    #
    #      \beta (p \times 1) vector of fixed effects coefficients
    #
    #      W_i (k \times q) design matrix for random effects for subject i
    #
    #      b_i (q \times 1) vector of random effects for subject i
    #
    #      \varepsilon (k \times 1) vector of errors for subject i
    #

    # seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]

    # model parameters
    n <- length(unique(obs))
    k <- length(Y) / n
    p <- dim(X)[2]
    q <- dim(W)[2]

    # check data conformability
    obs.temp <- as.matrix(obs)
    if (any(obs.temp[,1] != obs)) {
      cat("Error: obs is not a column vector.\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }   
    if(length(unique(tabulate(obs))) != 1) {
      cat("Error: Panel is not balanced [check obs vector].\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }
    Y.temp <- as.matrix(Y)
    if (any(Y.temp[,1] != Y)) {
      cat("Error: X matrix is not conformable [does not match Y].\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }
    if(dim(W)[1] != n * k) {
      cat("Error: W matrix is not conformable [does not match Y].\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }   
 
    # check iteration parameters
    check.mcmc.parameters(burnin, mcmc, thin)
    totiter <- mcmc + burnin

    # starting values for beta error checking
    beta.start <- NA
    ols.beta <- solve(t(X) %*% X) %*% t(X) %*% Y
    ols.sigma2 <-
      t(Y - X %*% ols.beta) %*% (Y - X %*% ols.beta) / (k*n - p - 1)
    ols.sigma2 <- as.double(ols.sigma2)
    if (is.na(beta.start)){ # use least squares estimates
      beta.start <- ols.beta
    }
    if(is.null(dim(beta.start))) {
      beta.start <- beta.start * matrix(1,p,1)  
    }
    if((dim(beta.start)[1] != p) || (dim(beta.start)[2] != 1)) {
      cat("Error: Starting value for beta not conformable.\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }
  
    # sigma2 starting values error checking
    if (is.na(sigma2.start)){
      sigma2.start <- ols.sigma2
    }   
    if(sigma2.start <= 0) {
      cat("Error: Starting value for sigma2 negative.\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }   

    # starting values for D error checking
    if (is.na(D.start)){ # use matrix of ones
      D.start <- .5 * ols.sigma2 * diag(q)
    }   
    if(is.null(dim(D.start))) {
      D.start <- D.start * diag(q)
    }
    if((dim(D.start)[1] != q) || (dim(D.start)[2] != q)) {
      cat("Error: Starting value for D not conformable.\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }

    # set up prior for beta
    if(is.null(dim(b0))) {
      b0 <- b0 * matrix(1,p,1)  
    }
    if((dim(b0)[1] != p) || (dim(b0)[2] != 1)) {
      cat("Error: N(b0,B0^-1) prior b0 not conformable.\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }  
    if(is.null(dim(B0))) {
      B0 <- B0 * diag(p)
    }
    if((dim(B0)[1] != p) || (dim(B0)[2] != p)) {
      cat("Error: N(b0,B0^-1) prior B0 not conformable.\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    } 
   
    # set up prior for sigma2
    if(nu0 <= 0) {
      cat("Error: G(nu0,delta0) prior nu0 less than or equal to zero.\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }
    if(delta0 <= 0) {
      cat("Error: G(nu0,delta0) prior delta0 less than or equal to zero.\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }
       
    # set up prior for D
    if(eta0 < q) {
      cat("Error: Wishart(eta0,R0) prior eta0 less than or equal to q.\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }   
    if(is.null(dim(R0))) {
      R0 <- R0 * diag(q)
    }
    if((dim(R0)[1] != q) || (dim(R0)[2] != q)) {
      cat("Error: Wishart(eta0,R0) prior R0 not comformable [q times q].\n")
      stop("Please respecify and call MCMCpanel() again.\n", call.=FALSE)
    }
         
    # set up big holder matrix
    sample <- matrix(0, mcmc/thin, p+q*q+1)
   
    # call C++ code to draw sample
    inv.obj <- .C("panelpost",
                  samdata = as.double(sample),
                  samrow = as.integer(nrow(sample)),
                  samcol = as.integer(ncol(sample)),
                  obsdata = as.double(obs),
                  obsrow = as.integer(length(obs)),
                  obscol = as.integer(1),   
                  ydata = as.double(Y),
                  yrow = as.integer(length(Y)),
                  ycol = as.integer(1),   
                  xdata = as.double(X),
                  xrow = as.integer(nrow(X)),
                  xcol = as.integer(ncol(X)),   
                  wdata = as.double(W),
                  wrow = as.integer(nrow(W)),
                  wcol = as.integer(ncol(W)),   
                  burnin = as.integer(burnin),
                  gibbs = as.integer(mcmc),
                  thin = as.integer(thin),
                  lecuyer = as.integer(lecuyer),
                  seedarray = as.integer(seed.array),
                  lecuyerstream = as.integer(lecuyer.stream),
                  verbose = as.integer(verbose),
                  bstartdata = as.double(beta.start),
                  bstartrow = as.integer(nrow(beta.start)),
                  bstartcol = as.integer(ncol(beta.start)),
                  sigma2start = as.double(sigma2.start),
                  Dstartdata = as.double(D.start),
                  Dstartrow = as.integer(nrow(D.start)),
                  Dstartcol = as.integer(ncol(D.start)),
                  b0data = as.double(b0),
                  b0row = as.integer(nrow(b0)),
                  b0col = as.integer(ncol(b0)),   
                  B0data = as.double(B0),
                  B0row = as.integer(nrow(B0)),
                  B0col = as.integer(ncol(B0)),  
                  nu0 = as.double(nu0),
                  delta0 = as.double(delta0),
                  eta0 = as.integer(eta0),    
                  R0data = as.double(R0),
                  R0row = as.integer(nrow(R0)),
                  R0col = as.integer(ncol(R0)),
                  n = as.integer(n),
                  k = as.integer(k),
                  p = as.integer(p),
                  q = as.integer(q),
                  PACKAGE="MCMCpack"
                  )     
    
    # put together matrix and build MCMC object to return
    sample <- matrix(inv.obj$samdata, inv.obj$samrow, inv.obj$samcol,
                     byrow=TRUE)   
  
    if (length(colnames(X))>0) {
        beta.names <- colnames(X)
    }
    else beta.names <- paste("beta", 1:p, sep = "")
    D.names <- paste("D", 1:(q*q), sep = "")
    sigma2.names <- "sigma2"
    names <- c(beta.names, D.names, sigma2.names)   

    output <- mcmc(data=sample, start=1, end=mcmc, thin=thin)
    varnames(output) <- names
    attr(output,"title") <- 
      "MCMCpack Linear Panel Model Posterior Density Sample"
    return(output)
  }

### sample from the posterior distribution of a Poisson regression
### model in R using linked C++ code in Scythe
###
### ADM 1/24/2003
## KQ 3/17/2003 [bug fix]
## Modified to meet new developer specification 7/15/2004 KQ
## Modified for new Scythe and rngs 7/26/2004 KQ

"MCMCpoisson" <-
  function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
           thin=1, tune=1.1, verbose = FALSE, seed = NA, beta.start = NA,
           b0 = 0, B0 = 0, ...) {
    
    ## checks
    check.offset(list(...))
    check.mcmc.parameters(burnin, mcmc, thin)
    
    ## seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]

    ## form response and model matrices
    holder <- parse.formula(formula, data)
    Y <- holder[[1]]
    X <- holder[[2]]
    xnames <- holder[[3]]    
    K <- ncol(X)  # number of covariates
        
    ## starting values and priors
    beta.start <- coef.start(beta.start, K, formula, family=poisson, data)
    mvn.prior <- form.mvn.prior(b0, B0, K)
    b0 <- mvn.prior[[1]]
    B0 <- mvn.prior[[2]]
    
    ## form the tuning parameter
    tune <- vector.tune(tune, K)
    V <- vcov(glm(formula=formula, data=data, family=poisson))
    
    ## test y non-negative
    if (sum(Y < 0) > 0) {
      cat("\n Elements of Y negative. ")
      stop("\n Check data and call MCMCpoisson() again. \n") 
    }
   
    ## define holder for posterior density sample
    sample <- matrix(data=0, mcmc/thin, dim(X)[2] )
  
    ## call C++ code to draw sample
    auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCpoisson",
                     sample=sample, Y=Y, X=X, burnin=as.integer(burnin),
                     mcmc=as.integer(mcmc), thin=as.integer(thin),
                     tune=tune, lecuyer=as.integer(lecuyer),
                     seedarray=as.integer(seed.array),
                     lecuyerstream=as.integer(lecuyer.stream),
                     verbose=as.integer(verbose), betastart=beta.start,
                     b0=b0, B0=B0, V=V) 
    
    ## put together matrix and build MCMC object to return
    output <- form.mcmc.object(posterior, names=xnames,
                               title="MCMCpoisson Posterior Density Sample")
    return(output)
    
  }


## sample from the posterior distribution of a probit
## model in R using linked C++ code in Scythe
##
## ADM and KQ 5/21/2002
## Modified to meet new developer specification 7/26/2004 KQ
## Modified for new Scythe and rngs 7/26/2004 KQ


"MCMCprobit" <-
  function(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
           thin = 1, verbose = FALSE, seed = NA, beta.start = NA,
           b0 = 0, B0 = 0, bayes.resid=FALSE, ...) {
    
    ## checks
    check.offset(list(...))
    check.mcmc.parameters(burnin, mcmc, thin)
    
    ## seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]

    ## form response and model matrices
    holder <- parse.formula(formula, data)
    Y <- holder[[1]]
    X <- holder[[2]]
    xnames <- holder[[3]]    
    K <- ncol(X)  # number of covariates
        
    ## starting values and priors
    beta.start <- coef.start(beta.start, K, formula,
                             family=binomial(link=probit), data)
    mvn.prior <- form.mvn.prior(b0, B0, K)
    b0 <- mvn.prior[[1]]
    B0 <- mvn.prior[[2]]

    ## residuals setup
    resvec <- NULL
    if (is.logical(bayes.resid) && bayes.resid==TRUE){
      resvec <- matrix(1:length(Y), length(Y), 1)
    }
    else if (!is.logical(bayes.resid)){
      resvec <- matrix(bayes.resid, length(bayes.resid), 1)
      if (min(resvec %in% 1:length(Y)) == 0){
        cat("Elements of bayes.resid are not valid row numbers.\n")
        stop("Check data and call MCMCprobit() again.\n") 
      }
    }
    
    ## y \in {0, 1} error checking
    if (sum(Y!=0 & Y!=1) > 0) {
      cat("Elements of Y equal to something other than 0 or 1.\n")
      stop("Check data and call MCMCprobit() again.\n") 
    }
     
    if (is.null(resvec)){
      ## define holder for posterior density sample
      sample <- matrix(data=0, mcmc/thin, dim(X)[2] )
  
      ## call C++ code to draw sample
      auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCprobit",
                       sample=sample, Y=Y, X=X, burnin=as.integer(burnin),
                       mcmc=as.integer(mcmc), thin=as.integer(thin),
                       lecuyer=as.integer(lecuyer),
                       seedarray=as.integer(seed.array),
                       lecuyerstream=as.integer(lecuyer.stream),
                       verbose=as.integer(verbose), betastart=beta.start,
                       b0=b0, B0=B0) 
      
      ## put together matrix and build MCMC object to return
      output <- form.mcmc.object(posterior, names=xnames,
                                 title="MCMCprobit Posterior Density Sample")

    }
    else{
      # define holder for posterior density sample
      sample <- matrix(data=0, mcmc/thin, dim(X)[2]+length(resvec) )

      ## call C++ code to draw sample
      auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCprobitres",
                       sample=sample, Y=Y, X=X, resvec=resvec,
                       burnin=as.integer(burnin),
                       mcmc=as.integer(mcmc), thin=as.integer(thin),
                       lecuyer=as.integer(lecuyer),
                       seedarray=as.integer(seed.array),
                       lecuyerstream=as.integer(lecuyer.stream),
                       verbose=as.integer(verbose), betastart=beta.start,
                       b0=b0, B0=B0) 
      
      ## put together matrix and build MCMC object to return
      xnames <- c(xnames, paste("epsilonstar", as.character(resvec), sep="") )
      
      output <- form.mcmc.object(posterior, names=xnames,
                                 title="MCMCprobit Posterior Density Sample")
      
    }
    return(output)
  
  }


# MCMCregress.R samples from the posterior distribution of a Gaussian
# linear regression model in R using linked C++ code in Scythe
#
# Original written by ADM and KQ 5/21/2002
# Updated with helper functions ADM 5/28/2004
# Modified to meet new developer specification 6/18/2004 KQ
# Modified for new Scythe and rngs 7/22/2004 ADM

"MCMCregress" <-
  function(formula, data=parent.frame(), burnin = 1000, mcmc = 10000,
   thin=1, verbose = FALSE, seed = NA, beta.start = NA,
   b0 = 0, B0 = 0, c0 = 0.001, d0 = 0.001, ...) {
    
    # checks
    check.offset(list(...))
    check.mcmc.parameters(burnin, mcmc, thin)
    
    # seeds
    seeds <- form.seeds(seed) 
    lecuyer <- seeds[[1]]
    seed.array <- seeds[[2]]
    lecuyer.stream <- seeds[[3]]

    # form response and model matrices
    holder <- parse.formula(formula, data)
    Y <- holder[[1]]
    X <- holder[[2]]
    xnames <- holder[[3]]    
    K <- ncol(X)  # number of covariates
        
    # starting values and priors
    beta.start <- coef.start(beta.start, K, formula, family=gaussian, data)
    mvn.prior <- form.mvn.prior(b0, B0, K)
    b0 <- mvn.prior[[1]]
    B0 <- mvn.prior[[2]]
    check.ig.prior(c0, d0)
   
    # define holder for posterior density sample
    sample <- matrix(data=0, mcmc/thin, K+1)

    # call C++ code to draw sample
    auto.Scythe.call(output.object="posterior", cc.fun.name="MCMCregress", 
                     sample=sample, Y=Y, X=X, burnin=as.integer(burnin),
                     mcmc=as.integer(mcmc), thin=as.integer(thin),
                     lecuyer=as.integer(lecuyer), 
                     seedarray=as.integer(seed.array),
                     lecuyerstream=as.integer(lecuyer.stream),
                     verbose=as.integer(verbose), betastart=beta.start,
                     b0=b0, B0=B0, c0=as.double(c0), d0=as.double(d0))
                  
    # pull together matrix and build MCMC object to return
    output <- form.mcmc.object(posterior,
                               names=c(xnames, "sigma2"),
                               title="MCMCregress Posterior Density Sample")
    return(output)
  }
## this function automates the Scythe C++ call making book-keeping
## much easier
##
##    output.object:  name of posterior sample that will be placed
##                    in the parent environment (string)
##
##    cc.fun.name:    name of the C++ function to be called (string)
##
##    package:        name of package (string, "MCMCpack" by default)
##
##    developer:      option that determines whether the R call and
##                    C++ template are echoed or whether the Scythe
##                    call is made (logical)
##
##    help.file:      option that determines whether a template of a
##                    helpfile for the calling R function should be
##                    generated (logical)
##
##    cc.file:        the file used to store the C++ template
##                    (string, output to the screen if "")
##
##    R.file:         the file used to store the clean R template
##                    (string, output to the screen if "")
##
##    ...:            list of objects passed to C++ 
##                    NOTE: this will only take integers (which have
##                    to be coerced), doubles, and matrices.  They
##                    should all be of the form "X = X," with the first
##                    part the C++ name and the second part the R name.
##                    Remember that C++ names cannot have periods in them.
##                    Matrices will be, for example, Xdata, Xrow, and
##                    Xcol.
##
##       This also build a skeleton C++ template and clean R template
##       for MCMCpack if developer=TRUE.  All are arguments constants
##       except the matrix named "sample".
##

"auto.Scythe.call" <-
  function(output.object, cc.fun.name, package="MCMCpack",
           developer=FALSE, help.file=FALSE, cc.file="", R.file="", ...) {

   # pull stuff from function call
   objects <- list(...)
   K <- length(objects)
   c.names <- names(objects)
   if (max(c.names=="sample") != 1){
     stop("argument named 'sample' must be specified in call to auto.Scythe.call()\n")
   }
   R.names <- strsplit(toString(match.call()), ",")[[1]]
   R.names <- R.names[(length(R.names)-K+1):length(R.names)]

   ## put default values in for burnin, mcmc, thin, and verbose
   ## if not explicitly supplied
   burnin.exist <- mcmc.exist <- thin.exist <- seed.exist <-
     verbose.exist <- FALSE
   for (k in 1:K){
     if(c.names[[k]]=="burnin" & is.integer(objects[[k]]))
       burnin.exist <- TRUE
     if(c.names[[k]]=="mcmc" & is.integer(objects[[k]]))
       mcmc.exist <- TRUE
     if(c.names[[k]]=="thin" & is.integer(objects[[k]]))
       thin.exist <- TRUE
     if(c.names[[k]]=="verbose" & is.null(dim(objects[[k]])))
       verbose.exist <- TRUE        
   }
   if (!burnin.exist){
     objects <- c(objects, burnin=as.integer(5000))
     R.names[length(objects)] <- "as.integer(5000)"
   }
   if (!mcmc.exist){
     objects <- c(objects, mcmc=as.integer(25000))
     R.names[length(objects)] <- "as.integer(25000)"     
   }
   if (!thin.exist){
     objects <- c(objects, thin=as.integer(1))
     R.names[length(objects)] <- "as.integer(1)"     
   }
   if (!verbose.exist){
     objects <- c(objects, verbose=as.integer(FALSE))
     R.names[length(objects)] <- "as.integer(FALSE)"     
   }
   K <- length(objects)
   c.names <- names(objects)

   ## check parameters
   check.mcmc.parameters(objects$burnin, objects$mcmc, objects$thin)
   
   ## write out a template R help file
   callfun <- strsplit(toString(sys.call(which=-1)),",")[[1]][1]
   if (help.file){
     prompt(callfun,file=paste(callfun, ".template.Rd", sep=""))
   }
   
   ##
   ## pull together R call
   ##

   # strings for R call
   start <- paste(".C('", cc.fun.name, "',", sep="")
   end <- paste("PACKAGE='", package, "')", sep="") 
   middle <- NULL
   
   for(k in 1:K) {
     if(is.double(objects[[k]]) & is.null(dim(objects[[k]]))) {
       if (regexpr("as.double", R.names[[k]])==-1){
         holder <- paste(c.names[[k]], " = as.double(", R.names[[k]], "),",
                         sep="")
         middle <- c(middle, holder)
       }
       else {
         holder <- paste(c.names[[k]], " =", R.names[[k]], ",",
                         sep="")
         middle <- c(middle, holder)
       }
     }
     else if(is.integer(objects[[k]]) & is.null(dim(objects[[k]]))) {
       if (regexpr("as.integer", R.names[[k]])==-1){
         holder <- paste(c.names[[k]], " = as.integer(", R.names[[k]], "),",
                         sep="")
         middle <- c(middle, holder)
       }
       else {
         holder <- paste(c.names[[k]], " =", R.names[[k]], ",",
                         sep="")
         middle <- c(middle, holder)
       }
     }
     else if(is.matrix(objects[[k]])) {
       holder.data <- paste(c.names[[k]], "data", " = as.double(",
                            R.names[[k]], "),", sep="")
       holder.row <- paste(c.names[[k]], "row", " =", " nrow(",
                           R.names[[k]], "),", sep="")
       holder.col <- paste(c.names[[k]], "col", " =", " ncol(",
                           R.names[[k]], "),", sep="")
       middle <- c(middle, holder.data, holder.row, holder.col)
     }
     else {
       stop("Integers, doubles, or matrices only to auto.Scythe.call().")
     }
   }


   
   # clean up and return R call
   middle <- paste(middle, sep=" ", collapse=" ")
   call <- paste(start, middle, end, sep=" ")
   call <- gsub('\\( ', '\\(', call)  

   ##   
   ## pull together C++ call
   ##

   # strings for C++ call
   c.start <- paste("void ", cc.fun.name, "(", sep="")
   c.end <- ")"   
   c.middle <- NULL
   together.call <- NULL

   for(k in 1:K) {
      if(is.double(objects[[k]]) & is.null(dim(objects[[k]]))) {
         holder <- paste("const double *", c.names[[k]], ",", sep="")
         c.middle <- c(c.middle, holder)
      }
      else if(is.integer(objects[[k]]) & is.null(dim(objects[[k]]))) {
         holder <- paste("const int *", c.names[[k]], ",", sep="")
         c.middle <- c(c.middle, holder)
      }
      else if(is.matrix(objects[[k]])) {
         
         # pull together Scythe call [note sample]
         if(c.names[[k]]=="sample") {
            holder.data <- paste("double *", c.names[[k]], "data,", sep="")
            scythe <- NULL
            together.call <- paste(together.call, scythe, sep="")
            sample.block <- paste("     const int size = *", c.names[[k]],
               "row * *", c.names[[k]],
               "col;\n     for(int i = 0; i < size; ++i)\n",
               "       ", c.names[[k]], "data[i] = STORAGEMATRIX[i];\n", 
               sep="")
         }
         else {
           holder.data <- paste("const double *", c.names[[k]],
                                "data,", sep="")
           scythe <- paste("     Matrix <double> ", c.names[[k]],
                           " = r2scythe(*",
                           c.names[[k]], "row, *", c.names[[k]], "col, ",
                           c.names[[k]], "data);\n", sep="")
                            
           together.call <- paste(together.call, scythe, sep="")
         }
         holder.row <- paste("const int *", c.names[[k]], "row,", sep="")
         holder.col <- paste("const int *", c.names[[k]], "col,", sep="")
         c.middle <- c(c.middle, holder.data, holder.row, holder.col)
      }
   }
   
   # clean up and print C++ function call
   c.middle <- paste(c.middle, sep=" ", collapse=" ")
   c.call <- paste(c.start, c.middle, c.end, sep="")
   c.call <- gsub(',)', ')', c.call)

   # if developer dump Scythe code to file, R function to screen, and evaluate
   if(developer) {
      comment.block <- paste("// ", cc.file, " DESCRIPTION HERE\n//\n// The initial version of this file was generated by the\n// auto.Scythe.call() function in the MCMCpack R package\n// written by:\n//\n// Andrew D. Martin\n// Dept. of Political Science\n// Washington University in St. Louis\n// admartin@wustl.edu\n//\n// Kevin M. Quinn\n// Dept. of Government\n// Harvard University\n// kevin_quinn@harvard.edu\n// \n// This software is distributed under the terms of the GNU GENERAL\n// PUBLIC LICENSE Version 2, June 1991.  See the package LICENSE\n// file for more information.\n//\n// Copyright (C) ", substring(date(),21,24), " Andrew D. Martin and Kevin M. Quinn\n// \n// This file was initially generated on ", date(), "\n// REVISION HISTORY\n\n", sep="")
      
      includes.block <- '#include "matrix.h"\n#include "distributions.h"\n#include "stat.h"\n#include "la.h"\n#include "ide.h"\n#include "smath.h"\n#include "MCMCrng.h"\n#include "MCMCfcds.h"\n\n#include <R.h>           // needed to use Rprintf()\n#include <R_ext/Utils.h> // needed to allow user interrupts\n\nusing namespace SCYTHE;\nusing namespace std;\n\n'
      
      main.block <- 'extern "C" {\n\n   // BRIEF FUNCTION DESCRIPTION\n'
      
      function.call <- paste('   ', c.call, ' {\n', sep="")
      
      together.block <- "   \n     // pull together Matrix objects\n     // REMEMBER TO ACCESS PASSED ints AND doubles PROPERLY\n"

      constants.block <- "\n     // define constants\n     const int tot_iter = *burnin + *mcmc;  // total number of mcmc iterations\n     const int nstore = *mcmc / *thin;      // number of draws to store\n     const int NUMBER_OF_PARAMETERS = ????; // YOU NEED TO FILL THIS IN\n"

      storage.block <-   "\n     // storage matrix or matrices\n     Matrix<double> STORAGEMATRIX(nstore, NUMBER_OF_PARAMETERS);\n"      

      seed.block <- "\n     // initialize rng stream\n     rng *stream = MCMCpack_get_rng(*lecuyer, seedarray, *lecuyerstream);\n"

      startval.block <- "\n     // set starting values\n     PARAMETER_BLOCK1 = ????;\n     PARAMETER_BLOCK2 = ????;\n     ETC.;\n"
      
      sample.call <- paste('\n     ///// MCMC SAMPLING OCCURS IN THIS FOR LOOP\n     for(int iter = 0; iter < tot_iter; ++iter){\n\n       // sample the parameters\n       PARAMETER_BLOCK1 = ????;\n       PARAMETER_BLOCK2 = ????;\n       ETC;\n\n       // store draws in storage matrix (or matrices)\n       if(iter >= *burnin && (iter % *thin == 0)){\n        // PUT DRAWS IN STORAGEMATRIX HERE\n       }\n\n       // print output to stdout\n       if(*verbose == 1 && iter % 500 == 0){\n         Rprintf("\\n\\n', cc.fun.name, ' iteration %i of %i \\n", (iter+1), tot_iter);\n         // ADD ADDITIONAL OUTPUT HERE IF DESIRED\n       }\n\n       void R_CheckUserInterrupt(void); // allow user interrupts\n     } // end MCMC loop\n\n     delete stream; // clean up random number stream\n\n     // load draws into sample array\n', sep="")
      
      end.block <- paste("\n   } // end", cc.fun.name,  '\n} // end extern "C"\n')

      if (cc.file == ""){
        cat("\n\nThe C++ template file is:\n")
      }
      cat(comment.block, includes.block, main.block, function.call,
          together.block, together.call, constants.block, storage.block,
          seed.block, startval.block, sample.call,
          sample.block, end.block, sep="", file=cc.file)
      if (cc.file != "") {
        cat("\nCreated file named '", cc.file, "'.\n", sep="")
        cat("Edit the file and move it to the appropriate directory.\n")
      }

      if (R.file == "") {
        cat("\n\nThe clean R template file is:\n")
      }
      dump(callfun, file=R.file)
      if (R.file != "") {
        cat("\nCreated file named '", R.file, "'.\n", sep="")
        cat("Edit the file and move it to the appropriate directory.\n")
        cat("Do not forget to edit the MCMCpack NAMESPACE file if\n")
        cat("installing new functions as part of MCMCpack.\n")
      }

      cat("\nThe call to .C is:\n")      
      draw.sample.call <- parse(text=paste(output.object, " <- ", call))
      print(draw.sample.call)
      cat("\nAUTOMATIC TEMPLATE FILE CREATION SUCCEEDED.\n\n")
      invokeRestart("abort")
    }

   # if not developer evaluate call leaving output.object in
   # parent environment
   if(!developer) {
      draw.sample.call <- parse(text=paste(output.object, " <- ", call))
      eval(draw.sample.call, envir=parent.frame(1))
   }
}

########## Density Functions and Random Number Generators ##########

##
## Wishart
##

# rwish delivers a pseudo-random Wishart deviate
#
# USAGE:
#
#   A <- rwish(v, S)
#
# INPUT:
#
#   v    degrees of freedom
#
#   S    Scale matrix
#
# OUTPUT:
#
#  A     a pseudo-random Wishart deviate
#
# Based on code originally posted by Bill Venables to S-news
# on 6/11/1998
#
# KQ on 2/5/2001

"rwish" <-
  function(v, S) {
    if (!is.matrix(S))
      S <- matrix(S)
    if (nrow(S) != ncol(S)) {
      stop(message="S not square in rwish().\n")
    }
    if (v < nrow(S)) {
      stop(message="v is less than the dimension of S in rwish().\n")
    }
    p <- nrow(S)
    CC <- chol(S) 
    Z <- matrix(0, p, p)
    diag(Z) <- sqrt(rchisq(p, v:(v-p+1)))
    if(p > 1) {
      pseq <- 1:(p-1)
      Z[rep(p*pseq, pseq) + unlist(lapply(pseq, seq))] <- rnorm(p*(p-1)/2)
    }
    return(crossprod(Z %*% CC))
  }

# dwish evaluations the Wishart pdf at positive definite matrix W.
# note: uses the Gelman, et. al. parameterization.
#
# USAGE:
#
#   x <- dwish(W, v, S)
#
# INPUT:
#
#   W    positive definite matrix at which to evaluate PDF
#
#   v    degrees of freedom
#
#   S    Scale matrix
#
# OUTPUT:
#
#   x    the PDF evaluated (scalar)
#
# ADM 8/16/2002

"dwish" <-
  function(W, v, S) {
    if (!is.matrix(S))
      S <- matrix(S)
    if (nrow(S) != ncol(S)){
      stop(message="W not square in dwish()\n\n")
    }
    if (!is.matrix(W))
      S <- matrix(W)
    if (nrow(W) != ncol(W)){
      stop(message="W not square in dwish()\n\n")
    }   
    if(nrow(S) != ncol(W)){
      stop(message="W and X of different dimensionality in dwish()\n\n")
    }
    if (v < nrow(S)){
      stop(message="v is less than the dimension of S in  dwish()\n\n")
    }    
    k <- nrow(S)
  
    # denominator
    gammapart <- 1
    for(i in 1:k) {
      gammapart <- gammapart * gamma((v + 1 - i)/2)
    } 
    denom <- gammapart *  2^(v * k / 2) * pi^(k*(k-1)/4)
  
    # numerator
    detS <- det(S)
    detW <- det(W)
    hold <- solve(S) %*% W
    tracehold <- sum(hold[row(hold) == col(hold)])  
    num <- detS^(-v/2) * detW^((v - k - 1)/2) * exp(-1/2 * tracehold)

    return(num / denom)
  }

##
## Inverse Wishart
##

# riwish generates a draw from the inverse Wishart distribution
# (using the Wishart generator)  

"riwish" <-
  function(v, S) {
    return(solve(rwish(v,S)))
  }

# diwish evaluates the inverse Wishart pdf at positive definite
# matrix W.  note: uses the Gelman, et. al. parameterization.
#
# USAGE:
#
#   x <- diwish(W, v, S)
#
# INPUT:
#
#   W    positive definite matrix at which to evaluate PDF
#
#   v    degrees of freedom
#
#   S    Scale matrix
#
# OUTPUT:
#
#   x    the PDF evaluated (scalar)
#
# ADM 8/16/2002
 
"diwish" <-
  function(W, v, S) {
    if (!is.matrix(S))
      S <- matrix(S)
    if (nrow(S) != ncol(S)){
      stop("W not square in diwish().\n")
    }
    if (!is.matrix(W))
      S <- matrix(W)
    if (nrow(W) != ncol(W)){
      stop("W not square in diwish().\n")
    }   
    if(nrow(S) != ncol(W)){
      stop("W and X of different dimensionality in diwish().\n")
    }
    if (v < nrow(S)){
      stop("v is less than the dimension of S in  diwish().\n")
    }
    
    k <- nrow(S)   

    # denominator
    gammapart <- 1
    for(i in 1:k) {
      gammapart <- gammapart * gamma((v + 1 - i)/2)
    } 
    denom <- gammapart *  2^(v * k / 2) * pi^(k*(k-1)/4)
  
    # numerator
    detS <- det(S)
    detW <- det(W)
    hold <- S %*% solve(W)
    tracehold <- sum(hold[row(hold) == col(hold)])  
    num <- detS^(v/2) * detW^(-(v + k + 1)/2) * exp(-1/2 * tracehold)

    return(num / denom)
  }

##
## Inverse Gamma
##

# evaluate the inverse gamma density
#
# Kevin Rompala 5/6/2003

"dinvgamma" <-
  function(x, shape, rate = 1) {

    # error checking
    if(shape <= 0 | rate <=0) {
      stop("Shape or rate parameter negative in dinvgamma().\n")
    }
   
    alpha <- shape
    beta <- rate
    return(beta^alpha / gamma(alpha) * x^(-1*(alpha + 1)) * exp(-beta/x))
  }

# generate draws from the inverse gamma density (using
# the gamma simulator)
#
# Kevin Rompala 5/6/2003

"rinvgamma" <-
  function(n, shape, rate = 1) {
    return(1 / rgamma(n, shape, rate))
  }

##
## Dirichlet (Multivariate Beta)
##

# ddirichlet evaluates the density of the Dirichlet at
# vector x given shape parameter vector (or matrix)
# alpha.
#
# note: this code is taken verbatim from the R-package
# "Greg's Miscellaneous Functions" maintained by
# Gregory R. Warnes <Gregory_R_Warnes@groton.pfizer.com>
#
# Kevin Rompala 5/6/2003

"ddirichlet" <-
  function(x, alpha) {

    dirichlet1 <- function(x, alpha) {
      logD <- sum(lgamma(alpha)) - lgamma(sum(alpha))
      s <- sum((alpha-1)*log(x))
      exp(sum(s)-logD)
    }

    # make sure x is a matrix
    if(!is.matrix(x))
      if(is.data.frame(x))
        x <- as.matrix(x)
      else
        x <- t(x)
    if(!is.matrix(alpha))
      alpha <- matrix( alpha, ncol=length(alpha), nrow=nrow(x), byrow=TRUE)

    if( any(dim(x) != dim(alpha)) )
      stop("Mismatch between dimensions of x and alpha in ddirichlet().\n")

    pd <- vector(length=nrow(x))
    for(i in 1:nrow(x))
      pd[i] <- dirichlet1(x[i,],alpha[i,])

    # Enforce 0 <= x[i,j] <= 1, sum(x[i,]) = 1
    pd[ apply( x, 1, function(z) any( z <0 | z > 1)) ] <- 0
    pd[ apply( x, 1, function(z) all.equal(sum( z ),1) !=TRUE) ] <- 0
    return(pd)
  }


# rdirichlet generates n random draws from the Dirichlet at
# vector x given shape parameter vector (or matrix)
# alpha.
#
# note: this code is taken verbatim from the R-package
# "Greg's Miscellaneous Functions" maintained by
# Gregory R. Warnes <Gregory_R_Warnes@groton.pfizer.com>
#
# Kevin Rompala 5/6/2003

"rdirichlet" <-
  function(n, alpha) {
    l<-length(alpha);
    x<-matrix(rgamma(l*n,alpha),ncol=l,byrow=TRUE);
    sm<-x%*%rep(1,l);
    return(x/as.vector(sm));
  }

##
## Non-Central Hypergeometric
##

# code to evaluate the noncentral hypergeometric density (at a single point
# or at all defined points).
#
# parameters:
#
#    n1, n2 -- number of subjects in group 1 and 2
#
#    Y1, Y2 -- number of subjects with positive outcome [unobserved]
#
#    psi -- odds ratio
#
#    m1 -- sum of observed values of Y1 and Y2 (Y1 + Y2)
#   
# output:
#
#   pi -- Pr(Y1 = x | Y1 + Y2 = m1) x=ll,...,uu
#
#   for ll = max(0, m1-n2) and uu = min(n1, m1)
#  
# if x is NA, then a matrix is returned, with the first column the possible
# values of x, and the second columns the probabilities.  if x is a scalar, 
# the density is evaluated at that point.
#
# ADM on 5/8/2003
#
# note: code adapted from R code published in conjunction with:
#
# Liao, J.G. And Rosen, O. (2001) Fast and Stable Algorithms for Computing and 
# Sampling from the Noncentral Hypergeometric Distribution.  The American
# Statistician 55, 366-369.
#

"dnoncenhypergeom" <-
  function (x = NA, n1, n2, m1, psi) {

    ##
    ## AUXILIARY FUNCTIONS
    ##

    mode.compute <- function(n1, n2, m1, psi, ll, uu) {
      a <- psi - 1
      b <- -( (m1+n1+2)*psi + n2-m1 )     
      c <- psi*(n1+1)*(m1+1)
      q <- b + sign(b)*sqrt(b*b-4*a*c)
      q <- -q/2
                         
      mode <- trunc(c/q) 
      if(uu>=mode && mode>=ll) return(mode)
      else return( trunc(q/a) )      
    }

    r.function <- function(n1, n2, m1, psi, i) {
      (n1-i+1)*(m1-i+1)/i/(n2-m1+i)*psi
    }

    ##
    ## MAIN FUNCTION
    ##

    # upper and lower limits for density evaluation
    ll <- max(0, m1-n2)
    uu <- min(n1, m1)

    # check parameters
    if(n1 < 0 | n2 < 0) {
       stop("n1 or n2 negative in dnoncenhypergeom().\n")
    }
    if(m1 < 0 | m1 > (n1 + n2)) {
       stop("m1 out of range in dnoncenhypergeom().\n")
    }
    if(psi <=0) {
       stop("psi [odds ratio] negative in dnoncenhypergeom().\n")
    }
    if(!is.na(x) & (x < ll | x > uu)) {
       stop("x out of bounds in dnoncenhypergeom().\n")
    }
    if(!is.na(x) & length(x) > 1) {
       stop("x neither missing or scalar in dnoncenhypergeom().\n")
    }  
  
    # evaluate density using recursion (from mode)
    mode <- mode.compute(n1, n2, m1, psi, ll, uu)
    pi <- array(1, uu-ll+1)
    shift <- 1-ll
    
    if(mode<uu) { # note the shift of location
      r1 <- r.function( n1, n2, m1, psi, (mode+1):uu )       
      pi[(mode+1 + shift):(uu + shift)] <- cumprod(r1)       
    }
   
    if(mode>ll) {
       r1 <- 1/r.function( n1, n2, m1, psi, mode:(ll+1) )
       pi[(mode-1 + shift):(ll + shift)] <- cumprod(r1)
    }
            
    pi <- pi/sum(pi)
    if(is.na(x)) return(cbind(ll:uu,pi))
    else return(pi[x + shift])
}

# code to generate random deviates from the noncentral hypergeometric density 
#
# parameters:
#
#    n -- the number of draws to make
#
#    n1, n2 -- number of subjects in group 1 and 2
#
#    Y1, Y2 -- number of subjects with positive outcome [unobserved]
#
#    psi -- odds ratio
#
#    m1 -- sum of observed values of Y1 and Y2 (Y1 + Y2)
#   
# output:
#
#   output -- a list of length n of random deviates
#  
#
# ADM on 5/9/2003
#
# note: code adapted from R code published in conjunction with:
#
# Liao, J.G. And Rosen, O. (2001) Fast and Stable Algorithms for Computing and 
# Sampling from the Noncentral Hypergeometric Distribution.  The American
# Statistician 55, 366-369.
#

"rnoncenhypergeom" <-
  function (n, n1, n2, m1, psi) {

    ##
    ## AUXILIARY FUNCTIONS
    ##  

    mode.compute <- function(n1, n2, m1, psi, ll, uu) {
      a <- psi - 1
      b <- -( (m1+n1+2)*psi + n2-m1 )     
      c <- psi*(n1+1)*(m1+1)
      q <- b + sign(b)*sqrt(b*b-4*a*c)
      q <- -q/2
                         
      mode <- trunc(c/q) 
      if(uu>=mode && mode>=ll) return(mode)
      else return( trunc(q/a) )
      
    }  
  
    sample.low.to.high <- function(lower.end, ran, pi, shift) { 
      for(i in lower.end:uu) {                                
        if(ran <= pi[i+shift]) return(i)
        ran <- ran - pi[i+shift]
        }                                
    }
       
    sample.high.to.low <- function(upper.end, ran, pi, shift) {
      for(i in upper.end:ll) {                              
        if(ran <= pi[i+shift]) return(i)
        ran <- ran - pi[i+shift]
      } 
    }
    
    single.draw <- function(n1, n2, m1, psi, ll, uu, mode, pi) {
      ran <- runif(1)
      shift <- 1-ll  
      if(mode==ll) return(sample.low.to.high(ll, ran, pi, shift))            
      if(mode==uu) return(sample.high.to.low(uu, ran, pi, shift))                                         
      if(ran < pi[mode+shift]) return(mode)             
      ran <- ran - pi[mode+shift]
      lower <- mode - 1                                                                            
      upper <- mode + 1
          
      repeat {           
        if(pi[upper + shift] >= pi[lower + shift]) {              
          if(ran < pi[upper+shift]) return(upper)
          ran <- ran - pi[upper+shift]
          if(upper == uu) return( sample.high.to.low(lower, ran, pi, shift) )
          upper <- upper + 1                            
        }
        else {
          if(ran < pi[lower+shift]) return(lower)
          ran <- ran - pi[lower+shift]
          if(lower == ll) return( sample.low.to.high(upper, ran, pi, shift) )
          lower <- lower - 1                   
        }     
      }
    }
  
    ##
    ## MAIN FUNCTION
    ##

    # upper and lower limits for density evaluation
    ll <- max(0, m1-n2)
    uu <- min(n1, m1)
    
    # check parameters
    if(n1 < 0 | n2 < 0) {
       stop("n1 or n2 negative in rnoncenhypergeom().\n")
    }
    if(m1 < 0 | m1 > (n1 + n2)) {
       stop("m1 out of range in rnoncenhypergeom().\n")
    }
    if(psi <=0) {
       stop("psi [odds ratio] negative in rnoncenhypergeom().\n")
    }


    # get density and other parameters
    mode <- mode.compute(n1, n2, m1, psi, ll, uu) 
    pi <- dnoncenhypergeom(NA, n1, n2, m1, psi)[,2]
    
    output <- array(0,n)
    for(i in 1:n) output[i] <- single.draw(n1, n2, m1, psi, ll, uu, mode, pi)    
    return(output)
  }
########## hidden functions to help in model implementation ##########

## NOTE: these are not exported to the user and should always be
##       used in model functions.  As such, fixing problems here
##       fixes them in all functions simultaneously.
##
## updated by ADM 7/22/04
## re-organized (alphabetical) by adm 7/28/04

## create an agreement score matrix from a vote matrix
## subjects initially on rows and items on cols of X
"agree.mat" <- function(X){
  X <- t(X) # put subjects on columns
  n <- ncol(X)
  X[is.na(X)] <- -999
  A <- matrix(NA, n, n)
  for (i in 1:n){
    A[i,] <- apply(X[,i] == X, 2, sum)
  }
  return(A/nrow(X))
}

## create constraints for measurement models
"build.factor.constraints" <-
  function(lambda.constraints, X, K, factors){

    ## build initial constraint matrices and assign var names
    Lambda.eq.constraints <- matrix(NA, K, factors)
    Lambda.ineq.constraints <- matrix(0, K, factors)    
    if (is.null(colnames(X))){
      X.names <- paste("V", 1:ncol(X), sep="")
    }
    else {
      X.names <- colnames(X)
    }    
    rownames(Lambda.eq.constraints) <- X.names
    rownames(Lambda.ineq.constraints) <- X.names

    ## setup the equality and inequality contraints on Lambda
    if (length(lambda.constraints) != 0){
      constraint.names <- names(lambda.constraints)  
      for (i in 1:length(constraint.names)){
        name.i <- constraint.names[i]
        lambda.constraints.i <- lambda.constraints[[i]]
        col.index <- lambda.constraints.i[[1]]
        replace.element <- lambda.constraints.i[[2]]
        if (is.numeric(replace.element)){
          Lambda.eq.constraints[rownames(Lambda.eq.constraints)==name.i,
                                col.index] <- replace.element
        }
        if (replace.element=="+"){
          Lambda.ineq.constraints[rownames(Lambda.ineq.constraints)==name.i,
                                  col.index] <- 1 
        }
        if (replace.element=="-"){
          Lambda.ineq.constraints[rownames(Lambda.ineq.constraints)==name.i,
                                  col.index] <- -1
        }
      }
    }
    
    testmat <- Lambda.ineq.constraints * Lambda.eq.constraints
    
    if (min(is.na(testmat))==0){
      if ( min(testmat[!is.na(testmat)]) < 0){
        cat("Constraints on factor loadings are logically inconsistent.\n")
        stop("Please respecify and call ", calling.function(), " again.\n")
      }
    }
    Lambda.eq.constraints[is.na(Lambda.eq.constraints)] <- -999
    
    return( list(Lambda.eq.constraints, Lambda.ineq.constraints, X.names))    
  }

# return name of the calling function
"calling.function" <-
   function(parentheses=TRUE) {
     calling.function <- strsplit(toString(sys.call(which=-3)),",")[[1]][1]
     if (parentheses){
       calling.function <- paste(calling.function, "()", sep="")
     }
     return(calling.function)
   }

# check inverse Gamma prior
"check.ig.prior" <-
   function(c0, d0) {
     
    if(c0 <= 0) {
      cat("Error: IG(c0/2,d0/2) prior c0 less than or equal to zero.\n")
      stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)
    }
    if(d0 <= 0) {
      cat("Error: IG(c0/2,d0/2) prior d0 less than or equal to zero.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
        call.=FALSE)    
    }
    return(0)
  }

# check mcmc parameters
"check.mcmc.parameters" <-
  function(burnin, mcmc, thin) {
  
    if(mcmc %% thin != 0) {
      cat("Error: MCMC iterations not evenly divisible by thinning interval.\n")
      stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
    }
    if(mcmc < 0) {
      cat("Error: MCMC iterations negative.\n")
      stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE) 
    }
    if(burnin < 0) {
      cat("Error: Burnin iterations negative.\n")
      stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
    }
    if(thin < 0) {
      cat("Error: Thinning interval negative.\n")
      stop("Please respecify and call ", calling.function(), " again.",
         call.=FALSE)
    }    
    return(0)
  }

# check to see if an offset is passed
"check.offset" <-
   function(args) {
      if(sum(names(args)=="offset")==1) {
         cat("Error: Offsets are currently not supported in MCMCpack.\n")
         stop("Please respecify and call ", calling.function(), " again.\n")
      }
   return(0)
   }

# put together starting values for coefficients
# NOTE: This can be used for any GLM model by passing the right family
#       or for another model by passing default starting values to
#       the function   
"coef.start" <-
   function(beta.start, K, formula, family, data, defaults=NA) {
          
     if (is.na(beta.start)[1] & is.na(defaults)[1]){ # use GLM estimates
       beta.start <- matrix(coef(glm(formula, family=family, data=data)), K, 1)
     }
     else if(is.na(beta.start)[1] & !is.na(defaults)[1]){ # use passed values
        beta.start <- matrix(defaults,K,1)     
     }
     else if(is.null(dim(beta.start))) {
       beta.start <- beta.start * matrix(1,K,1)  
       }
     else if(!all(dim(beta.start) == c(K,1))) {
       cat("Error: Starting value for beta not conformable.\n")
       stop("Please respecify and call ", calling.function(), " again.\n",
         call.=FALSE)
       }
     return(beta.start)
   }

## generate starting values for a factor loading matrix
"factload.start" <-
  function(lambda.start, K, factors, Lambda.eq.constraints,
           Lambda.ineq.constraints){

    Lambda <- matrix(0, K, factors)
    if (is.na(lambda.start)){# sets Lambda to equality constraints & 0s
      for (i in 1:K){
        for (j in 1:factors){
          if (Lambda.eq.constraints[i,j]==-999){
            if(Lambda.ineq.constraints[i,j]==0){
              Lambda[i,j] <- 0
            }
            if(Lambda.ineq.constraints[i,j]>0){
              Lambda[i,j] <- .5
            }
            if(Lambda.ineq.constraints[i,j]<0){
              Lambda[i,j] <- -.5
            }          
          }
          else Lambda[i,j] <- Lambda.eq.constraints[i,j]
        }
      }    
    }
    else if (is.matrix(lambda.start)){
      if (nrow(lambda.start)==K && ncol(lambda.start)==factors)
        Lambda  <- lambda.start
      else {
        cat("lambda.start not of correct size for model specification.\n")
        stop("Please respecify and call ", calling.function(), " again.\n")
      }
    }
    else if (length(lambda.start)==1 && is.numeric(lambda.start)){
      Lambda <- matrix(lambda.start, K, factors)
      for (i in 1:K){
        for (j in 1:factors){
          if (Lambda.eq.constraints[i,j] != -999)
            Lambda[i,j] <- Lambda.eq.constraints[i,j]
        }
      }    
    }
    else {
      cat("lambda.start neither NA, matrix, nor scalar.\n")
      stop("Please respecify and call ", calling.function, " again.\n")
    }

    return(Lambda)
  }




## based on code originally written by Keith Poole
## takes a subject by subject agreement score matrix as input
"factor.score.eigen.start" <- function(A, factors){

  AA <- A
  arow <- matrix(NA, nrow(A), 1)
  acol <- matrix(NA, ncol(A), 1)
  
  for (i in 1:nrow(A)){
    arow[i] <- mean(A[i,])
  }
  for (i in 1:ncol(A)){
    acol[i] <- mean(A[,i])
  }
  
  matrixmean <- mean(acol)
  
  for (i in 1:nrow(A)){
    for (j in 1:ncol(A)){
      AA[i,j] <- (A[i,j]-arow[i]-acol[j]+matrixmean)/(2)
    }
  }
  
  ev <- eigen(AA)
  scores <- matrix(NA, nrow(A), factors)
  for (i in 1:factors){
    scores[,i] <- ev$vec[,i]*sqrt(ev$val[i])
    scores[,i] <- (scores[,i] - mean(scores[,i]))/sd(scores[,i])
  }
  return(scores)
}


## check starting values of factor scores or ability parameters
## subjects on rows of X
"factor.score.start.check" <- function(theta.start, X, prior.mean,
                                       prior.prec, eq.constraints,
                                       ineq.constraints, factors){

  N <- nrow(X)
  
  ## set value of theta.start
  if (is.na(theta.start)) {
    theta.start <- factor.score.eigen.start(agree.mat(X), 1)
    for (i in 1:factors){
      theta.start[,i] <- prior.mean[i] + theta.start[,i] *
        sqrt(1/prior.prec[i,i])
    }
  }
  else if(is.numeric(theta.start) && is.null(dim(theta.start))) {
    theta.start <- theta.start * matrix(1, N, 1)  
  }
  else if((dim(theta.start)[1] != N) ||
           (dim(theta.start)[2] != factors)) {
    cat("Starting value for theta not appropriately sized.\n")
    stop("Please respecify and call", calling.function(), " again.\n",
         call.=FALSE)
  }
  else {
    cat("Inappropriate value of theta.start passed.\n")
    stop("Please respecify and call", calling.function(), " again.\n",
         call.=FALSE)      
  }

  ## check value of theta.start
  prev.bind.constraints <- rep(0, factors)
  for (i in 1:N){
    for (j in 1:factors){
      if (eq.constraints[i,j]==-999){
        if(ineq.constraints[i,j]>0 && theta.start[i,j] < 0){
          if (prev.bind.constraints[j]==0){
            theta.start[,j] <- -1*theta.start[,j]
          }
          else {
            cat("Parameter constraints logically inconsistent.\n")
            stop("Please respecify and call ", calling.function(), " again.",
                 call.=FALSE)                
          }
          prev.bind.constraints[j] <- prev.bind.constraints[j] + 1
        }
        if(ineq.constraints[i,j]<0 && theta.start[i,j] > 0){
          if (prev.bind.constraints[j]==0){
            theta.start[,j] <- -1*theta.start[,j]
          }
          else {
            cat("Parameter constraints logically inconsistent.\n")
            stop("Please respecify and call ", calling.function(), " again.",
                 call.=FALSE)                
          }
          prev.bind.constraints[j] <- prev.bind.constraints[j] + 1
        }
      }
      else {
        if ((theta.start[i,j] * eq.constraints[i,j]) > 0){
          theta.start[i,j] <- eq.constraints[i,j]
        }
        else {
          if (prev.bind.constraints[j]==0){
            theta.start[,j] <- -1*theta.start[,j]
            theta.start[i,j] <- eq.constraints[i,j]
          }
          else {
            cat("Parameter constraints logically inconsistent.\n")
            stop("Please respecify and call ", calling.function(), " again.",
                 call.=FALSE)                
          }
          prev.bind.constraints[j] <- prev.bind.constraints[j] + 1
        }
      }      
    }
  }    
  return(theta.start)
}

## get starting values for factor uniqueness matrix (Psi)
"factuniqueness.start" <-
  function(psi.start, X){
    
    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=TRUE)

       output <- mcmc(data=holder, start=1,
                      end=posterior.object$mcmc,
                      thin=posterior.object$thin)
       varnames(output) <- names
       attr(output,"title") <- title
       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))) {
       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)
     }
     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)
  } 
  return(diag(as.double(mcmc.tune)))
}

########## Scythe Inter-Operation Functions ##########


# writes a matrix out to an ASCII file that can be read by Scythe.
# it puts the number of rows and columns in the first row
# followed by the data.
#
# ADM 1/29/2003

"write.Scythe" <-
  function(outmatrix, outfile = NA, overwrite=FALSE) {
    outmatrix <- as.matrix(outmatrix)
   
    if(is.na(outfile)) {
      stop("Please specify a file name in the write.Scythe() call.\n")
    }
    if(overwrite==FALSE & file.exists(outfile)) {
      cat("File already exists in the write.Scythe() call.\n")
      stop("Either delete the file, or flip the overwrite switch.\n")
    }
   
    outfile <- file(outfile, "w")
    cat(dim(outmatrix), "\n", file=outfile)
    write.table(outmatrix, file=outfile,
                row.names=FALSE, col.names=FALSE, quote=FALSE)
    close(outfile)
    return(0)
  } 


# reads in a matrix from an ASCII file written by Scythe.
# the number of rows and columns should be in the first row followed
# by the data.
#
# Kevin Rompala 5/1/2003
# fixed by ADM 7/25/2004

"read.Scythe" <-
  function(infile = NA) {
    
    if(is.na(infile)) {
      stop("Please specify a file name in the read.Scythe() call.\n")
    }
    if(!file.exists(infile)) {
      stop("Specified source file does not exist in read.Scythe() call.\n")
    }

    infile <- file(infile, "r")
    dimensions <- scan(file=infile,n=2)
    inputdata <- scan(file=infile)
    close(infile)
    hold <- matrix(data=inputdata,
           nrow=dimensions[1], ncol=dimensions[2], byrow=TRUE)
    return(hold) 
  }
########## Tomography Plots for Ecological Inference ##########

# produces tomography plots (see King, 1997, A Solution to the
# Ecological Inference Problem, Princeton University Press)
#
# KQ 11/9/2002

"tomogplot" <-
  function(r0, r1, c0, c1, 
           xlab="fraction of r0 in c0 (p0)",
           ylab="fraction of r1 in c0 (p1)",
           bgcol="white", ...) {
    if (length(r0) != length(r1)) {
      stop("r0 and r1 different lengths in tomogplot().\n")
    }
    if (length(r0) != length(c0)) {
      stop("r0 and c0 different lengths in tomogplot().\n")
    }
    if (length(r0) != length(c1)) {
      stop("r0 and c1 different lengths in tomogplot().\n")
    }
    
    intercept <-  c0/r1
    slope <- -1 * r0/r1
    N <- length(r0)
    
    par(pty="s")
    plot(0:1, 0:1, type="n", main="", xlab=xlab, ylab=ylab)
    rect(0, 0, 1, 1, col=bgcol, lty=0)
    
    for (year in 1:N) {
      abline(intercept[year], slope[year])
    }
    
    rect(-0.05, -0.05, 1.05, 0, col="white", lty=0)
    rect(-0.05, -0.05, 0, 1.05, col="white", lty=0)
    rect(-0.05, 1, 1.05, 1.05, col="white", lty=0)
    rect(1, -0.05, 1.05, 1.05, col="white", lty=0)
    box()
    par(pty="m")
    return(0)
  }

# produces temporally organized tomography plots
# (see King, 1997, A Solution to the Ecological Inference
# Problem, Princeton University Press)
#
# KQ 11/9/2002

"dtomogplot" <-
  function(r0, r1, c0, c1, time.vec=NA, delay=0, 
           xlab="fraction of r0 in c0 (p0)",
           ylab="fraction of r1 in c0 (p1)",
           color.palette=heat.colors,
           bgcol="black", ...) {
    if (length(r0) != length(r1)){
      stop("r0 and r1 different lengths in dtomogplot().\n")
    }
    if (length(r0) != length(c0)){
      stop("r0 and c0 different lengths in dtomogplot().\n")
    }
    if (length(r0) != length(c1)){
      stop("r0 and c1 different lengths in dtomogplot().\n")
    }
    if (length(r0) != length(time.vec) & !is.na(time.vec)[1]){
      stop("r0 and time.vec different lengths in dtomogplot().\n")
    }
  
    intercept <-  c0/r1
    slope     <- -1 * r0/r1
    N <- length(r0)
    if (is.na(time.vec)[1])
      time.vec <- 1:N
    col.vec <- color.palette(N)

    mar.orig <- (par.orig <- par(c("mar", "las", "mfrow")))$mar
    on.exit(par(par.orig))
    w <- (3 + mar.orig[2]) * par("csi") * 2.54
    layout(matrix(c(2,1), nc=2), widths=c(1,lcm(w)))
    par(las=1)
    mar <- mar.orig
    mar[4] <- mar[2]
    mar[2] <- 1
    par(mar=mar)
    par(pty="m")
    plot.new()
    plot.window(xlim=c(0,1), ylim=range(time.vec), xaxs="i",
                yaxs="i")
    rect(0, time.vec[-length(time.vec)], 1, time.vec[-1], col=col.vec)
    axis(4)
    box()
    mar <- mar.orig
    mar[4] <- 1
    par(mar=mar)
    par(pty="s")   
    plot(0:1, 0:1, type="n", main="", xlab=xlab, ylab=ylab)
    rect(0, 0, 1, 1, col=bgcol, lty=0)
  
    for (year in 1:N) {
      time.last <- proc.time()[3]
      time.next <- proc.time()[3]
      while ( (time.next - time.last) < delay ){
        time.next <- proc.time()[3]        
      }
      abline(intercept[year], slope[year], col=col.vec[year])
    }

    rect(-0.05, -0.05, 1.05, 0, col="white", lty=0)
    rect(-0.05, -0.05, 0, 1.05, col="white", lty=0)
    rect(-0.05, 1, 1.05, 1.05, col="white", lty=0)
    rect(1, -0.05, 1.05, 1.05, col="white", lty=0)
    box()  
    par(pty="m")
    return(0)
  }
########## Utility Functions ##########

# takes a symmetric matrix x and returns lower diagonal
# note: does not check for symmetry
#
# ADM 4/18/2003 

"vech" <-
  function (x) {
    x <- as.matrix(x)
    if (dim(x)[1] != dim(x)[2]) {
      stop("Non-square matrix passed to vech().\n")
    }
    output <- x[lower.tri(x, diag = TRUE)]
    dim(output) <- NULL
    return(output)
  }

# takes vector x and returns an nrow times nrow symmetric matrix
# this will recycle the elements of x as needed to fill the matrix
#
# ADM 4/18/2003
# ADM 11/13/2003 [bug fix]

"xpnd" <-
  function (x, nrow) {
    dim(x) <- NULL
    output <- matrix(0, nrow, nrow)
    output[lower.tri(output, diag = TRUE)] <- x
    hold <- output
    hold[upper.tri(hold, diag=TRUE)] <- 0
    output <- output + t(hold)    
    return(output)
  }
.onAttach <- function(...) {
   cat("##\n## Markov Chain Monte Carlo Package (MCMCpack)\n")
   cat("## Copyright (C) 2003, 2004, Andrew D. Martin and Kevin M. Quinn\n")
   cat("##\n## Support provided by the U.S. National Science Foundation\n")
   cat("## (Grants SES-0350646 and SES-0350613)\n##\n")
   require(coda, quietly=TRUE)
   require(MASS, quietly=TRUE)
}

.onUnload <- function(libpath) {
    library.dynam.unload("MCMCpack", libpath)
}
