# (c) 2004, Jonathan Taylor, Stanford University

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

mvnorm.model <- function(cov.3=cov.TOTAL.3, cov.2=cov.TOTAL.2) {
  poiss.mean.YES <- numeric(ncenter)
  poiss.mean.NO <- numeric(ncenter)
  for (j in 1:ncenter) {
    mean <- mu.TOTAL[j]
    count <- nmachine[j]
    p <- mu.YES[j] / mean
    if (count > 1) {
      poiss.mean.YES[j] <- mvnorm.model.prob.tie(mean, count, p, cov.2=cov.2, cov.3=cov.3)
      poiss.mean.NO[j] <- mvnorm.model.prob.tie(mean, count, 1-p, cov.2=cov.2, cov.3=cov.3)
    }
  }
  mean.YES <- sum(poiss.mean.YES)
  mean.NO <- sum(poiss.mean.NO)
  return(data.frame(mean.YES, mean.NO))
}

print(mvnorm.model())

# A fairly close model to the one reported in The Economist

min.var = min(diag(cov.TOTAL.3))
print(mvnorm.model(cov.3=diag(rep(min.var,3)), cov.2=diag(rep(min.var,2))))