Beginner question on dice roll statistical test I have 2 cohorts of 10,000 people. I give each person in cohort 1 a die $a$, and ask them to roll it. I record the values as a tally. I then repeat this process with cohort 2 and give them die $b$.
I now have a table that looks like this:
cohort 1: 1=500, 2=1500, 3=5000, 4=2000, 5=1000
cohort 2: 1=100, 2=1000, 3=5000, 4=2500, 5=1400
I want to try and work out whether these results are statistically different. More specifically I want to understand if one die is more likely to produce a higher number.
While I can think of methods that do a lot of multiple comparisons such as $\chi^2$ tests I can't help think there must be a much better way to do this.
 A: Whats wrong with the chi-squared test? Just look at the $t = \sum \frac{(E-O)^2}{E}$ on $df = (5-1)(2-1)$.
Doing the calculations
countA = c(500,1500,5000,2000,1000)
countB = c(100,1000,5000,2500,1400)
t = sum((countA-countB)^2/countA)
1-pchisq(t,df = 4)

Doing this you get a p-value <<<< 0.05 So you know that the 2 are different. From here if you want to find the one most likely to give a higher number (over many throws) look for the one with the highest expect value.
resA = c(rep(1,500), rep(2,1500), rep(3,5000), rep(4,2000), rep(5,1000))
resB = c(rep(1,100), rep(2,1000), rep(3,5000), rep(4,2500), rep(5,1400))
mean(resA)
[1] 3.15
mean(resB)
[1] 3.41

So knowing they are statistically different and knowing that die b has a higher average you can say that die b is more likely to give a higher number over repeated trials.
You could also use the frequencies to calculate the probabilities of one dice throwing higher than the other.
$$P(B>A) = \sum_{N=2}^5(B=N~and~ A<(N-1))$$
Since the throws are independent 
$$P(B>A) = \sum_{N=2}^5[P(B=N) P(A<(N-1))]$$
$$P(B>A) = \sum_{N=2}^5[P(B=N) \sum_{n=1}^NP(A=n)]$$
Here we get $P(B>A) = 0.406$ and $P(A>B) = 0.2645$ Here we can see that dice B has a higher chance of having a larger roll than A.
A: Suppose one die $X$ can produce the values $x_1 \lt x_2 \lt \cdots \lt x_m$ with probabilities $p_1,p_2,\ldots, p_m$ (respectively) and the other die $Y$ can produce values $y_1\lt y_2 \lt \cdots \lt y_n$ with probabilities $q_1, q_2, \ldots, q_n$ (respectively).  Consider a game in which $X$ and $Y$ are independently tossed.  For independent throws, the chance $X$ exceeds $Y$ is
$$\Pr(X\gt Y) = \sum_{i=1}^m \sum_{j\mid x_i \gt y_j}^n p_iq_j.$$
From your (many) observations you can estimate the chances $p_i$ and $q_j$ and thereby estimate this probability.  
To assess the uncertainty in the probability estimate, consider bootstrapping the data.  One bootstrap iteration consists of constructing dice $X^{*}$ and $Y^{*}$ having the observed probabilities and rolling the pair ten thousand times--that is, for as many times as you have observations--and tallying the frequencies of the events $X^{*} \gt Y^{*},$ $X^{*} \lt Y^{*},$ or $X^{*}=Y^{*}.$  After several bootstrap iterations--50 may do, 500 is good, and an efficient algorithm with a modern computer will permit many more--you can study the distributions of these frequencies.  They reveal the likely uncertainties in the inferences you will be drawing from your data.
Allow me to illustrate with a smaller dataset having more uncertainty.  I have divided your counts by 100, producing two sets of 100 observations.  Here are the results of 10,000 bootstrap iterations:
 
For example, the lefthand plot headed "$P(x \gt y)$" shows the two bootstrapped dice tended to tie about 26% of the time, but the rates of ties varied from 15% to 40%.  Proceeding to the right, we see $X^{*}$ was less than $Y^{*}$ about 41% of the time ($X^{*}$ "lost" and $Y^{*}$ "won"), but this frequency varied from 25% to 55%; and $X^{*}$ won about 33% of the time, varying from 25% to 44%.  Finally, at the far right is a histogram of the difference in winning frequencies.  This difference averaged about -0.15 but extended from -0.4 to +0.1 within the bootstrap.  The vertical red line is the threshold of zero: only 3% of all the results are positive.  Put another way, in only 3% of the bootstrap iterations did $X^{*}$ win more often than $Y^{*}.$
The 3% number is a bootstrap "p-value" for testing whether $X$ has a lower rate of winning than $Y.$  As such, it is "one sided."  At the outset you likely wanted to test only whether $X$ and $Y$ have different winning rates.  That p-value for that "two sided" test is obtained by doubling the one-sided p-value to 6%.
The two-sided bootstrap p-value of 6% is small enough to give some evidence that $X$ will tend to win less often than $Y.$  It often would not be considered "significant" evidence where only values of 5% or less are considered sufficient evidence.
With the large amount data you actually have, the results are clear: on average $X^{*}$ wins 26.5% of the time, $Y^{*}$ wins 40.5% of the time, (and the remaining times they tie).  The p-value is less than $10^{-4},$ which is so tiny we may infer with high confidence that $Y$ wins more often than $X:$ about 14% more often.  (Notice the nature of this inference: from observed behaviors of $X^{*}$ and $Y^{*}$ we are making deductions about $X$ and $Y.$  That's the essential idea of the bootstrap.)  In this particular bootstrap calculation (again with 10,000 iterations) the difference in winning frequencies ranges from -17% to -11%, reflecting some remaining uncertainty about how much inferior $X$ is to $Y:$ but inferior it is.
Finally, note that the right hand histogram provides information not afforded by, say, a chi-squared test: its quantitative display of plausible alternative estimates shows you the amount of risk you undertake when making an inference about winning rates from your data.

The following R code performed the bootstrapping and (after first dividing x and y by 100) created the figure.
x <- c(500, 1500, 5000, 2000, 1000); names(x) <- 1:5
y <- c(100, 1000, 5000, 2500, 1400); names(y) <- 1:5
#
# Compute chances of winning, losing, and tieing given two dice.
#
p <- function(x, y) {
  p <- outer(x, y) / (sum(x) * sum(y))
  x.gt.y <- outer(as.numeric(names(x)), as.numeric(names(y)), `>`)
  x.lt.y <- outer(as.numeric(names(x)), as.numeric(names(y)), `<`)
  x.eq.y <- ! (x.gt.y | x.lt.y)
  c(`P(x > y)`=sum(p[x.gt.y]), `P(x < y)`=sum(p[x.lt.y]), `P(x==y)`=sum(p[x.eq.y]))
}
#
# Bootstrap the data.
#
n <- 1e4
x.boot <- rmultinom(n, sum(x), x) # This summarizes n * sum(x) = 10^8 throws of `X`
y.boot <- rmultinom(n, sum(y), y) # This summarizes n * sum(y) = 10^8 throws of `Y`
bootstrap <- sapply(seq.int(n), function(i) p(x.boot[,i], y.boot[,i])) # Compare the dice
bootstrap <- rbind(bootstrap, 
                   `P(x>y)-P(x<y)`=bootstrap["P(x > y)", ] - bootstrap["P(x < y)", ])
#
# Plot the results.
#
par(mfrow=c(1,4))
sapply(rownames(bootstrap), function(s) hist(bootstrap[s,], main=s, xlab="Value"))
abline(v=0, col="Red", lwd=2)
par(mfrow=c(1,1))
#
# Compute and report a one-sided p-value.
#
(p.value <- mean(bootstrap["P(x>y)-P(x<y)",] > 0))

