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

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)
}

proportion <- function(urn) {
  return(urn['black']/sum(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
print(p.value)
table(draws)