# (c) 2004, Jonathan Taylor, Stanford University

source('http://www-stat.stanford.edu/~jtaylo/venezuela/load-data.R')

total.YES <- round(mu.YES * nmachine)
total.NO <- round(mu.NO * nmachine)

multinomial.model <- function(count.YES=total.YES, count.NO=total.NO, p.assign=rep(1/3,3), num.machines=nmachine) {
  prob.YES.multinom <- numeric(ncenter)
  prob.NO.multinom <- numeric(ncenter)
  sample.center <- sample(1:ncenter)
  for (j in 1:ncenter) {
    prob.YES.multinom[j] <- multinomial.model.prob.tie(count.YES[j], num.machines[j], p.3=p.assign)
    prob.NO.multinom[j] <- multinomial.model.prob.tie(count.NO[j], num.machines[j], p.3=p.assign)
  }

  mean.YES <- sum(prob.YES.multinom)
  mean.NO <- sum(prob.NO.multinom)

  return(data.frame(mean.YES, mean.NO))
}

print(multinomial.model(p.assign=rep(1/3,3)))
print(multinomial.model(p.assign=c(0.3,0.3,0.4)))