bsc <- function (counts, surveyed, order = T, R = 10000) {
  library(boot)
  if (length(counts) != length(surveyed)) {stop
    ("Vectors must be of equal length.")
  }
  binom <- unlist(mapply(c,
                         lapply(counts, function(x) rep(1, x)),
                         lapply(surveyed - counts, function(x) rep(0, x)),SIMPLIFY=F)
  )
  snumb <- paste("s", 1:length(surveyed), sep="")
  ident <- rep(snumb, surveyed)
  
  scsample <- data.frame(binom,ident)
  
  if (order == TRUE) {n <- factorial(length(counts))}
  if (order == FALSE) {n <- 1}
  
  bsc <- function(x, id) {
    sc1 <- tapply(x[id,1], x[id,2], mean)
    sc <- prod(sc1) / n
    return(sc)
  }
  
  boot.scsample <- boot(scsample, statistic = bsc, R, strata = scsample[, 2]
  )
  return(boot.scsample)
}


bsc2 <- function(bsc.alt1a, bsc.alt2a, order = T, R = 10000){
  library(boot)
  bsc.alt1 <- bsc.alt1a$data
  bsc.alt2 <- bsc.alt2a$data
  bsc.alt1$scid <- "first"
  bsc.alt2$scid <- "second"
  bsc.diff.df <- rbind(bsc.alt1,bsc.alt2)
  bsc.diff.df$comb <- as.factor(paste(bsc.diff.df$scid,bsc.diff.df$ident, sep = ""))
  
  bsc.diff.df$scid <- NULL
  bsc.diff.df$ident <- NULL
  
  if (order == TRUE) {  n1 <- factorial(length(unique(bsc.alt1$ident)))
  n2 <- factorial(length(unique(bsc.alt2$ident)))}
  if (order == FALSE) { n1 <- 1
  n2 <- 1}
  
  l <- length(unique(bsc.alt1$ident))
  m <- length(unique(bsc.alt2$ident))
  
  bsc.diff <- function(x, id) {
    sc1 <- tapply(x[id,1], x[id,2], mean)
    sca <- (prod(sc1[1:l]) / n1)
    scb <- (prod(sc1[(l+1):(l+m)]) / n2)
    sc <- sca - scb
    return(sc)
  }
  
  boot.diff <- boot(bsc.diff.df, statistic = bsc.diff, R, strata = bsc.diff.df[, 2]
  )
  return(boot.diff)
}


summary.bsc <- function (bsc.alt) {
  bsc.ci.alt <- boot.ci(bsc.alt, type="bca")
  title <- "BOOTSTRAPPING SOUND CHANGES"
  prob <- paste("Estimated P =", round(bsc.alt$t0*100, digits = 5), "%")
  bca <- paste("Estimated 95 % BCa CI = [", round(bsc.ci.alt$bca[4]*100, digits = 4),"%,",
               round(bsc.ci.alt$bca[5]*100, digits = 4),"%]")
  #rnsc <- paste(pasteR, n.sc.paste, countsp, surveyed, sep = "\n")
  probbca <- paste(prob, bca, sep = "\n")
  cat(title, probbca, sep = "\n\n")
}


summary.bsc2 <- function (bsc2.alt) {
  bsc2.ci.alt <- boot.ci(bsc2.alt, type="bca")
  title <- "BOOTSTRAPPING SOUND CHANGES - COMPARE"
  prob <- paste("Estimated",expression(Delta), "P =", round(bsc2.alt$t0*100, digits = 5), "%")
  bca <- paste("Estimated 95 % BCa CI = [", round(bsc2.ci.alt$bca[4]*100, digits = 4),"%,",
               round(bsc2.ci.alt$bca[5]*100, digits = 4),"%]")
  if (bsc2.ci.alt$bca[4] > 0 & bsc2.ci.alt$bca[5] > 0) {
    sig <- "P(A1) is significantly higher than P(A2)."
  }
  else if (bsc2.ci.alt$bca[4] < 0 & bsc2.ci.alt$bca[5] < 0) {
    sig <- "P(A1) is significantly lower than P(A2)."
  } else {
    sig <- "P(A1) and P(A2) are not significantly different."
  }
  probbca <- paste(prob, bca, sep = "\n")
  cat(title, probbca,sig, sep = "\n\n")
}