Here is a heuristic argument for the mean of the order statistics $s_1,s_2,\ldots,s_n$.
Let's define a random variable $X_{ij}$ as follows
$$X_{ij}=\begin{cases} 1 & \text{, if player $i$ beats player $j$,}\\ 0 & \text{, if player $i$ loses to player $j$.}\end{cases}$$
Then, if we assume all the players are equally skilled, the $X_{ij}\sim \text{Ber}(0.5)$, i.e. they are Bernoulli random variables with probability of success $0.5$. However, they are not independent. Indeed, we can more precisely state that
$$\begin{eqnarray}X_{ij}\sim\text{Ber}(0.5) & \;\; , & \;\; \text{ for } i<j \\ X_{ij}=1-X_{ji} & \;\; , &\;\; \text{ for } i>j \\ X_{ii}=0 & \;\; , &\;\; \text{ players don't play against themselves } \\ \end{eqnarray} $$
This defines an $n\times n$ matrix whose components are the individual match results. The rowsums corresponds to the total amount of wins, the colomn sums to the total amount of losses.
$$W_i=\sum_{j=1}^n X_{ij} \;\; \text{ and } \;\; L_j=\sum_{i=1}^n X_{ij}$$
Of course, the losses and wins of a player count up to the total amount of matches played
$$W_i+L_i=n-1 \; .$$
We are interested in the order statistics $s_1=W_{(1)},s_2=W_{(2)},\ldots,s_n=W_{(n)}$. $W_{(i)}$ is the i-th element of the ascendingly sorted vector of $(W_1,W_2,\ldots,W_n)$.
Solving this problem is a pretty hard nut to crack. So, I resorted to a heuristic argument that only covers the expectation of the order statistics. What is easy to find is the marginal distribution of a player $i$. These will all be the same and binomial: $W_i\sim \text{Bin}(n-1,0.5)$.
Now, I reasoned that since an individual player could rank anywhere in the hierarchy, the possible outcomes of $W_i$ encode precisely that. Moreover, one could interpret the quantiles of the distribution as follows: if player $i$ ends up to be the player that wins most matches, this means that he is somewhere in the $1/n$th upper quantile of the distribution. If he loses the most matches, he'll be in the $1/n$th lower quantile. Etc, for other quantiles. So, I reasoned that essentially, the function for the means of the $s_i$ is just the quantile distribution of an individual player's wins.
So I did a small simulation in R and found that this indeed gives good results. I simulated 10000 roundrobin tournaments of 100 players. (Don't mind the probably inefficient code.)
n<-100;
mtres<-rep(0,n)
nsim<-10000
for (t in 1:nsim) {
tourres<-matrix(rep(0,n*n),nrow=n,ncol=n)
for (k in 2:n) {
for (l in 1:(k-1)) {
tourres[k,l]<-sample(c(0,1),1)
tourres[l,k]<-1-tourres[k,l]
}
}
mtres<-mtres+sort(rowSums(tourres))
}
Then I plotted the mean of the $s_1,s_2,\ldots,s_n$ against the quantiles of the binomial (I added a shift of 0.5 continuity correction) and also the normal approximation of the binomial.
plot(mtres/nsim,type='l')
me<-(n-1)/2
se<-sqrt((n-1)/4)
curve(qnorm(x/n,me,se),from=0,to=n,add=T,col='red')
curve(qbinom((x-0.5)/n,n-1,1/2),from=0,to=n,add=T,col='blue')
This is the result:
You almost don't see the difference between the normal approximation and the simulated result. (Although in the tails the difference will be larger)
Here's also the comparison between the simulated results and the binomial quantiles:
> mtres/nsim
[1] 37.0529 38.8178 39.8103 40.5246 41.1004 41.5688 41.9787 42.3477 42.6726
[10] 42.9847 43.2639 43.5249 43.7751 44.0082 44.2256 44.4402 44.6479 44.8379
[19] 45.0221 45.2040 45.3828 45.5611 45.7286 45.8976 46.0504 46.2097 46.3602
[28] 46.5153 46.6663 46.8137 46.9510 47.0826 47.2206 47.3658 47.5138 47.6548
[37] 47.7893 47.9121 48.0304 48.1495 48.2822 48.4198 48.5623 48.6912 48.8136
[46] 48.9286 49.0407 49.1549 49.2819 49.4216 49.5650 49.7027 49.8218 49.9305
[55] 50.0414 50.1550 50.2857 50.4243 50.5592 50.7022 50.8349 50.9580 51.0779
[64] 51.1992 51.3318 51.4758 51.6228 51.7688 51.9080 52.0362 52.1782 52.3212
[73] 52.4726 52.6264 52.7842 52.9417 53.0956 53.2592 53.4346 53.6118 53.7920
[82] 53.9683 54.1521 54.3493 54.5514 54.7600 54.9891 55.2233 55.4744 55.7365
[91] 56.0231 56.3181 56.6480 57.0165 57.4239 57.9007 58.4657 59.1880 60.1857
[100] 61.9515
> qbinom(((1:n)-0.5)/n,n-1,1/2)
[1] 37 39 40 41 41 42 42 42 43 43 43 44 44 44 44 44 45 45 45 45 45 46 46 46 46
[26] 46 46 47 47 47 47 47 47 47 48 48 48 48 48 48 48 48 49 49 49 49 49 49 49 49
[51] 50 50 50 50 50 50 50 50 51 51 51 51 51 51 51 51 52 52 52 52 52 52 52 53 53
[76] 53 53 53 53 54 54 54 54 54 55 55 55 55 55 56 56 56 57 57 57 58 58 59 60 62
>
If you round the simulation results to the nearest integer, the difference is zero.
A more able statistician than me can probably make this rigorous and add error margins to this result. But I thought it was already neat enough to post this answer.
EDIT : I wanted to add some additional findings. For small values of $n$, the above approximations will be abysmally bad. But, in those cases, it is not to hard to compute by hand the distributions of $(s_1,s_2,\ldots,s_n)$.
Take the case $n=3$ for instance. There's really only two situations possible, being $(1,1,1)$ and $(0,1,2)$. By looking at the matrices I defined before, it's not difficult to see that you can realize the first situation in $1/4$ cases and the second in $3/4$ cases. From there you can compute means, variance, all you want.
Take the case $n=4$. Then you have 4 different possibilities: $(0,1,2,3)$ with probability $4!/2^6$, $(1,1,1,3)$ with probability $8/2^6$, $(0,2,2,2)$ with probability $8/2^6$ and finally $(1,1,2,2)$ with probability $4!/2^6$. You can check that all the realization are really partitions of the number $6=4(4-1)/2$ with the restriction that no element of the partition is larger than 3, plus additional restrictions. The mean values are easily computed to be $(0.5,1.125,1.875,2.5)$. And indeed the simulation agrees with this result.
I wouldn't be surprised if the problem can be recast into the theory of random partitions.
EDIT2 I wrote a little program that computes the exact probability distribution. But it's not very efficient memory or time wise. I had done the computation by hand for a 5 player tournament, and then I tried for a 6 player but decided it would be wiser to program the rules. The result is this
n<-6
m<-n*(n-1)/2
store<-c()
for (k in 0:(2^m-1)) {
x<-as.numeric(intToBits(k))[1:m]
y<-1-x
vec1<-rep(0,n)
for (ell in (n-1):1) {
vec1[n-ell]<-sum(x[1:ell])
x<-x[(ell+1):length(x)]
}
vec2<-rep(0,n)
for (ell in (n-1):1) {
temp<-seq(1,1+ell*(n-ell-1),ell)
vec2[n-ell+1]<-sum(y[temp])
y<-y[-temp]
}
store<-c(store,toString(sort(vec1+vec2)))
}
table(store)
table(store)/2^m
barplot(table(store)/2^m)
I tried to run it for $n=7$ (which is about 64 times more computations and more storage space needed than for $n=6$) but I did not have the patience to wait for the end result. So, Monte Carlo simulation is more efficient to determine round robin tournament results. Or my heuristic if you look at tournaments with very high numbers of participants.
Best Answer
The MathOverflow post has several alternative links in the comments, and in particular Berge's proof (Theorem 6 in section 10.2 of Graphs and Hypergraphs) seems fairly nice.
First, Berge proves a theorem about complements. (Here, if $G$ is a directed graph, $\overline{G}$ is the directed graph which has all possible directed edges $G$ does not have; if $G$ has no edge in either direction between $u$ and $v$, then $\overline G$ has an edge in both directions.) I will write $h(G)$ for the number of Hamiltonian paths in $G$.
Let $s(H)$ denote the number of ways to number the vertices of $H$ as $1$ through $n$ such that all edges are of the form $(i,i+1)$. Then $s(H)=0$ unless $H$ is a union of paths; if $H$ has $k$ path components, then $s(H) = k!$. Therefore $s(H) \equiv 1 \pmod 2$ only if $H$ is a directed path. In particular, $$h(G) \equiv \sum_{E(H) \subseteq E(G)} s(H) \pmod 2$$ where the sum is over all spanning subgraphs of $G$.
Meanwhile, if we take a Hamiltonian path in $\overline G$ and number the vertices of $G$ along that path, then no edge of $G$ is of the form $(i,i+1)$, because all such edges are in $\overline G$. So $n! - h(\overline G)$ counts the number of labelings in which at least one edge of $G$ is of the form $(i,i+1)$. By the principle of inclusion-exclusion, we get $$ n! - h(\overline G) = \sum_{E(H) \subseteq E(G)} (-1)^{|E(H)|} s(H) \equiv \sum_{E(H) \subseteq E(G)} s(H) \pmod 2 $$ and as we saw earlier, this means $n! - h(\overline G) \equiv h(G) \pmod 2$. For $n>1$, this gives $h(G) \equiv h(\overline G) \pmod 2$.
Now we are ready to prove:
The idea is to check that the parity of $h$ does not change if we reverse an edge. Thus, the parity is the same as the number of Hamiltonian paths in a transitive tournament, which is $1$.
Let $e$ be an edge of $T$; let $f$ be the reverse of $T$. We define four directed graphs from $T$: $G_0 = T - e$, $G_e = T = G_0 + e$, $G_f = G_0 + f$, and $G_{ef} = G_0 + e + f$. So these agree everywhere except about which of edges $e$ and $f$ they have.
Note that $\overline{G_{ef}}$ gives us $G_0$ with all of its edges reversed, and reversing all the edges does not change the number of Hamiltonian paths. So by Theorem 1, $h(G_{ef}) \equiv h(G_0) \pmod 2$.
We have $h(G_{ef}) - h(G_0) = (h(G_e) - h(G_0)) + (h(G_f) - h(G_0))$. Here, the LHS counts Hamiltonian paths that use edge $e$ or $f$, and the RHS has two terms counting each of those cases separately. Therefore $$ h(G_e) + h(G_f) = h(G_{ef}) + h(G_0) \equiv 0 \pmod 2 $$ from which we get $h(G_e) \equiv h(G_f) \pmod 2$.