TrueSkill: The probability of a predicted or observed draw This has me somewhat stumped as even the TrueSkill authors put it into the too-hard basket in a sense. I refer, of course, to:
R. Herbrich, T. Minka, and T. Graepel, “TrueSkill(TM): A Bayesian Skill Rating System,” Microsoft Res., Jan. 2007, Accessed: May 10, 2018. [Online]. Available: https://www.microsoft.com/en-us/research/publication/trueskilltm-a-bayesian-skill-rating-system/.
where on page 2 the first footnote reads:

The transitive relation "1 draws with 2" is not modelled exactly by
the relation $\left\lvert t_1 - t_2\right\rvert \le \epsilon$, which
is non-transitive. If $\left\lvert t_1 - t_2\right\rvert \le \epsilon$
and $\left\lvert t_2 - t_3\right\rvert \le \epsilon$ the the model
generates a draw among three teems despite the possibility that
$\left\lvert t_1 - t_3\right\rvert \le \epsilon$.

I admit I don't understand what they mean by transitive and non-transitive here, though I understand the point they are making and the meaning the difference between transitive and intransitive verbs. And clarification of the jargon here would please me, but is not the problem I'm stuck on.
To clarify, the extract above refers to $t_i$ which is the performance of team $i$ in an observed of predicted ranking which is modelled as a random variable with a Gaussian distribution.
If we consider performance thusly modelled:
$$
P_i \sim \mathcal{N}(p_i \mid \mu_{pi},\sigma_{pi}^2)
$$
and consider the relationship of two performers $P_1$ and P_2$ we can write:
$$
P_\delta \sim \mathcal{N}(p_\delta \mid \mu_{p\delta},\sigma_{p\delta}^2)
$$
where:
$$
\begin{align*}
P_\delta &= P_1 - P_2\\
\mu_{\delta} &= \mu_{1}-\mu_{2}\\
\sigma_{\delta}^2 &= \sigma_{1}^2+\sigma_{2}^2
\end{align*}
$$
$P_\delta$ as the difference between two Gaussian distributions is itself a Gaussian distribution as described, and describes the difference in performance between performer 1 and performer 2.
If we define a draw margin $\epsilon$ such that if $P_\delta \le \epsilon$ defines a draw and $P_\delta > \epsilon$ defines a victory to $P_1$ (remember $P_1 > P_2$ in the ranking), then it's not hard to show (from the Normal CDF $\Phi$) that:
$$
\begin{align*}
\Pr(\mathrm{win}) &= \Pr(P_i-P_{i+1}>\epsilon) = \Phi\left(\frac{\mu_\delta-\epsilon}{\sigma_\delta}\right) \\
\Pr(\mathrm{draw}) &= \Pr(\left\lvert P_i-P_{i+1}\right\rvert\le\epsilon) = \Phi\left(\frac {\epsilon-\mu_\delta} {\sigma_\delta}\right) - \Phi\left(\frac {-\epsilon-\mu_\delta} {\sigma_\delta}\right)\\\\
\end{align*}
$$
The thing that has me stumped is precisely what I cited from the paper above, how to calculate the probability of a 3 way draw (or more generally n-way). These simple conclusions are comparing only two performances.
How can I represent and calculate say the probability that $P_1$, $P_2$ and $P_3$ would all draw or did draw.
An obvious formulation is to walk the line, calculate the it as the joint porbability of each pair drawing that is formulate it as:
$$
\Pr(\mathrm{draw}) = \Pr(P_1-P_2\le\epsilon) \Pr(P_2-P_3\le\epsilon)
$$
This immediately confronts the problem the TrueSkill authors lived with and admitted to that in this draw it is possible that $P_1-P_3>\epsilon$.
Of course that might be fixed naively (as I am sure this wrong) as:
$$
\Pr(\mathrm{draw}) = \Pr(P_1-P_2\le\epsilon) \Pr(P_2-P_3\le\epsilon) \Pr(P_1-P_3\le\epsilon)
$$
Now I am sure that's wrong but not sure why I feel so sure. I have an inkling that the multiplication of probabilities here presumes independence, that $\Pr(P_1-P_2\le\epsilon)$ is independent of $\Pr(P_2-P_3\le\epsilon)$. But $\Pr(P_1-P_3\le\epsilon)$ is categorically not independent of the first two.
And yet there's something in it. There is some way of describing the true draw probability, namely that all performances (1, 2 and 3 in this example, but an arbitrary number in practice) are within $\epsilon$ of all other performances.
But suddenly I'm at the end of my savvy and facing some pretty deep reading with not even any guarantee that there is a solution to this, given TrueSkill's authors parked it in the too-hard basket. Or that if there is one, that it's not easy to express or calculate (but a as the number of tied performances grows blows out in some nasty way).
But the optimist in me hopes there is a relatively lucid formulation. A gut feel tells me that as long as a new random variable (like $P_\delta$ above) can be defined as a linear combination of the tied performances that can (perforce would) use $\epsilon$ can be tested in relation to some condition (being negative, positive, less than or greater than some constant) that there is hope. I'm lacking inspiration in that space today mind you ;-).
 A: The condition that all 3 performances are within $\epsilon$ can be written as follows:
$$
P_2 \in [P_1-\epsilon, P_1+\epsilon] \\
P_3 \in [P_1-\epsilon, P_1+\epsilon] \\
P_3 \in [P_2-\epsilon, P_2+\epsilon]
$$
The last two conditions can be merged into the one condition:
$$
P_3 \in [\max(P_1,P_2)-\epsilon, \min(P_1,P_2)+\epsilon]
$$
Thus the probability that 3 teams are tied is:
$$
Pr(P_1+\epsilon \ge P_2 \ge P_1-\epsilon, \min(P_1,P_2)+\epsilon \ge P_3 \ge \max(P_1,P_2)-\epsilon)
$$
This formula can be simplified through a series of steps.
The first step is to eliminate the min/max by summing over the two cases $P_1 \ge P_2$ and $P_1 < P_2$:
$$
Pr(P_1 \ge P_2 \ge P_1-\epsilon, P_2+\epsilon \ge P_3 \ge P_1-\epsilon)
+ \\ Pr(P_1+\epsilon \ge P_2 > P_1, P_1+\epsilon \ge P_3 \ge P_2-\epsilon)
$$
The first term can be simplified as follows:
$$
Pr(P_1 \ge P_2 \ge P_1-\epsilon, P_2+\epsilon \ge P_3 \ge P_1-\epsilon)
= \\
Pr(P_1 \ge P_2 \ge P_1-\epsilon, P_2+\epsilon \ge P_3)
- Pr(P_1 \ge P_2 \ge P_1-\epsilon, P_3 < P_1-\epsilon)
$$
The first term in this subtraction can be further simplified:
$$
Pr(P_1 \ge P_2 \ge P_1-\epsilon, P_2+\epsilon \ge P_3)
= \\
Pr(P_1 \ge P_2, P_2+\epsilon \ge P_3)
- Pr(P_2 < P_1-\epsilon, P_2+\epsilon \ge P_3)
= \\
Pr(P_1 \ge P_2 \ge P_3-\epsilon)
- Pr(P_1-\epsilon > P_2 \ge P_3-\epsilon)
$$
The term $Pr(P_1 \ge P_2 \ge P_3-\epsilon)$ can be evaluated using the bivariate normal cdf, because it is equivalent to $Pr(x \le 0, y \le 0)$ where $x= P_2-P_1, y=P_3-\epsilon-P_2$ are correlated Gaussian random variables.  The other term $Pr(P_1-\epsilon > P_2 \ge P_3-\epsilon)$ can also be evaluated using the bivariate normal cdf.  Similarly, we can simplify the other term above:
$$
Pr(P_1 \ge P_2 \ge P_1-\epsilon, P_3 < P_1-\epsilon)
= \\
Pr(P_1 \ge P_2, P_3 < P_1-\epsilon)
- Pr(P_2 < P_1-\epsilon, P_3 < P_1-\epsilon)
$$
These two terms can be evaluated using the bivariate normal cdf.
In this way, all of the probabilities above can be simplified down to sums and differences of bivariate normal cdf terms, with 8 terms in all.
