# Matching, Optimal Transport and Statistical Tests

To explain the “optimal transport” problem, we usually start with Gaspard Monge’s “Mémoire sur la théorie des déblais et des remblais“, where the the problem of transporting a given distribution of matter (a pile of sand for instance) into another (an excavation for instance). This problem is usually formulated using distributions, and we seek the “optimal” transport from one distribution to the other one. The formulation, in the context of distributions has been formulated in the 40’s by Leonid Kantorovich, e.g. from the distribution on the left to the distribution on the right. Consider now the context of finite sets of points. We want to transport mass from points $\{A_1,\cdots,A_4\}$ to points $\{B_1,\cdots,B_4\}$. It is a complicated combinatorial problem. For 4 points, there are only 24 possible transfer to consider, but it exceeds 20 billions with 15 points (on each side). For instance, the following one is usually seen as inefficient while the following is usually seen as much better Of course, it depends on the cost of the transport, which depends on the distance between the origin and the destination. That cost is usually either linear or quadratic.

There are many application of optimal transport in economics, see eg Alfred’s book Optimal Transport Methods in Economics. And there are also applications in statistics, that what I’ve seen while I was discussing with Pierre while I was in Boston, in June. For instance if we want to test whether some sample were drawn from the same distribution,

set.seed(13) npoints <- 25 mu1 <- c(1,1) mu2 <- c(0,2) Sigma1 <- diag(1, 2, 2) Sigma2 <- diag(1, 2, 2) Sigma2[2,1] <- Sigma2[1,2] <- -0.5 Sigma1 <- 0.4 * Sigma1 Sigma2 <- 0.4 *Sigma2 library(mnormt) X1 <- rmnorm(npoints, mean = mu1, Sigma1) X2 <- rmnorm(npoints, mean = mu2, Sigma2) plot(X1[,1], X1[,2], ,col="blue") points(X2[,1], X2[,2], col = "red") Here we use a parametric model to generate our sample (as always), and we might think of a parametric test (testing whether mean and variance parameters of the two distributions are equal). or we might prefer a nonparametric test. The idea Pierre mentioned was based on optimal transport. Consider some quadratic loss

ground_p <- 2 p <- 1 w1 <- rep(1/npoints, npoints) w2 <- rep(1/npoints, npoints) C <- cost_matrix_Lp(t(X1), t(X2), ground_p) library(transport) library(winference) a <- transport(w1, w2, costm = C^p, method = "shortsimplex")

then it is possible to match points in the two samples

nonzero <- which(a$mass != 0) from_indices <- a$from[nonzero] to_indices <- a$to[nonzero] for (i in from_indices){ segments(X1[from_indices[i],1], X1[from_indices[i],2], X2[to_indices[i], 1], X2[to_indices[i],2]) } Here we can observe two things. The total cost can be seen as rather large > cost=function(a,X1,X2){ nonzero <- which(a$mass != 0) naa=a[nonzero,] d=function(i) (X1[naa$from[i],1]-X2[naa$to[i],1])^2+(X1[naa$from[i],2]-X2[naa$to[i],2])^2 sum(Vectorize(d)(1:npoints)) } > cost(a,X1,X2)  9.372472

and the angle of the transport direction is alway in the same direction (more or less)

> angle=function(a,X1,X2){ nonzero <- which(a$mass != 0) naa=a[nonzero,] d=function(i) (X1[naa$from[i],2]-X2[naa$to[i],2])/(X1[naa$from[i],1]-X2[naa$to[i],1]) atan(Vectorize(d)(1:npoints)) } > mean(angle(a,X1,X2))  -0.3266797 > library(plotrix) > ag=(angle(a,X1,X2)/pi)*180 > ag[ag<0]=ag[ag<0]+360 > dag=hist(ag,breaks=seq(0,361,by=1)-.5) > polar.plot(dag$counts,seq(0,360,by=1),main=”Test Polar Plot”,lwd=3,line.col=4)

(actually, the following plot has been obtain by generating a thousand of sample of size 25) In order to have a decent test, we need to see what happens under the null assumption (when drawing samples from the same distribution), see Here is the optimal matching Here is the distribution of the total cost, when drawing a thousand samples,

VC=rep(NA,1000) VA=rep(NA,1000*npoints) for(s in 1:1000){ X1a <- rmnorm(npoints, mean = mu1, Sigma1) X1b <- rmnorm(npoints, mean = mu1, Sigma2) ground_p <- 2 p <- 1 w1 <- rep(1/npoints, npoints) w2 <- rep(1/npoints, npoints) C <- cost_matrix_Lp(t(X1a), t(X1b), ground_p) ab <- transport(w1, w2, costm = C^p, method = "shortsimplex") VC[s]=cout(ab,X1a,X1b) VA[s*npoints-(0:(npoints-1))]=angle(ab,X1a,X1b) } plot(density(VC) So our cost of 9 obtained initially was not that high. Observe that when drawing from the same distribution, there is now no pattern in the optimal transport

ag=(VA/pi)*180 ag[ag<0]=ag[ag<0]+360 dag=hist(ag,breaks=seq(0,361,by=1)-.5) polar.plot(dag\$counts,seq(0,360,by=1),main="Test Polar Plot",lwd=3,line.col=4) Nice isn’t it? I guess I will spend some time next year working on those transport algorithm, since we have great R packages, and hundreds of applications in economics…

# Jeu de mains, jeu de malins ?

Ce week-end – merci à Serge qui m’a raconté cette histoire – beaucoup de mes illusions de jeunesse se sont envolées…. Mais gardons un peu de suspens pour que mes lecteurs lisent le billet en entier…. • Pierre, ciseaux ou feuille

Le jeu “pierre, ciseaux, feuille” semble être né en Chine (soit sous le joli nom 五雜俎, alias wuzazu, datant de la fin de période Ming, soit 手勢令 i.e.shoushiling) ou au Japon (soit sous la forme 一二三 , i.e. Hi-fu-mu, c’est à dire un, deux, trois,soit じゃんけんぽん ou plus simplement じゃんけん, alias janken  , lui même basé sur deux anciens jeux, 数拳 – su ken – et 三すくみ拳 -saan sukumi ken). Pour rappel, la pierre (ぐう) bat les ciseaux (ちょき) qui bat à la feuille (ぱあ) qui à son tour bat la pierre. Bref, ce jeux vient d’orient !
Pour preuve, en avril 2005, Takashi Hashiyama, alors président de Maspro Denkoh Corporation, une entreprise d’électronique située à Nagoya, au Japon, n’arrivait pas à savoir chez qui (entre Christie’s et Sotheby’s) il allait vendre sa collection d’art (estimée à plus de 20 millions de dollars). N’arrivant pas à choisir, le choix s’est fait par “pierre, ciseaux, feuille” (ici).

• Le théorème min-max

Pour comprendre ce jeu, il faut revenir un peu sur des résultats classiques de théorie des jeux. La matrice de gain pour le joueur X est ici

 GAINS X joue « pierre » X joue « feuille » X joue « ciseaux » Y joue « pierre » 0 1 -1 Y joue « feuille » -1 0 1 Y joue « ciseaux » 1 -1 0

et de manière symétrique, on peut définir la matrice de gain pour le joueur Y. Notons que l’on suppose ici un “jeu à somme nulle” (i.e. ce que perd X,  Y le gagne, et réciproquement). Cette seconde matrice de gains est alors l’opposé de celle-ci.
La stratégie de base du jeu a été énoncée par Italo Calvino sous la forme suivante, “tu sais que ce que tu peux espérer de mieux est d’éviter le pire”.
Intuitivement, supposons que X souhaite jouer “pierre” on peut penser que Y reconstitue les anticipations de X, et donc joue “feuille”, mais à son tour X doit modifier ses anticipations et donc jouer “ciseaux”, plutôt que “feuille”, etc. Bref, rien ne se stabilise et aucune stratégie pure ne s’impose.
Mais s’il n’y a pas de stratégie pure, tentons des stratégies mixtes, i.e. choisissons intelligemment une mesure de probabilités sur l’ensemble des stratégies pures. En pratique, on jouera alors aléatoirement, en pondérant chaque tirage suivant la mesure de probas choisie. Une stratégie mixte (pour X) est alors caractérisée par un vecteur du simplexe, La  stratégie mixte consiste à décider au  hasard quel coup on jouera avec équiprobabilité.
Formellement, notons que le tirage équiprobable – par les deux joueurs – est un équilibre de Nash. En effet, si A désigne la matrice de gain pour X (on est sur des jeux à somme nulle), John von Neumann a montré ce que l’on appelle aujourd’hui le théorème du minimax (en fait, Emile Borel avait également obtenu ce résultat) en 1926, sous la forme suivante: il y a un équilibre de Nash associé à A si et seulement si L’intuition est assez évidente (mais la preuve de ce théorème du point fixe est un peu technique): est un équilibre de Nash si aucun des jours n’a intérêt de dévier dans sa stratégie, i.e. et On recherche alors un point selle dans la matrice des gains.
Dans le cas de notre jeu de “pierre, ciseaux, feuille“, notons que si alors , et on note que ainsi que . Aussi, pour toute  stratégie mixte, et de même Bref, cette stratégie est une équilibre de Nash.

• Pierre, ciseaux, feuille, ou comment trouver une stratégie

Bref, comme on vient de le noter, la meilleure stratégie à adopter quand on jour à “pierre, caillou, feuille” c’est de ne pas avoir de stratégie… ce qui revient à dire que l’issue du jeu est purement aléatoire. Pour Kenneth Binmore (ici) la méthode du minimax consiste à dire que face à un “bon” joueur, il n’y a rien de mieux à faire de que choisir la stratégie mixte du théorème de von Neumann.
Bon, si on réfléchit deux minutes, on sait qu’on ne joue pas une fois à ce jeux, mais que l’on répète le jeux quelques fois (avant de se lasser). C’est ce qu’a décrit Julia Robinson en 1950. L’idée était d’avoir une approche que l’on peut qualifier de bayésienne, consistant à utiliser pour la stratégie de l’adversaire les probabilités observées.

• Un champion du monde de “pierre, ciseaux, feuille” ? Maintenant que j’ai fini avec les rappels de théorie, revenons à la discussion de samedi soir…. J’ai appris qu’il existait un champion du monde de “pierre, ciseaux, feuille” ! Serge prétendait – samedi soir – qu’un joueur avait été vainqueur deux années de suite (ce qui m’a vraiment fait douté sur le côté aléatoire du jeu… malheureusement pour mon histoire, ce n’est qu’une légende, car les champions du monde entre 2002 et 2008 ont été Peter Lovering, Rob Krueger, Lee Rammage, Andrew Bergel, Bob Cooper, Andrea Farina et Monica Martinez).
Comment peut-il y avoir un champion du monde à un jeu qui repose exclusivement sur du hasard ? C’est comme si on inventait un champion du monde de lancer de pièces (ou de dés), qui serait récompensé pour sa régularité1 !
En creusant un peu, j’ai aussi découvert qu’il existait un livre qui apprend à gagner au jeu (je ne l’ai encore lu, mais promis, si la BU le commande, je le dévorerais et ferais un compte rendu détaillé ici même).