new.urn <- function(n.balls=2,black.proportion=0.5) {
  return(c('black'=n.balls*black.proportion,
           'white'=n.balls*(1-black.proportion)))
}

proportion <- function(urn) {
  return(urn['black']/sum(urn))
}

draw.ball <- function(urn) {
  prob.black <- proportion(urn)
  drawn.ball <- c('white','black')[1*(runif(1)<prob.black)+1]
  urn[drawn.ball] <- urn[drawn.ball] + 1
  return(urn)
}

# Initial proportion of black balls
theta <- 0.5
n.balls.to.start <- 10
n.draws <- 100
precision <- 10000

draws <- numeric(precision)
for ( i in 1:precision ) {
  urn <- new.urn(n.balls.to.start,theta)
  for ( j in 1:n.draws ) {
    urn <- draw.ball(urn)
  }
  draws[i] <- proportion(urn)
}

observed.test.statistic <- 0.60
p.value <- sum(draws>=observed.test.statistic) / precision
ci <- p.value + c(-1,1)*qnorm(1-0.05/2)*sqrt(p.value*(1-p.value)/precision)
cat("p-value is ",p.value,", with a 95% confidence interval of (",ci[1],",",ci[2],")\n",sep="")
#table(draws)

hist(draws)

# What is the rejection region?
# How could you do a power calculation?