Reproduce Cross Correlation results in Python

Question

I'm trying and failing to reproduce the results found in this paper, where the cross correlation between two distributions is calculated:

p(XY) = 0.5 * P(X = Y) + 1 * P(X > Y)

Here' the python code (x and y are taken from the first table on page 5):

import numpy as np
x = [-0.01, -0.05, -0.13, -0.02, -0.17, -0.09, 0.0, -0.04, -0.03, 0.06, -1.37, -0.03, 0.01, -0.57, 0.04, -0.09, -0.04, -5.56, -0.02, 0.0, 0.0, 0.0, 0.0, -0.03, -0.55, -2.6, -0.42, -1.35, 0.0, 0.43, -0.74, -0.47, 0.0, -10.25, -11.18, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

y = [0.0, -0.16, -0.11, -0.04, -0.1, 0.02, -0.04, -0.13, -0.03, -0.04, 0.0, -0.01, 0.01, -0.38, 0.13, -0.26, -0.11, -0.09, 0.0, 0.0, 0.18, 0.0, -0.02, 0.03, -0.17, 0.21, 0.17, -0.42, 0.0, 0.87, 0.15, 0.0, 0.47, 3.73, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
bins = np.arange(min(x + y), max(x + y), 0.01)

# x histogram (normalized to sum 1.0)
hx, dx = np.histogram(x, bins)
hx = 1.0*hx/np.sum(hx)  # max(hx) = 0.3170731707317073

# y histogram (normalized to sum 1.0)
hy, dy = np.histogram(y, bins)
hy = 1.0*hy/np.sum(hy) # max(hy) = 0.34999999999999998

# Now comes the hard part, the result doesn't match with the paper (0.345)
# this gives 1.01219512
np.correlate(hx, hy)*0.5 + sum(np.correlate(hx, hy, 'same')) 

# this gives 1.06097561
np.correlate(hx, hy)*0.5 + sum(np.correlate(hx, hy, 'full'))

What i'm doing wrong? Maybe i'm summing up all the lags? (positive and negative) Thanks for you help!

Answer 1

Background

The referenced paper is about evaluating the strength of chess players. The data consist of computer evaluations of the strengths of individual moves (for one game in this example but, in practice, over multiple games), with $x_1, \ldots, x_n$ the values of the moves by the first player and $y_1, \ldots, y_m$ the values of the moves by the second player. These are being considered as if they were each random samples of move-strength distributions $X$ and $Y$ associated with the players. (This assumption is debatable but won't be questioned in the following analysis.)

It is argued that the expected outcome of a game between those two players (which is one point for a win by $X$, half a point for a draw, and no points for a loss by $X$) is the expectation of $g(X,Y)$ where

$$g(X,Y)=1\, (X\gt Y);\quad g(X,Y)=1/2\, (X=Y);\quad g(X,Y)=0\, (X\lt Y).$$

Although this is not the expected outcome--the outcome is determined by the sums of the move values, which have a different distribution than the move values themselves--it may be a decent approximation to the expected outcome.

The question concerns how to estimate $\mathbb{E}(g(X,Y))$ from the data.

Solution

To make this estimate, use the empirical distribution functions of $X$ and $Y$, assuming $X$ and $Y$ are independent. This means the estimate is

$$\frac{1}{nm} \sum_{i=1}^n\sum_{j=1}^m g(x_i,y_i).$$

This can be computed exactly as written (as a double loop) or as an outer product. The latter in pseudocode (which will actually execute in R, for instance) could be written

z <- outer(x, -y, '+')
(sum(z > 0) + sum(z==0)/2)/length(z)

Its value for the sample data is $0.3450327$. It agrees with the number $0.345$ given in the paper to its full precision, suggesting this calculation correctly reflects the one intended in the paper.

---

For rapid calculation the expectation can be expressed in terms of the ranks of the combined data, thereby requiring $O((n+m)\log(n+m))$ effort instead of the $O(nm)$ effort implied by the outer product calculation. For instance, when there are no ties between any of the $x_i$ and $y_j$, let $r(i)$ be the rank of $x_i$ when the combined dataset is sorted in increasing order. Then the number of the $y_j$ less than $x_i$ is the number of combined data less than $x_i$ (equal to $r(i)-1$), minus the number of $x_{i^\prime}$ less than $x_i$ (equal to $i-1$), which is just $r(i)-i$, whence

$$\eqalign{ \sum_{i=1}^n\sum_{j=1}^m g(x_i,y_i) &= \sum_{x_i \gt y_j}1 \\ &= \sum_{i=1}^{n} \left(r(x_i) - i\right) \\ &= \sum_{i=1}^n r(x_i) - \frac{1}{2}n(n+1). }$$

When there are ties, a more complicated version of this formula will adjust for them. It would be equivalent to assigning the midrank to all groups of ties, as with this R code, whose value is the same as before:

(sum(rank(c(x,y), ties.method="average")[1:n]) - n*(n+1)/2) / (m*n)

Answer 2

Just as a note, here is the Python solution based on @whuber code:

import numpy as np
x = [-0.01, -0.05, -0.13, -0.02, -0.17, -0.09, 0.0, -0.04, -0.03, 0.06, -1.37, -0.03, 0.01, -0.57, 0.04, -0.09, -0.04, -5.56, -0.02, 0.0, 0.0, 0.0, 0.0, -0.03, -0.55, -2.6, -0.42, -1.35, 0.0, 0.43, -0.74, -0.47, 0.0, -10.25, -11.18, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
y = [0.0, -0.16, -0.11, -0.04, -0.1, 0.02, -0.04, -0.13, -0.03, -0.04, 0.0, -0.01, 0.01, -0.38, 0.13, -0.26, -0.11, -0.09, 0.0, 0.0, 0.18, 0.0, -0.02, 0.03, -0.17, 0.21, 0.17, -0.42, 0.0, 0.87, 0.15, 0.0, 0.47, 3.73, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

cr = np.subtract.outer(x, y)
p = (np.sum(cr > 0) + np.sum(cr == 0)*0.5)/cr.size
print p # 0.34503271861986912

Or without numpy

cr = [[xi - yi for yi in y] for xi in x]
cr = [item for sublist in cr for item in sublist]
p = (sum([g > 0 for g in cr]) + sum([g == 0 for g in cr])*0.5)/len(cr)
print p # 0.3450327186198691