<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="http://www.r-bloggers.com/a-simpson-paradox-of-sorts/"></div></div><div><a data-url="http://www.r-bloggers.com/a-simpson-paradox-of-sorts/" data-counturl="http://www.r-bloggers.com/a-simpson-paradox-of-sorts/" data-count="horizontal" data-text="a Simpson paradox of sorts" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://xianblog.wordpress.com/2016/05/06/a-simpson-paradox-of-sorts/"> R – Xi'an's Og</a></strong>, and kindly contributed to <a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><strong>T</strong>he riddle from <a target="_blank" href="http://fivethirtyeight.com/features/the-perplexing-puzzle-of-the-proud-partygoers/">The Riddler</a> this week is about finding an undirected graph with N nodes and no isolated node such that the number of nodes with more connections than the average of their neighbours is maximal. A representation of a connected graph is through a matrix X of zeros and ones, on which one can spot the nodes satisfying the above condition as the positive entries of the vector (X<strong>1</strong>)^2-(X^2<strong>1</strong>), if <strong>1</strong> denotes the vector of ones. I thus wrote an R code aiming at optimising this target</p>
<pre>
targe &lt;- function(F){
  sum(F%*%F%*%rep(1,N)/(F%*%rep(1,N))^2&lt;1)}
</pre>
<p>by mere simulated annealing:</p>
<pre>
rate &lt;- function(N){ 
# generate matrix F
# 1. no single 
F=matrix(0,N,N) 
F[sample(2:N,1),1]=1 
F[1,]=F[,1] 
for (i in 2:(N-1)){ 
if (sum(F[,i])==0) 
F[sample((i+1):N,1),i]=1 
F[i,]=F[,i]} 
if (sum(F[,N])==0) 
F[sample(1:(N-1),1),N]=1 
F[N,]=F[,N] 
# 2. more connections 
F[lower.tri(F)]=F[lower.tri(F)]+
  sample(0:1,N*(N-1)/2,rep=TRUE,prob=c(N,1)) 
F[F>1]=1
F[upper.tri(F)]=t(F)[upper.tri(t(F))]
#simulated annealing
T=1e4
temp=N
targo=targe(F)
for (t in 1:T){
  #1. local proposal
  nod=sample(1:N,2)
  prop=F
  prop[nod[1],nod[2]]=prop[nod[2],nod[1]]=
     1-prop[nod[1],nod[2]]
  while (min(prop%*%rep(1,N))==0){
    nod=sample(1:N,2)
    prop=F
    prop[nod[1],nod[2]]=prop[nod[2],nod[1]]=
     1-prop[nod[1],nod[2]]}
  target=targe(prop)
  if (log(runif(1))*temp&lt;target-targo){ 
    F=prop;targo=target} 
#2. global proposal 
  prop=F prop[lower.tri(prop)]=F[lower.tri(prop)]+
   sample(c(0,1),N*(N-1)/2,rep=TRUE,prob=c(N,1)) 
prop[prop>1]=1
  prop[upper.tri(prop)]=t(prop)[upper.tri(t(prop))]
  target=targe(prop)
  if (log(runif(1))*temp&lt;target-targo){
      F=prop;targo=target}
   temp=temp*.999
   }
return(F)}
</pre>
<p><a target="_blank" href="https://xianblog.files.wordpress.com/2016/05/simpson.png"><img height="264" alt="Eward Simpson" width="222" class="wp-image-30561" src="https://xianblog.files.wordpress.com/2016/05/simpson.png?w=222&h=264"></a>This code returns quite consistently (modulo the simulated annealing uncertainty, which grows with N) the answer N-2 as the number of entries above average! Which is rather surprising in a Simpson-like manner since all entries but two are above average. (Incidentally, I found out that Edward Simpson recently wrote a paper in <a target="_blank" href="http://onlinelibrary.wiley.com/doi/10.1111/j.1740-9713.2010.00424.x/full">Significance</a> about the Simpson-Yule paradox and him being a member of the <a target="_blank" href="http://onlinelibrary.wiley.com/doi/10.1111/j.1740-9713.2010.00424.x/full">Bletchley Park Enigma</a> team. I must have missed out the connection with the Simpson paradox when reading the paper in the first place…)</p>
<p>Filed under: <a target="_blank" href="https://xianblog.wordpress.com/category/books/">Books</a>, <a target="_blank" href="https://xianblog.wordpress.com/category/kids/">Kids</a>, <a target="_blank" href="https://xianblog.wordpress.com/category/pictures/">pictures</a>, <a target="_blank" href="https://xianblog.wordpress.com/category/statistics/r-statistics/">R</a> Tagged: <a target="_blank" href="https://xianblog.wordpress.com/tag/bletchley-park/">Bletchley Park</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/edward-simpson/">Edward Simpson</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/enigma-code-machine/">Enigma code machine</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/graph/">graph</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/mathematical-puzzle/">mathematical puzzle</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/significance/">Significance</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/simpsons-paradox/">Simpson’s paradox</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/simulated-annealing/">simulated annealing</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/the-riddler/">The Riddler</a>, <a target="_blank" href="https://xianblog.wordpress.com/tag/yule/">Yule</a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/xianblog.wordpress.com/30557/"><img alt="" src="http://feeds.wordpress.com/1.0/comments/xianblog.wordpress.com/30557/" border="0"></a> <img height="1" alt="" width="1" src="https://pixel.wp.com/b.gif?host=xianblog.wordpress.com&blog=5051449&%23038;post=30557&%23038;subd=xianblog&%23038;ref=&%23038;feed=1" border="0"></p>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="http://www.r-bloggers.com/a-simpson-paradox-of-sorts/"></div></div><div><a data-url="http://www.r-bloggers.com/a-simpson-paradox-of-sorts/" data-counturl="http://www.r-bloggers.com/a-simpson-paradox-of-sorts/" data-count="vertical" data-text="a Simpson paradox of sorts" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://xianblog.wordpress.com/2016/05/06/a-simpson-paradox-of-sorts/"> R – Xi'an's Og</a></strong>.</div>
<hr>
<a rel="nofollow" href="http://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="http://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="http://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="http://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="http://www.r-bloggers.com/?s=Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="http://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="http://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="http://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="http://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="http://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="http://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="http://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="http://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="http://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="http://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="http://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="http://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="http://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="http://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/RduWu3OOyS8"></a></div>