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

dnd.BF=function(s1,s2,
  alpha=NULL,priorinfluence=1,
  nbin=min(length(s1),length(s2))-1,bins=NULL,
  nsamp1=1e4,nsamp=5e3, 
  fixedsample=F,
  ci=.95,
  fasttest=F,sig=.95,
  allconverge=F,prec=.01,
  maxN=3e4,
  verbose=F,addattr=F) {
  # Bayes factor test of samples s1 and s2 from distributions Y and Z 
  # respectively, with three choices Y>Z, Y<Z and Y<>Z
 
  # PARAMETERS
  # s1, s2: vector samples of values from distributions Y and Z being tested
  #
  # PRIOR
  # alpha: dirichlet prior parameter
  # priorinfluence: if alpha not specified set so that prior has
  #  influence equal to this number of samples
  #
  # DATA HISTOGRAM
  # nbin: number of used by test, default is a fine division
  # bins: cutpoints for multinomial bins, if NULL set to nbin-1
  #   equally spaced quantiles of the joint (s1 U s2) sample  
  #
  # MC SAMPLE SIZE PER CYCLE
  # nsamp1: number of prior or posterior samples for first pass
  # nsamp: numer of prior or posteriOR sampleS in later passes
  #
  # HOW MANY SAMPLES?
  # fixedsample=T: ignore convergence tests and take a sample >= maxN
  # maxN: minimum number of fixed smaples or maximum number of samples 
  #   to take before giving up on convergence
  #
  # TESTS FOR WHEN TO STOP SAMPLING
  # ci: 100*ci% credible interval width used in convergence tests 
  #
  # DEFAULT TEST
  # fasttest=T pull out when no ci's cross value "sig"
  # sig: used by fasttest, this posterior model probability not within a ci   
  # STRICTER TEST
  # allconverge=T (ignored if fasttest=T): pull out when no ci widths 
  #   greater than prec
  # prec: used by allconverge, required minimum precision, if allconverge=F
  #  then prec is minimum precision of largest posterior model probability
  # 
  # OUTPUT CONTROL
  # verbose: print counter values and ci's on each cycle of sampling 
  # addattr: if true an attribute "mc" added to BF element of output list
  
  # VALUE
  # a list with two three named elements:
  #  choice: name of choice with highest BF
  #    Y>Z  (s1 dominates s2)
  #    Y<Z  (s2 dominates s1)
  #    Y<>Z (non-dominance)
  # ln.BF and BFp: vector of ln(Bayes factors) and posterior model 
  #   probabilities named as in "choice"
  # If addattr BFp given attribute "mc", a list containing:
  #   npost:  number of posterior samples of three types
  #   nprior: number of prior samples of three types
  #   nsamp: number of samples from prior and from posterior
  #   cum1: cumulative number of samples in each bin of Y data
  #   cum2: cumulative number of samples in each bin of Z data
  #   alpha: dirichlet prior parameter
  #   bins: bin boundaries used
 
  isdom=function(p1,p2) {
  # number of dominant samples 1>2, 2>1 and 
  # remainder (2<>1, non-dominance)  
    db12 <- p1[,1]>p2[,1]
    db21 <- p1[,1]<p2[,1]
    P1 <- p1[,1]
    P2 <- p2[,1]
    nbin1 <- dim(p1)[2]-1
    i=2
    repeat {
      if (i==nbin1 || (!any(c(db12,db21)))) break
      p1[db12,1] <- p1[db12,1]+p1[db12,i]
      p2[db12,1] <- p2[db12,1]+p2[db12,i]
      db12[db12] <- p1[db12,1]>p2[db12,1]
      P1[db21] <-P1[db21]+p1[db21,i]
      P2[db21] <- P2[db21]+p2[db21,i]
      db21[db21] <- P1[db21]<P2[db21]
      i=i+1
    }
    out <- c(sum(db12),sum(db21))
    c(out,dim(p1)[1]-sum(out))
  }

  tabnodrop=function(samp,bins) {
  # tabulates without dropping levels
    nout <- length(bins)-1
    out <- rep(0,nout)
    names(out) <- 1:nout
    tab <- table(cut(samp,bins,labels=F))
    out[names(tab)] <- as.vector(tab)
    out
  }
 
  cis=function(ci,npos,nprs,N) {
  # upper and lower 100*ci% credible interval values 
  # for posterior model probabilities 

    cip=function(ci,npo,npr,N) {
    # One posterior model probability credible interval

      calcci=function(ci,n,N) {
      # credible interval for a proportion with Beta(1,1) prior
        c(qbeta((1-ci)/2,n+1,N-n+1),qbeta(1-(1-ci)/2,n+1,N-n+1))
      }

      popci <- calcci(ci,npo,N)
      prpci <- calcci(ci,npr,N)
      c(popci[1]/(prpci[2]+popci[1]),popci[2]/(prpci[1]+popci[2]))
    }
 
    out <- vector(length=0)
    for (i in 1:length(npos)) 
      out <- cbind(out,cip(ci,npos[i],nprs[i],N))
    out
  }
  
  rdirichlet=function (n, alpha) {
  # modified from MCMCpack to never return a zero 
  # (sets zeros to arond min double precision = 3e-324
    l <- length(alpha)
    x <- matrix(rgamma(l * n, alpha), ncol = l, byrow = TRUE)
    sm <- x %*% rep(1, l)
    out=x/as.vector(sm)
    out[out<3e-324] <- 3e-324
    return(out)
  }
  
  
  # MAIN BODY
  
  # get histograms
  if (is.null(bins)) 
    bins <- unique(c(-Inf,quantile(c(s1,s2),
      seq(1/nbin,(nbin-1)/nbin,by=1/nbin)),Inf))
  n1 <- tabnodrop(s1,bins)       # s1 samples
  n2 <- tabnodrop(s2,bins)       # s2 samples
  
  # dirichlet prior
  if (is.null(alpha)) 
    alpha <- rep(priorinfluence/nbin,nbin)

  # zero counters
  nposts=0   # number of posterior samples
  npriors=0  # number of prior samples
  N=0        # number of samples
  ns=nsamp1  # number of samples for current pass
  converge=F
  
  repeat {
  # repeat sampling until estimates sufficiently precise
    priorp1 <- rdirichlet(ns,alpha)
    priorp2 <- rdirichlet(ns,alpha)
    postp1 <- rdirichlet(ns,n1+alpha)
    postp2 <- rdirichlet(ns,n2+alpha)
    npost <- isdom(postp1,postp2)
    nprior <- isdom(priorp1,priorp2)
    nposts <- nposts+npost
    npriors <- npriors+nprior
    N <- N+ns
    if (!fixedsample) { # Test for convergence
      dcis <- cis(ci,nposts,npriors,N) # credible interval on proportions
      if (fasttest) {
        converge <- !any(dcis[1,]<sig & dcis[2,]>sig)
      } else {
        if (allconverge) 
          converge <- all(apply(dcis,2,diff)<prec) else { 
            BF <- (nposts+1)/(npriors+1)
            converge <- apply(dcis,2,diff)[BF==max(BF)]<prec
          }
      }
    }
    if (verbose) {
      cat(paste("Sample counts (out of ",N,")\n",sep=""))
      names(npriors)=c("prior: Y>Z","Z>Y","Z<>Y")
      print(npriors)
      names(nposts)=c("post:  Y>Z","Z>Y","Z<>Y")
      print(nposts)
      if (!fixedsample) {
      	cat(paste("Posterior model probability ",round(100*ci),
      	  "% credible interval\n",sep=""))
      	dimnames(dcis)=list(c("Lower","Upper"),c("Y>Z","Z>Y","Z<>Y"))
      	print(round(dcis,4))
      }
      cat("\n")
    }
    if (converge | N>=maxN) break
    ns <- nsamp
  }
  BF <- (nposts+1)/(npriors+1)     # Counting Bayes factor
  names(BF) <- c("Y>Z","Z>Y","Z<>Y")
  BFp <- BF/sum(BF) # posterior model probabilities
  if (addattr) attr(BFp,"mc") <- 
    list(npost=npost,nprior=nprior,
         cum1=cumsum(n1),cum2=cumsum(n2),
         nsamp=N,alpha=alpha,bins=bins)
  list(choice=names(BF)[BF==max(BF)][1],ln.BF=log(BF),p=BFp)
}