<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/matching-optimal-transport-and-statistical-tests/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://freakonometrics.hypotheses.org/50958"> R-english – Freakonometrics</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>To explain the “optimal transport” problem, we usually start with Gaspard Monge’s “<i>Mémoire sur la théorie des déblais et des remblais</i>“, where the the problem of transporting a given distribution of matter (a pile of <span>sand</span> 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.</p>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_131.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Consider now the context of finite sets of points. We want to transport mass from points <img alt="https://latex.codecogs.com/gif.latex?%5C%7BA_1%2C%5Ccdots%2CA_4%5C%7D" src="https://latex.codecogs.com/gif.latex?%5C%7BA_1%2C%5Ccdots%2CA_4%5C%7D"> to points <img alt="https://latex.codecogs.com/gif.latex?{B_1,cdots,B_4}" src="https://latex.codecogs.com/gif.latex?%5C%7BB_1,%5Ccdots,B_4%5C%7D">. 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</p>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_130.png?resize=322%2C255&ssl=1" data-recalc-dims="1"></p>
<p>while the following is usually seen as much better</p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_129.png?resize=323%2C259&ssl=1" data-recalc-dims="1"></p>
<p>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.</p>
<p>There are many application of optimal transport in economics, see eg Alfred’s book <a rel="nofollow" href="http://press.princeton.edu/titles/10870.html" target="_blank">Optimal Transport Methods in Economics.</a> And there are also applications in statistics, that what I’ve seen while I was discussing with <a rel="nofollow" href="https://sites.google.com/site/pierrejacob/" target="_blank">Pierre</a> while I was in Boston, in June. For instance if we want to test whether some sample were drawn from the same distribution,</p>
<p><code>set.seed(13)<br>
npoints &lt;- 25<br>
mu1 &lt;- c(1,1)<br>
mu2 &lt;- c(0,2)<br>
Sigma1 &lt;- diag(1, 2, 2)<br>
Sigma2 &lt;- diag(1, 2, 2)<br>
Sigma2[2,1] &lt;- Sigma2[1,2] &lt;- -0.5<br>
Sigma1 &lt;- 0.4 * Sigma1<br>
Sigma2 &lt;- 0.4 *Sigma2<br>
library(mnormt)<br>
X1 &lt;- rmnorm(npoints, mean = mu1, Sigma1)<br>
X2 &lt;- rmnorm(npoints, mean = mu2, Sigma2)<br>
plot(X1[,1], X1[,2], ,col=&quot;blue&quot;)<br>
points(X2[,1], X2[,2], col = &quot;red&quot;)</code></p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_112.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>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).</p>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_113.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>or we might prefer a nonparametric test. The idea Pierre mentioned was based on optimal transport. Consider some quadratic loss</p>
<p><code>ground_p &lt;- 2<br>
p &lt;- 1<br>
w1 &lt;- rep(1/npoints, npoints)<br>
w2 &lt;- rep(1/npoints, npoints)<br>
C &lt;- cost_matrix_Lp(t(X1), t(X2), ground_p)<br>
library(transport)<br>
a &lt;- transport(w1, w2, costm = C^p, method = &quot;shortsimplex&quot;)</code></p>
<p>then it is possible to match points in the two samples</p>
<p><code>nonzero &lt;- which(a$mass != 0)<br>
from_indices &lt;- a$from[nonzero]<br>
to_indices &lt;- a$to[nonzero]<br>
for (i in from_indices){<br>
segments(X1[from_indices[i],1], X1[from_indices[i],2], X2[to_indices[i], 1], X2[to_indices[i],2])<br>
}</code></p>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_114.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Here we can observe two things. The total cost can be seen as rather large</p>
<p><code>> cost=function(a,X1,X2){<br>
nonzero &lt;- which(a$mass != 0)<br>
naa=a[nonzero,]<br>
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<br>
sum(Vectorize(d)(1:npoints))<br>
}<br>
> cost(a,X1,X2)<br>
[1] 9.372472</code></p>
<p>and the angle of the transport direction is alway in the same direction (more or less)</p>
<p><code>> angle=function(a,X1,X2){<br>
nonzero &lt;- which(a$mass != 0)<br>
naa=a[nonzero,]<br>
d=function(i) (X1[naa$from[i],2]-X2[naa$to[i],2])/(X1[naa$from[i],1]-X2[naa$to[i],1])<br>
atan(Vectorize(d)(1:npoints))<br>
}<br>
> mean(angle(a,X1,X2))<br>
[1] -0.3266797</code></p>
<p>> library(plotrix)<br>
> ag=(angle(a,X1,X2)/pi)*180<br>
> ag[ag&lt;0]=ag[ag&lt;0]+360<br>
> dag=hist(ag,breaks=seq(0,361,by=1)-.5)<br>
> polar.plot(dag$counts,seq(0,360,by=1),main=”Test Polar Plot”,lwd=3,line.col=4)</p>
<p>(actually, the following plot has been obtain by generating a thousand of sample of size 25)</p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_120.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>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</p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_116.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Here is the optimal matching</p>
<p><img alt="" src="https://i1.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_115.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>Here is the distribution of the total cost, when drawing a thousand samples,</p>
<p><code>VC=rep(NA,1000)<br>
VA=rep(NA,1000*npoints)<br>
for(s in 1:1000){<br>
X1a &lt;- rmnorm(npoints, mean = mu1, Sigma1)<br>
X1b &lt;- rmnorm(npoints, mean = mu1, Sigma2)<br>
ground_p &lt;- 2<br>
p &lt;- 1<br>
w1 &lt;- rep(1/npoints, npoints)<br>
w2 &lt;- rep(1/npoints, npoints)<br>
C &lt;- cost_matrix_Lp(t(X1a), t(X1b), ground_p)<br>
ab &lt;- transport(w1, w2, costm = C^p, method = &quot;shortsimplex&quot;)<br>
VC[s]=cout(ab,X1a,X1b)<br>
VA[s*npoints-(0:(npoints-1))]=angle(ab,X1a,X1b)<br>
}<br>
plot(density(VC)</code></p>
<p><img alt="" src="https://i0.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_118.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>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</p>
<p><code>ag=(VA/pi)*180<br>
ag[ag&lt;0]=ag[ag&lt;0]+360<br>
dag=hist(ag,breaks=seq(0,361,by=1)-.5)<br>
polar.plot(dag$counts,seq(0,360,by=1),main=&quot;Test Polar Plot&quot;,lwd=3,line.col=4)</code></p>
<p><img alt="" src="https://i2.wp.com/freakonometrics.hypotheses.org/files/2017/07/S%C3%A9lection_119.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p> </p>
<p>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…</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://freakonometrics.hypotheses.org/50958"> R-english – Freakonometrics</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/ENUXeueyWho" height="1"></a></div>