<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/couplings-of-normal-variables/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://statisfaction.wordpress.com/2018/08/24/couplings-of-normal-variables/"> R â€“ Statisfaction</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><img data-orig-file="https://statisfaction.files.wordpress.com/2018/08/2018-gganimate-coupling.gif?w=456" data-orig-size="480,480" data-medium-file="https://statisfaction.files.wordpress.com/2018/08/2018-gganimate-coupling.gif?w=456?w=400" data-attachment-id="7021" data-image-meta="{" data-image-description="" src="https://statisfaction.files.wordpress.com/2018/08/2018-gganimate-coupling.gif?w=456" data-permalink="https://statisfaction.wordpress.com/2018/08/24/couplings-of-normal-variables/2018-gganimate-coupling/" data-large-file="https://statisfaction.files.wordpress.com/2018/08/2018-gganimate-coupling.gif?w=456?w=480" alt="2018-gganimate-coupling" data-image-title="2018-gganimate-coupling" data-comments-opened="1" class="wp-image-7021" data-recalc-dims="1"></p>
<p>Hi,</p>
<p>Just to play a bit with the <a rel="nofollow" href="https://github.com/thomasp85/gganimate" target="_blank">gganimate</a> package, and to celebrate National Coupling Day, the above plot shows different couplings of two univariate Normal distributions, Normal(0,1) and Normal(2,1). That is, each point is a pair (x,y) where x follows a Normal(0,1) and y follows a Normal(2,1). Below I’ll recall briefly how each coupling operates, in the Normal case. The code is available at the end of the post.</p>
<p><span></span></p>
<p>To create the plots, x-values were drawn from a Normal(0,1) distribution and kept fixed. Only the y-values change from a coupling to the next.</p>
<p>The optimal transport coupling is the simplest here. The y-values are simply defined by y = x + 2, i.e. they are translated by 2, the difference between the means of Normal(0,1) and Normal(2,1). It leads to the smallest expected value for the distance between y and x, over all couplings of the given Normals. The multivariate Normal case is well explained in <a rel="nofollow" href="http://djalil.chafai.net/blog/2010/04/30/wasserstein-distance-between-two-gaussians/" target="_blank">this blog post</a>Ã‚Â byÃ‚Â Djalil ChafaÃƒÂ¯.</p>
<p>The reflection coupling is also simple. Each point is reflected with respect to the midpoint between the means. Here it just amounts to defining y = 2 – x. Reflection couplings can be used to study stochastic processses, for instance multivariate diffusionsÃ‚Â as in <a rel="nofollow" href="https://projecteuclid.org/euclid.aop/1176992442" target="_blank">this important 1986 paper by Lindvall and Rogers</a>, and <a rel="nofollow" href="https://link.springer.com/article/10.1007/s00440-015-0673-1" target="_blank">this more recent one by Eberle in 2016</a>.</p>
<p>A variant of the reflection coupling, termed “reflection-maximal” coupling in the above plot, appears in this <a rel="nofollow" href="https://arxiv.org/abs/1805.00452" target="_blank">recent study of Hamiltonian Monte Carlo</a>Ã‚Â by Bou-Rabee, Eberle and Zimmer. It’s a maximal coupling with a reflection element to it (see Eq. (15) in that article for a precise definition). Essentially, there is an attempt at accepting x as a draw from a Normal(2,1). If this does not succeed, then y is defined as 2 – x (a reflection). The procedure is a maximal coupling because the probability of the event {x = y} is maximal over all possible couplings.</p>
<p>Finally the coupling referred to as “maximal coupling” in the plot is obtained with the procedure described in <a rel="nofollow" href="https://statisfaction.wordpress.com/2017/09/06/sampling-from-a-maximal-coupling/" target="_blank">this earlier post</a>. It is also a maximal coupling (these are not unique), but the distribution of y given that “y is not x” is independent of x.</p>
<p>That’s it! The code is below.</p>
<pre>

library(ggplot2)
library(gganimate)
set.seed(1)

# normal means
mu1 &lt;- 0
mu2 &lt;- 2
# std deviation
sigma &lt;- 1
# number of samples
nsamples &lt;- 1000
# reflection-maximal coupling
reflmax_samples &lt;- matrix(nrow = nsamples, ncol = 2)
# draw x components first
reflmax_samples[,1] &lt;- rnorm(nsamples)
# this follows the notation of Bou-Rabee et al, 2018, roughly
z &lt;- mu1 - mu2
normz &lt;- sqrt(sum(z^2))
e &lt;- z / normz
utilde &lt;- runif(nsamples, 0, 1)
accepts &lt;- (log(utilde) &lt; (dnorm(e * reflmax_samples[,1] + normz, 0, 1, log = TRUE) - dnorm(e*reflmax_samples[,1], log = TRUE)))
eta &lt;- reflmax_samples[,1]
eta[accepts] &lt;- reflmax_samples[accepts,1] + z
eta[!accepts]&lt;- reflmax_samples[!accepts,1] - 2 * (e * reflmax_samples[!accepts,1]) * e
reflmax_samples[,2] &lt;- mu2 + eta

df &lt;- data.frame(coupling = rep(&quot;reflection-maximal&quot;, nsamples), x = reflmax_samples[,1], y = reflmax_samples[,2])

# reflection coupling
refl_samples &lt;- matrix(0, nrow = nsamples, ncol = 2)
refl_samples[,1] &lt;- reflmax_samples[,1]
refl_samples[,2] &lt;- (mu2-mu1)-refl_samples[,1]

df &lt;- rbind(df, data.frame(coupling = rep(&quot;reflection&quot;, nsamples), x = refl_samples[,1], y = refl_samples[,2]))

# optimal transport coupling
transport_samples &lt;- matrix(0, nrow = nsamples, ncol = 2)
transport_samples[,1] &lt;- reflmax_samples[,1]
transport_samples[,2] &lt;- mu2 - mu1 + transport_samples[,1]

df &lt;- rbind(df, data.frame(coupling = rep(&quot;optimal transport&quot;, nsamples), x = transport_samples[,1], y = transport_samples[,2]))

# max coupling
max_samples &lt;- matrix(0, nrow = nsamples, ncol = 2)
max_samples[,1] &lt;- reflmax_samples[,1]
dp &lt;- function(x) dnorm(x, mean = mu1, sd = 1, log = TRUE)
dq &lt;- function(x) dnorm(x, mean = mu2, sd = 1, log = TRUE)
rq &lt;- function(n) rnorm(n, mean = mu2, sd = 1)
for (isample in 1:nsamples){
x &lt;- max_samples[isample,1]
if (dp(x) + log(runif(1)) &lt; dq(x)){
max_samples[isample,2] &lt;- x
} else {
reject &lt;- TRUE
y &lt;- NA
while (reject){
y &lt;- rq(1)
reject &lt;- (dq(y) + log(runif(1)) &lt; dp(y))
}
max_samples[isample,2] &lt;- y
}
}
df &lt;- rbind(df, data.frame(coupling = rep(&quot;maximal&quot;, nsamples), x = max_samples[,1], y = max_samples[,2]))

## Scatter plots and marginals
# ggplot(df, aes(x=x, y=y, group = coupling, colour = factor(coupling))) +
# geom_point()+
# theme_minimal() + viridis::scale_color_viridis(discrete=T)
#
# ggplot(df, aes(x=x, group = coupling, fill = factor(coupling))) + geom_histogram(aes(y = ..density..), position = position_dodge()) +
# theme_minimal() + viridis::scale_fill_viridis(discrete=T)
#
# ggplot(df, aes(x=y, group = coupling, fill = factor(coupling))) + geom_histogram(aes(y = ..density..), position = position_dodge()) +
# theme_minimal() + viridis::scale_fill_viridis(discrete=T)

# gganimate

ggplot(df, aes(x = x, y = y)) + xlim(-4, 4) + ylim(-5, 6) +
geom_point()+ geom_text(data = data.frame(coupling = unique(df$coupling)), aes(label = coupling, x = -1, y = -4), size = 10) +
theme_minimal() +
transition_states(coupling, 3, 1)

</pre>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://statisfaction.wordpress.com/2018/08/24/couplings-of-normal-variables/"> R â€“ Statisfaction</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 alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/JfQ-C6-uLTo" height="1"></a></div>