wald.test <- function(x,n, p=0.5, pdiff=0, conf.level=0.95, alternative="two.sided") {
  ### Wald Test function
  ### 

#
  if( (length(x)==2) && (length(n)==1) )  n=c(n,n)
  dn = paste(deparse(substitute(x)),"out of", deparse(substitute(n)))


  #
  if(length(x)==1) {   ## One-Sample Wald Test
  	x=x[1];n=n[1]
    if( x<5 || n-x<5 ) warning("The sample size is not large enough (see Section 8.4)")
    test = "One-Sample Wald Test (Hawkes only!!)"
    nv   = "proportions"
    phat=x/n
    se = sqrt( phat*(1-phat)/n)
    
    sew = sqrt( p*(1-p)/n)
    ts = (phat-p)/sew
    c = 1-conf.level
    if(alternative=="two.sided") { 
      pval = 2*pnorm(-abs(ts))
      alternativeHypothesis = paste("true proportion is not equal to",p)
      c=(1-conf.level)/2
    } else if(alternative=="less") {
      pval = pnorm(ts)
      alternativeHypothesis = paste("true proportion is less than",p)
    } else if (alternative=="greater") {
      pval = 1-pnorm(ts)
      alternativeHypothesis = paste("true proportion is greater than",p)
    } else {
      stop("Alternative hypothesis is not one of the three required:\n\t\"less\", \"two.sided\", or \"greater\"")
    }
    
    E = -qnorm(c)*se
    lcl = phat - E
    ucl = phat + E
    
    est=phat
    pest=phat
    names(est)[1] <- paste( "proportion of",deparse(substitute(x1)) )
    expectedValue=phat
  }
  
  
  
  if(length(x)==2) { ## Two-sample Wald tets
  	x1=x[1];x2=x[2];n1=n[1];n2=n[2]
  	
    if( x1<5 || n1-x1<5 || x2<5 || n2-x2<5 ) warning("At least one sample size is not large enough (see Section 9.4)")
    test  = "Two-Sample Wald Test (Hawkes only!!)"
    nv    = "difference in proportions"
    phat1 = x1/n1; phat2=x2/n2
    se = sqrt(phat1*(1-phat1)/n1 + phat2*(1-phat2)/n2)
    
    p = (x1+x2)/(n1+n2)
    ts = ((phat1-phat2)-pdiff) / sqrt(p*(1-p)*(1/n1+1/n2))
    
    c=1-conf.level
    if(alternative=="two.sided") { 
      pval = 2*pnorm(-abs(ts))
      alternativeHypothesis = paste("true difference in proportions is not equal to",pdiff)
      c = c/2
    } else if(alternative=="less") {
      pval = pnorm(ts)
      alternativeHypothesis = paste("true difference in proportions less than",pdiff)
    } else if (alternative=="greater") {
      pval = 1-pnorm(ts)
      alternativeHypothesis = paste("true difference in proportions is greater than",pdiff)
    } else { 
      stop("Alternative hypothesis is not one of the three required:\n\t\"less\", \"two.sided\", or \"greater\"")
    }
    
    E = -qnorm(c) * se
    lcl = phat1-phat2 - E
    ucl = phat1-phat2 + E
    est=c(phat1,phat2)
    pest=phat1-phat2
    names(est)[1] <- "prop 1"
    names(est)[2] <- "prop 2"
    expectedValue = pdiff
  }
#



### Set up the return list
#
# Populate the results object
results <- list(
  statistic=ts, 
  parameter= Inf,
  p.value=pval,
  conf.int=c(lcl, ucl),
  estimate=est,
  null.value=expectedValue,
  alternative=alternativeHypothesis,
  method=test,
  data.name=dn
  )
#
attr(results$conf.int,"conf.level") <- conf.level
names(results$statistic) = "z"
names(results$parameter) = "df"
names(results$null.value) = nv
#
class(results) <- "htest"
return(results)

} # End wald.test function







dBenford <- function(d=1:9) { 
d = round(d)
## Check
 if(any(d<1) | any(d>9)) stop("The leading digit must be \n\tan integer between 1 and 9.")
log10( 1+1/(d) ) 
}


dGBenford <- function(turnout) {
 theta = log10(turnout)
 dist = numeric()

 for(d in 1:9) {
  pt0 = floor(theta)/theta * log10(1+1/d)
  pt1 = 0
  pt2 = (theta-floor(theta)+log10(1/d))/theta
  pt3 = log10(1+1/d)/theta

  add = pt2
  InitDigit = 10^(theta-floor(theta))

  if( InitDigit<d )     add = pt1
  if( d+1<=InitDigit )  add = pt3

  dist[d] = pt0+add
 }

dist
}


## The Benford Test
BenfordTest <- function(counts) {

# Determine the distribution of the initial digits
theDist = getInitDigitDistribution(counts)

# Perform the Chi-squared test
chisq.test(theDist, p=log10(1+1/(1:9)))

}



getInitDigit <- function(count) {
  as.numeric(substring(count,1,1))
}

getInitDigitDistribution <- function(count) {
ld = getInitDigit(count)
tabulate(ld,9)
}




benfordTest <- function(counts) {
 obs = getInitDigitDistribution(counts)
 exp = dBenford()
chisq.test(obs,p=exp)
}




findMultinomialProbabilities1 = function(turnout) {
 avgTurnout = mean(turnout)
 dGBenford(avgTurnout)
}


findMultinomialProbabilities2 = function(turnout) {
  n = length(turnout)
  p = matrix(NA, nrow=n,ncol=9)
  for(ed in 1:n) p[ed, ]=dGBenford(turnout[ed])
apply(p,2,mean)
}






getElectionLikelihood <- function(votes,turnout) {
 dig1s = getInitDigit(votes)
#theta = log10(turnout)
 nDivs = length(turnout)

 logprob=0
 for( division in 1:nDivs) {
  dgb  = dGBenford(turnout[division])
  vote = dig1s[division]
  logprob = logprob + log10(dgb[vote])
 }
logprob
} ##











### The Generalized Benford Test with 3 options

GBenfordTest <- function(votes,turnout, test="MA1",B=1e4 ) 
{

orig = options("warn")
options("warn"=-1)

###
if(test=="LL") {
message("Note:\tThis method uses simulation. As such, \n\tit may be rather slow. Please be aware of this.")
flush.console()
}

###
if(test=="MA1") {
  testName = "Generalized Benford Test using Multinomial Averaging, Type 1\n  Averaging the turnouts"
  dist = findMultinomialProbabilities1(turnout)
}

###
if(test=="MA2") {
  testName = "Generalized Benford Test using Multinomial Averaging, Type 2\n  Averaging the distributions"
  dist = findMultinomialProbabilities2(turnout)
}

if(test!="LL") {
 obs = getInitDigitDistribution(votes)
 tt=chisq.test(obs,p=dist)

cat("\n\n",testName,"\n\n", sep="")
cat("Estimated p-value:\n",tt$p.value,"\n\n")
if( min(sum(obs)*dist)<5 ) {
  cat("\n\nNote:\tBecause of the small number of divisions, the \n\tChi-squared approximation may be incorrect.\n\n")
}
flush.console()

res = list( p.value=tt$p.value)


}



### 

if(test=="LL") {

obs = getElectionLikelihood(votes,turnout)

# theta = log10(turnout)
 nDivs = length(turnout)

 lnlik = numeric()

# Estimate distribution
 for(elxn in 1:B) {

  # an election likelihood
  logprob=0

  for(division in 1:nDivs) {
    dgb  = dGBenford(turnout[division])
    vote = sample(1:9, size=1, prob=dgb)
    logprob = logprob + log10(dgb[vote])
  }
 lnlik[elxn] = logprob
 }
# Dist estimated


## pval
pval = min( mean(lnlik<=obs),mean(lnlik>=obs) )*2

testName = "Generalized Benford Test using Likelihood Estimation"

cat("\n\n",testName,"\n\n", sep="")
cat("Observed log-likelihood:\n",obs,"\n")
cat("Estimated p-value:\n",pval,"\n\n")
flush.console()

res = list( p.value=pval, obs.loglikelihood=obs, generated.loglikelihood=lnlik)
}


options("warn"=orig[[1]])

invisible(res)

}