Tag Archives: Rstats

Visualizing the Gender of US Senators With R and Highmaps

Collaborations, Maps, The World We Live In, , , ,

I wake up every morning in a house that was built by slaves (Michelle Obama)

Some days ago I was invited by the people of Highcharts to write a post in their blog. What I have done is a simple but revealing map of women senators of the United States of America. Briefly, this is what I’ve done to generate it:

  • read from the US senate website a XML file with senators info
  • clean and obtain gender of senators from their first names
  • summarize results by state
  • join data with a US geojson dataset to create the highmap

You can find details and R code here.

It is easy creating a highcharts using highcharter, an amazing library as genderizeR, the one I use to obtain gender names. I like them a lot.

Visualizing Stirling’s Approximation With Highcharts

Curiosities, Simulation, , ,

I said, “Wait a minute, Chester, you know I’m a peaceful man”, He said, “That’s okay, boy, won’t you feed him when you can” (The Weight, The Band)

It is quite easy to calculate the probability of obtaining the same number of heads and tails when tossing a coin N times, and N is even. There are 2^{N} possible outcomes and only C_{N/2}^{N} are favorable so the exact probability is the quotient of these numbers (# of favorable divided by # of possible).

There is another way to approximate this number incredibly well: to use the Stirling’s formula, which is 1/\sqrt{\pi\cdot N/2}

The next plot represents both calculations for N from 2 to 200. Although for small values of N, Stirling’s approximation tends to overestimate probability, you can see hoy is extremely precise as N becomes bigger:

James Stirling published this amazing formula in 1730. It simplifies the calculus to the extreme and also gives a quick way to obtain the answer to a very interesting question: how many tosses are needed to be sure that the probability of obtaining the same number of heads and tails is under any given threshold? Just solve the formula for N and you will obtain the answer. And, also, the formula is another example of the presence of pi in the most unexpected places, as happens here.

Just another thing: the more I use highcharter package the more I like it.

This is the code:

library(highcharter)
library(dplyr)
data.frame(N=seq(from=2, by=2, length.out = 100)) %>%
  mutate(Exact=choose(N,N/2)/2**N, Stirling=1/sqrt(pi*N/2))->data
hc <- highchart() %>% 
  hc_title(text = "Stirling's Approximation") %>% 
  hc_subtitle(text = "How likely is getting 50% heads and 50% tails tossing a coin N times?") %>% 
  hc_xAxis(title = list(text = "N: Number of tosses"), categories = data$N) %>% 
  hc_yAxis(title = list(text = "Probability"), labels = list(format = "{value}%", useHTML = TRUE)) %>% 
  hc_add_series(name = "Stirling", data = data$Stirling*100,  marker = list(enabled = FALSE), color="blue") %>% 
  hc_add_series(name = "Exact", data = data$Exact*100,  marker = list(enabled = FALSE), color="lightblue") %>% 
  hc_tooltip(formatter = JS("function(){return ('<b>Number of tosses: </b>'+this.x+'
<b>Probability: </b>'+Highcharts.numberFormat(this.y, 2)+'%')}")) %>%
  hc_exporting(enabled = TRUE) %>%
  hc_chart(zoomType = "xy")
hc

Women in Orchestras

The World We Live In, , , ,

I believe in the truth of fairy-tales more than I believe in the truth in the newspaper (Lotte Reiniger)

In my opinion, this graph is a visual demonstration that we live in a male chauvinist world.

Orchestras2

In this experiment I download the members of ten top orchestras of the world with the amazing rvest package. After cleaning texts, I obtain the gender of names with genderizeR package as I did here. Since I only take into account names genderized with high probability, these numbers cannot be exact. Apart of this, the plot speaks by itself.

setwd("YOUR WORKING DIRECTORY HERE")
library(rvest)
library(dplyr)
library(genderizeR)
read_html("http://www.berliner-philharmoniker.de/en/orchestra/") %&gt;%
html_nodes(".name") %&gt;%
html_text(trim=TRUE) %&gt;%
iconv("UTF-8") %&gt;%
gsub("[\r,\n]"," ", .) %&gt;%
gsub("\\s+", " ", .) %&gt;%
paste(collapse=" ") %&gt;%
findGivenNames() -&gt; berliner
saveRDS(berliner, file="berliner.RDS")
read_html("https://www.concertgebouworkest.nl/en/musicians") %&gt;%
html_nodes(".u-padding--b2") %&gt;%
html_text(trim=TRUE) %&gt;%
iconv("UTF-8") %&gt;%
gsub("\\s+", " ", .) %&gt;%
paste(collapse=" ") %&gt;%
findGivenNames() -&gt; rco
saveRDS(rco, file="rco.RDS")
read_html("http://www.philharmonia.spb.ru/en/about/orchestra/zkrasof/contents/") %&gt;%
html_nodes(".td") %&gt;%
html_text(trim=TRUE) %&gt;%
iconv("UTF-8") %&gt;%
gsub("[\r,\n]"," ", .) %&gt;%
gsub("\\s+", " ", .) %&gt;%
.[23] %&gt;%
findGivenNames() -&gt; spb
saveRDS(spb, file="spb.RDS")
read_html("http://ocne.mcu.es/conoce-a-la-ocne/orquesta-nacional-de-espana/componentes/") %&gt;%
html_nodes(".col-main") %&gt;%
html_text(trim=TRUE) %&gt;%
iconv("UTF-8") %&gt;%
gsub("[\r,\n]"," ", .) %&gt;%
gsub("\\s+", " ", .) %&gt;%
gsub("([[:lower:]])([[:upper:]][[:lower:]])", "\\1 \\2", .) %&gt;%
findGivenNames() -&gt; one
saveRDS(one, file="one.RDS")
read_html("http://www.gewandhausorchester.de/en/orchester/") %&gt;%
html_nodes("#content") %&gt;%
html_text(trim=TRUE) %&gt;%
iconv("UTF-8") %&gt;%
gsub("[\r,\n]"," ", .) %&gt;%
gsub("\\s+", " ", .) %&gt;%
findGivenNames() -&gt; leipzig
saveRDS(leipzig, file="leipzig.RDS")
read_html("http://www.wienerphilharmoniker.at/orchestra/members") %&gt;%
html_nodes(".ModSuiteMembersC") %&gt;%
html_text(trim=TRUE) %&gt;%
iconv("UTF-8") %&gt;%
gsub("[\r,\n,\t,*]"," ", .) %&gt;%
gsub("\\s+", " ", .) %&gt;%
gsub("([[:lower:]])([[:upper:]][[:lower:]])", "\\1 \\2", .) %&gt;%
paste(collapse=" ") %&gt;%
.[-18] %&gt;%
findGivenNames() -&gt; wiener
saveRDS(wiener, file="wiener.RDS")
read_html("http://www.laphil.com/philpedia/orchestra-roster") %&gt;%
html_nodes(".view-content") %&gt;%
html_text(trim=TRUE) %&gt;%
iconv("UTF-8") %&gt;%
gsub("\\s+", " ", .) %&gt;%
gsub("(?%
.[1] %&gt;%
findGivenNames() -&gt; laphil
saveRDS(laphil, file="laphil.RDS")
read_html("http://nyphil.org/about-us/meet/musicians-of-the-orchestra") %&gt;%
html_nodes(".resp-tab-content-active") %&gt;%
html_text(trim=TRUE) %&gt;%
iconv("UTF-8") %&gt;%
gsub("[\r,\n]"," ", .) %&gt;%
gsub("\\s+", " ", .) %&gt;%
gsub("(?%
findGivenNames() -&gt; nyphil
saveRDS(nyphil, file="nyphil.RDS")
urls=c("http://lso.co.uk/orchestra/players/strings.html",
"http://lso.co.uk/orchestra/players/woodwind.html",
"http://lso.co.uk/orchestra/players/brass.html",
"http://lso.co.uk/orchestra/players/percussion-harps-and-keyboards.html")
sapply(urls, function(x)
{
read_html(x) %&gt;%
html_nodes(".clearfix") %&gt;%
html_text(trim=TRUE) %&gt;%
iconv("UTF-8") %&gt;%
gsub("[\r,\n,\t,*]"," ", .) %&gt;%
gsub("\\s+", " ", .)
}) %&gt;% paste(., collapse=" ") %&gt;%
findGivenNames() -&gt; lso
saveRDS(lso, file="lso.RDS")
read_html("http://www.osm.ca/en/discover-osm/orchestra/musicians-osm") %&gt;%
html_nodes("#content-column") %&gt;%
html_text(trim=TRUE) %&gt;%
iconv("UTF-8") %&gt;%
gsub("[\r,\n]"," ", .) %&gt;%
gsub("\\s+", " ", .) %&gt;%
findGivenNames() -&gt; osm
saveRDS(osm, file="osm.RDS")
rbind(c("berliner", "Berliner Philharmoniker"),
c("rco", "Royal Concertgebouw Amsterdam"),
c("spb", "St. Petersburg Philharmonic Orchestra"),
c("one", "Orquesta Nacional de España"),
c("leipzig", "Gewandhaus Orchester Leipzig"),
c("wiener", "Wiener Philarmoniker"),
c("laphil", "The Los Angeles Philarmonic"),
c("nyphil", "New York Philarmonic"),
c("lso", "London Symphony Orchestra"),
c("osm", "Orchestre Symphonique de Montreal")) %&gt;% as.data.frame()-&gt; Orchestras
colnames(Orchestras)=c("Id", "Orchestra")
list.files(getwd(),pattern = ".RDS") %&gt;%
lapply(function(x)
readRDS(x) %&gt;% as.data.frame(stringsAsFactors = FALSE) %&gt;% cbind(Id=gsub(".RDS", "", x))
) %&gt;% do.call("rbind", .) -&gt; all
all %&gt;% mutate(probability=as.numeric(probability)) %&gt;%
filter(probability &gt; 0.9 &amp; count &gt; 15) %&gt;%
filter(!name %in% c("viola", "tuba", "harp")) %&gt;%
group_by(Id, gender) %&gt;%
summarize(Total=n())-&gt;all
all %&gt;% filter(gender=="female") %&gt;% mutate(females=Total) %&gt;% select(Id, females) -&gt; females
all %&gt;% group_by(Id) %&gt;% summarise(Total=sum(Total)) -&gt; total
inner_join(total, females, by = "Id") %&gt;% mutate(po_females=females/Total) %&gt;%
inner_join(Orchestras, by="Id")-&gt; df
library(ggplot2)
library(scales)
opts=theme(legend.position="none",
plot.background = element_rect(fill="gray85"),
panel.background = element_rect(fill="gray85"),
panel.grid.major.y=element_blank(),
panel.grid.major.x=element_line(colour="white", size=2),
panel.grid.minor=element_blank(),
axis.title = element_blank(),
axis.line.y = element_line(size = 2, color="black"),
axis.text = element_text(colour="black", size=18),
axis.ticks=element_blank(),
plot.title = element_text(size = 35, face="bold", margin=margin(10,0,10,0), hjust=0))
ggplot(df, aes(reorder(Orchestra, po_females), po_females)) +
geom_bar(stat="identity", fill="darkviolet", width=.5)+
scale_y_continuous(labels = percent, expand = c(0, 0), limits=c(0,.52))+
geom_text(aes(label=sprintf("%1.0f%%", 100*po_females)), hjust=-0.05, size=6)+
ggtitle(expression(atop(bold("Women in Orchestras"), atop("% of women among members", "")))) +
coord_flip()+opts

Playing With Julia (Set)

Drawings, Fractals, Image Processing, , , , ,

Viento, me pongo en movimiento y hago crecer las olas del mar que tienes dentro (Tercer Movimiento: Lo de Dentro, Extremoduro)

I really enjoy drawing complex numbers: it is a huge source of entertainment for me. In this experiment I play with the Julia Set, another beautiful fractal like this one. This is what I have done:

  • Choosing the function f(z)=exp(z3)-0.621
  • Generating a grid of complex numbers with both real and imaginary parts in [-2, 2]
  • Iterating f(z) over the grid a number of times so zn+1 = f(zn)
  • Drawing the resulting grid as I did here
  • Gathering all plots into a GIF with ImageMagick as I did in my previous post: each frame corresponds to a different number of iterations

This is the result:

julia

I love how easy is doing difficult things in R. You can play with the code changing f(z) as well as color palettes. Be ready to get surprised:

library(ggplot2)
library(dplyr)
library(RColorBrewer)
setwd("YOUR WORKING DIRECTORY HERE")
dir.create("output")
setwd("output")
f = function(z,c) exp(z^3)+c
# Grid of complex
z0 <- outer(seq(-2, 2, length.out = 1200),1i*seq(-2, 2, length.out = 1200),'+') %>% c()
opt <-  theme(legend.position="none",
              panel.background = element_rect(fill="white"),
              plot.margin=grid::unit(c(1,1,0,0), "mm"),
              panel.grid=element_blank(),
              axis.ticks=element_blank(),
              axis.title=element_blank(),
              axis.text=element_blank())
for (i in 1:35)
{
  z=z0
  # i iterations of f(z)
  for (k in 1:i) z <- f(z, c=-0.621) df=data.frame(x=Re(z0), y=Im(z0), z=as.vector(exp(-Mod(z)))) %>% na.omit() 
  p=ggplot(df, aes(x=x, y=y, color=z)) + 
    geom_tile() + 
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand=c(0,0))+
    scale_colour_gradientn(colours=brewer.pal(8, "Paired")) + opt
  ggsave(plot=p, file=paste0("plot", stringr::str_pad(i, 4, pad = "0"),".png"), width = 1.2, height = 1.2)
}
# Place the exact path where ImageMagick is installed
system('"C:\\Program Files\\ImageMagick-6.9.3-Q16\\convert.exe" -delay 20 -loop 0 *.png julia.gif')
# cleaning up
file.remove(list.files(pattern=".png"))

Zooming

Drawings, Image Processing, , , ,

You don’t have to be beautiful to turn me on (Kiss, Prince)

I discovered recently how easy is to create GIFs with R using ImageMagick and I feel like a kid with a new toy. To begin this new era of my life as R programmer I have done this:

zooming
First of all, read this article: it explains very well how to start doing GIFs from scratch. The one I have done is inspired in this previous post where I take a set of complex numbers to transform and color it using HSV technique. In this case I use this next transformation: f(z)= -Im(z)+(Re(z)+0.5*Im(z))*1i

Modifying the range of Real and Imaginary parts of complex numbers I obtain the zooming  effect. The code is very simple. Play with it changing the transformation or the animation options. Send me your creations, I would love to see them:

library(dplyr)
library(ggplot2)
dir.create("output")
setwd("output")
id=1 # label tO name plots
for (i in seq(from=320, to=20, length.out = 38)){
z=outer(seq(from = -i, to = i, length.out = 300),1i*seq(from = -i, to = i, length.out = 500),'+') %>% c()
z0=z
for (k in 1:100) z <- -Im(z)+(Re(z)+0.5*Im(z))*1i
df=data.frame(x=Re(z0),
y=Im(z0),
h=(Arg(z)<0)*1+Arg(z)/(2*pi), s=(1+sin(2*pi*log(1+Mod(z))))/2, v=(1+cos(2*pi*log(1+Mod(z))))/2) %>% mutate(col=hsv(h,s,v))
ggplot(df, aes(x, y)) +
geom_tile(fill=df$col)+
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0))+
labs(x=NULL, y=NULL)+
theme(legend.position="none",
panel.background = element_rect(fill="white"),
plot.margin=grid::unit(c(1,1,0,0), "mm"),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
ggsave(file=paste0("plot",stringr::str_pad(id, 4, pad = "0"),".png"), width = 1, height = 1)
id=id+1
}
system('"C:\\Program Files\\ImageMagick-6.9.3-Q16\\convert.exe" -delay 10 -loop 0 -duplicate 1,-2-1 *.png zooming.gif')
# cleaning up
file.remove(list.files(pattern=".png"))

The Hype Bubble Map for Dog Breeds

Going Bananas, Google Trends, , , , ,

In the whole history of the world there is but one thing that money can not buy… to wit the wag of a dog’s tail (Josh Billings)

In this post I combine several things:

  • Simple webscraping to read the list of companion dogs from Wikipedia. I love rvest package to do these things
  • Google Trends queries to download the evolution of searchings of breeds during last 6 months. I use gtrendsR package to do this and works quite well
  • A dinamic Highchart visualization using the awesome highcharter package

The experiment is based on a simple idea: what people search on the Internet is what people do. Can be Google Trends an useful tool to know which breed will become fashionable in the future? To be honest, I don’t really know but I will make my own bet.

What I have done is to extract last 6 months of Google trends of this list of companion breeds. After some simple text mining, I divide the set of names into 5-elements subsets because Google API doesn’t allow searchings with more than 5 items. The result of the query to Google trends is a normalized time series, meaning the 0 – 100 values are relative, not absolute, measures. This is done by taking all of the interest data for your keywords and dividing it by the highest point of interest for that date range. To make all 5-items of results comparable I always include King Charles Spaniel breed in all searchings (as a kind of undercover agent I will use to compare searching levels). The resulting number is my “Level” Y-Axis of the plot. I limit searchings to code=”0-66″ which is restrict results to Animals and pets category. Thanks, Philippe, for your help in this point. I also restrict rearchings To the United States of America.

There are several ways to obtain an aggregated trend indicator of a time series. My choice here was doing a short moving average order=2 to the resulting interest over time obtained from Google. The I divide the weekly variations by the smoothed time series. The trend indicator is the mean of these values. To obtains a robust indicator, I remove outliers of the original time series. This is my X-axis.

I discovered recently a wonderful package called highcharter which allows you to create incredibly cool dynamic visualizations. I love it and I could not resist to use it to do the previous plot with the look and feel of The Economist:

Inspired by Gartner’s Hype Cycle of Emerging Technologies I distinguish two sets of dog breeds:

  • Plateau of Productivity Breeds (succesful breeds with very high level indicator and possitive trend): Golden Retriever, Pomeranian, Chihuahua, Collie and Shih Tzu.
  • Innovation Trigger Breeds (promising dog breeds with very high trend indicator and low level): Mexican Hairless Dog, Keeshond, West Highland White Terrier and German Spitz.

And here comes my prediction. After analyzing the set Innovation Trigger Breeds, my bet is Keeshond will increase its popularity in the nearly future: don’t you think it is lovely?

640px-Little_Puppy_Keeshond
Photo by Terri Brown – Flickr: IMG_4723, CC BY 2.0

Here you have the code:

library(gtrendsR)
library(rvest)
library(dplyr)
library(stringr)
library(forecast)
library(outliers)
library(highcharter)
library(ggplot2)
library(scales)
 
x="https://en.wikipedia.org/wiki/Companion_dog"
read_html(x) %>% 
  html_nodes("ul:nth-child(19)") %>% 
  html_text() %>% 
  strsplit(., "\n") %>% 
  unlist() -> breeds
breeds=iconv(breeds[breeds!= ""], "UTF-8")
usr <- "YOUR GOOGLE ACCOUNT"
psw <- "YOUR GOOGLE PASSWORD"
gconnect(usr, psw)
#Reference
ref="King Charles Spaniel"
#New set
breeds=setdiff(breeds, ref)
#Subsets. Do not worry about warning message
sub.breeds=split(breeds, 1:ceiling(length(breeds)/4))
results=list()
for (i in 1:length(sub.breeds))
{
  res <- gtrends(unlist(union(ref, sub.breeds[i])), 
          start_date = Sys.Date()-180,
          cat="0-66",
          geo="US")
  results[[i]]=res
}
trends=data.frame(name=character(0), level=numeric(0), trend=numeric(0))
for (i in 1:length(results))
{
  df=results[[i]]$trend
  lr=mean(results[[i]]$trend[,3]/results[[1]]$trend[,3])
  for (j in 3:ncol(df))
  {
    s=rm.outlier(df[,j], fill = TRUE)
    t=mean(diff(ma(s, order=2))/ma(s, order=2), na.rm = T)
    l=mean(results[[i]]$trend[,j]/lr)
    trends=rbind(data.frame(name=colnames(df)[j], level=l, trend=t), trends)
  }
}
 
trends %>% 
  group_by(name) %>% 
  summarize(level=mean(level), trend=mean(trend*100)) %>% 
  filter(level>0 & trend > -10 & level<500) %>% 
  na.omit() %>% 
  mutate(name=str_replace_all(name, ".US","")) %>% 
  mutate(name=str_replace_all(name ,"[[:punct:]]"," ")) %>% 
  rename(
    x = trend,
    y = level
  ) -> trends
trends$y=(trends$y/max(trends$y))*100
#Dinamic chart as The Economist
highchart() %>% 
  hc_title(text = "The Hype Bubble Map for Dog Breeds") %>%
  hc_subtitle(text = "According Last 6 Months of Google Searchings") %>% 
  hc_xAxis(title = list(text = "Trend"), labels = list(format = "{value}%")) %>% 
  hc_yAxis(title = list(text = "Level")) %>% 
  hc_add_theme(hc_theme_economist()) %>%
  hc_add_series(data = list.parse3(trends), type = "bubble", showInLegend=FALSE, maxSize=40) %>% 
  hc_tooltip(formatter = JS("function(){
                            return ('<b>Trend: </b>' + Highcharts.numberFormat(this.x, 2)+'%' + '<br><b>Level: </b>' + Highcharts.numberFormat(this.y, 2) + '<br><b>Breed: </b>' + this.point.name)
                            }"))

The Coaster Maker by Shiny

Drawings, Shiny, , ,

The word you invented is well formed and could be used in the Italian language (The Accademia della Crusca regarding to the word “Petaloso”, recently invented by an eight-year-old boy)

Are you tired of your old coasters? Do you like to make things by your own? Do you have a PC and a printer at home? If you answered yes to all these questions, just follow these simple instructions:

  • Install R and RStudio in your PC
  • Open RStudio and create a new Shiny Web App multiple file (ui.R/server.R)
  • Substitute sample code of each file by the code below
  • Press Run App
  • Press buttom Get your coaster! until you obtain a image you like
  • Print the image
  • Cut out the image
  • Place on the coaster your favorite drinking

These are some examples:

This is the code of ui.R

#
# This is the user-interface definition of a Shiny web application. You can
# run the application by clicking 'Run App' above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#
library(shiny)
shinyUI(fluidPage(
  titlePanel("The coaster maker"),
  sidebarLayout(
    sidebarPanel(
      #helpText(),
 
      # adding the new div tag to the sidebar
      tags$div(class="header", checked=NA,
               tags$p("This coasters are generated by hypocycloid curves.The curve is formed by the locus of a point,
                      attached to a circle, that rolls on the inside of another circle.
                      In the curve's equation the first part denotes the relative position between the two circles,
                      the second part denotes the rotation of the rolling circle.")),
      tags$div(class="header", checked=NA,
               HTML(
"More info <a href=\"http://www.2dcurves.com/roulette/rouletteh.html#rhodon\">here</a>")),
      actionButton('rerun','Get your coaster!')
    ),
    mainPanel(
      plotOutput("HarmPlot")
    )
  )
))

This is the code of server.R

# This is the server logic of a Shiny web application. You can run the
# application by clicking 'Run App' above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#
library(shiny)
library(ggplot2)
CreateDS = function ()
{
  t=seq(-31*pi, 31*pi, 0.002)
  a=sample(seq(from=1/31, to=29/31, by=2/31), 1)
  b=runif(1, min = 1, max = 3)
  data.frame(x=(1-a)*cos(a*t)+a*b*cos((1-a)*t), y=(1-a)*sin(a*t)-a*b*sin((1-a)*t))
}
shinyServer(function(input, output) {
  dat<-reactive({if (input$rerun) dat=CreateDS() else dat=CreateDS()})
    output$HarmPlot<-renderPlot({
    ggplot(dat())+
      geom_point(data=data.frame(x=0,y=0), aes(x,y), color=rgb(rbeta(1, .5, .5), rbeta(1, .5, .5), rbeta(1, .5, .5)) , shape=19, fill="yellow", size=220)+
      geom_polygon(aes(x, y), fill=rgb(rbeta(1, 2, 2), rbeta(1, 2, 2), rbeta(1, 2, 2))) +
      theme(legend.position="none",
            panel.background = element_rect(fill="white"),
            panel.grid=element_blank(),
            axis.ticks=element_blank(),
            axis.title=element_blank(),
            axis.text=element_blank())
  }, height = 500, width = 500)
})

Sunflowers

Drawings, , , , ,

The world is full of wonderful things, like sunflowers (Machanguito, my islander friend)

Sunflower seeds are arranged following a mathematical pattern where golden ratio plays a starring role. There are tons of web sites explaining this amazing fact. In general, the arrangement of leaves on a plant stem are ruled by spirals. This fact is called phyllotaxis, and I did this experiment about it some time ago. Voronoi tessellation originated by points arranged according the golden angle spiral give rise to this sunflowers:sunflowers

I know this drawing will like to my friend Machanguito because he loves sunflowers. He also loves dancing, chocolate cookies, music and swimming in the sea. Machanguito loves life, it is just that simple. He is also a big defender of renewable energy and writes down his thoughts on recycled papers. You can follow his adventures here.

This is the code:

library(deldir)
library(ggplot2)
library(dplyr)
opt = theme(legend.position  = "none",
            panel.background = element_rect(fill="red4"),
            axis.ticks       = element_blank(),
            panel.grid       = element_blank(),
            axis.title       = element_blank(),
            axis.text        = element_blank())
CreateSunFlower <- function(nob=500, dx=0, dy=0) {   data.frame(r=sqrt(1:nob), t=(1:nob)*(3-sqrt(5))*pi) %>%
    mutate(x=r*cos(t)+dx, y=r*sin(t)+dy)
}
g=seq(from=0, by = 45, length.out = 4)
jitter(g, amount=2) %>%
  expand.grid(jitter(g, amount=2)) %>%
  apply(1, function(x) CreateSunFlower(nob=round(jitter(220, factor=15)), dx=x[1], dy=x[2])) %>%
  do.call("rbind", .) %>% deldir() %>% .$dirsgs -> sunflowers
ggplot(sunflowers) +
  geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), color="greenyellow") +
  scale_x_continuous(expand=c(0,0))+
  scale_y_continuous(expand=c(0,0))+
  opt

A Silky Drawing and a Tiny Experiment

Drawings,

It is a capital mistake to theorize before one has data (Sherlock Holmes, A Scandal in Bohemia)

One of my favorite entertainments is drawing things: crazy curves, imaginary flowers, celestial bodies, fractalic acacias … but sometimes I wonder myself if these drawings result interesting to whom arrive to my blog. One way to define interesting could be wanting to reproduce the drawing. I know some people do it because they sometimes share with me their creations. So, how many people appreciate the code I write? I manage some a priori for this number (which I will maintain for myself) but I want to refine my estimation with the next experiment. I have done this drawing, which shows that simple mathematics can produce very nice patterns:

silky

To estimate how many people is really interested in this plot, at the end of the post I will publish all the code except for a line. If you want the line, you will have to ask it to me. How? It is very easy: you will have to send me a direct message in Twitter. If you don’t follow me, do it here and I will follow you back. If you already follow me but I don’t, tweet something mentioning me and I will follow you back. Then you will be able to send me the direct message. If you prefer, you can send me an email. You can find my email address here.

I know this experiment can be quite biased, but I am also pretty sure that the resulting estimation will be much better than the one I manage nowadays. This is the kidnapped code:

library(magrittr)
library(ggplot2)
opt = theme(legend.position  = "none",
panel.background = element_rect(fill="violetred4"),
axis.ticks       = element_blank(),
panel.grid       = element_blank(),
axis.title       = element_blank(),
axis.text        = element_blank())
seq(from=-10, to=10, by = 0.05) %>%
expand.grid(x=., y=.) %>%
#HERE COMES THE KIDNAPPED LINE
geom_point(alpha=.1, shape=20, size=1, color="white") + opt

The Unbereable Insolence of Prime Numbers or (Playing to be Ulam)

Curiosities, , ,

So rock me mama like a wagon wheel, rock me mama anyway you feel (Wagon Wheel, Old Crow Medicine Show)

This is the third iteration of Hilbert curve. I placed points in its corners. Since the curve has beginning and ending, I labeled each vertex with the order it occupies:hilbert_primes3_1Dark green vertex are those labeled with prime numbers and light ones with non-prime. This is the sixth iteration colored as I described before (I removed lines and labels):hilbert_primes6_2

Previous plot has 4.096 points. There are 564 primes lower than 4.096. What If I color 564 points randomly instead coloring primes? This is an example:
hilbert_primes6_2rand
Do you see any difference? I do. Let me place both images together (on the left, the one with primes colored):
hilbert_primes6_3

The dark points are much more ordered in the first plot. The second one is more noisy. This is my particular tribute to Stanislaw Ulam and its spiral: one of the most amazing fruits of boredom in the history of mathematics.

This is the code:

library(reshape2)
library(dplyr)
library(ggplot2)
library(pracma)
opt=theme(legend.position="none",
          panel.background = element_rect(fill="white"),
          panel.grid=element_blank(),
          axis.ticks=element_blank(),
          axis.title=element_blank(),
          axis.text=element_blank())
hilbert = function(m,n,r) {
  for (i in 1:n)
  {
    tmp=cbind(t(m), m+nrow(m)^2)
    m=rbind(tmp, (2*nrow(m))^r-tmp[nrow(m):1,]+1)
  }
  melt(m) %>% plyr::rename(c("Var1" = "x", "Var2" = "y", "value"="order")) %>% arrange(order)}
iter=3 #Number of iterations
df=hilbert(m=matrix(1), n=iter, r=2)
subprimes=primes(nrow(df))
df %>%  mutate(prime=order %in% subprimes,
               random=sample(x=c(TRUE, FALSE), size=nrow(df), prob=c(length(subprimes),(nrow(df)-length(subprimes))), replace = TRUE)) -> df
#Labeled (primes colored)
ggplot(df, aes(x, y, colour=prime)) +
  geom_path(color="gray75", size=3)+
  geom_point(size=28)+
  scale_colour_manual(values = c("olivedrab1", "olivedrab"))+
  scale_x_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
  scale_y_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
  geom_text(aes(label=order), size=8, color="white")+
  opt
#Non labeled (primes colored)
ggplot(df, aes(x, y, colour=prime)) +
  geom_point(size=5)+
  scale_colour_manual(values = c("olivedrab1", "olivedrab"))+
  scale_x_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
  scale_y_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
  opt
#Non labeled (random colored)
ggplot(df, aes(x, y, colour=random)) +
  geom_point(size=5)+
  scale_colour_manual(values = c("olivedrab1", "olivedrab"))+
  scale_x_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
  scale_y_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
  opt