Twitter incessantly produces copious amount of data. The locations of tweets can help with some interesting questions. One that comes to mind, and one that I plan to do when the time is right is- What part of the world is interested in The Champions League final vs The World Cup final. The rationale for this is the amount of debate currently happening on this topic.
Basically, this post will answer "where in the world are people searching for [something]?" Also, all explanation will be in the comments itself.

Code

library(ROAuth)
library(twitteR)
library(ggplot2)
library(maps)
library(dismo)

# Check to see if R is connected to twitter
registerTwitterOAuth(twitCred)

# searchString parameter for Twitter API
key <- "#UCLfinal"

# requesting Twitter API
tag <- searchTwitter(key, n = 2000, lang= "en") 

# Tweets data frame
# At this stage it is quite possible to get rid of all the non-geotagged tweets.
# However, a very very small portion of users geotag tweets. Therefore, another approach
# is used here. In the next step, the location in their profile description will
# be extracted.
df <- do.call("rbind", lapply(tag, as.data.frame))

# User data frame
userInfo <- do.call("rbind", lapply(lookupUsers(df$screenName), as.data.frame))   

# geocoding all users with some sort of location identification
# also creating "Interpreted Place"
# All locations with invalid location will be dropped after this step
# Package dismo used here
# Although using oneRecord=T decreases the size, it produces more
# reliable location data frame
locations <- geocode(userInfo$location, progress="text", oneRecord=T)

# getting rid of all rows with na
locations <- locations[complete.cases(locations),]

# also getting rid of all tweets that only have country name as locations
# For example, it is good to avoid the center of Australia (which is very very sparse) 
# showing massive number of tweets.
# The easiest way to do this is to get rid of all rows that do not have a comma
locations <- locations[grep("\\,",locations$interpretedPlace),]

# Map of the world
# It is also possible to do the same with country/places maps.
# However, it is necessary to make sure that the coordinates are correct
result <- ggplot(map_data("world")) + geom_path(aes(x = long, y = lat, group = group))

# Adding Tweet locations
result <- result + geom_point(data = locations, aes(x = longitude, y = latitude),
                              color = "red", alpha = .2, size = 3)
result <- result + ggtitle(key) + theme_minimal() + theme(axis.text=element_blank(),
                                                          axis.ticks=element_blank(),
                                                          axis.title=element_blank())  

# Time Stamping the file name
filename <- paste(format(Sys.Date(),"%d%m%y"),format(Sys.time(), "%H%M%S"),".png",sep="")

# Saving the file
ggsave(filename, units="in", width=8.15, height=5.20, dpi=300)

After coupling the result with the wordcloud code (Click here to go to that blog post) and Photoshop, these are the products.

#ChampionshipPlayOffFinal

This is an austere implementation of mining tweets using R. I will work on making it better sometime in the future (once I engulf R data structures more).

Connecting to Twitter API

The first step was to get consumer key and consumer secret from Twitter. This is used for authentication purposes. To get these, I just created a new app on Twitter.
Apparently, on Windows system the authentication requires an extra step. Thanks to this blog post for clarifying that. See below.

library("ROAuth")
library("twitteR")

#necessary step for Windows
download.file(url="http://curl.haxx.se/ca/cacert.pem", destfile="cacert.pem")

#to get your consumerKey and consumerSecret see the twitteR documentation for instructions
credentials <- OAuthFactory$new
                        (consumerKey='#################...',
                         consumerSecret='#################...',
                         requestURL='https://api.twitter.com/oauth/request_token',
                         accessURL='https://api.twitter.com/oauth/access_token',
                         authURL='https://api.twitter.com/oauth/authorize')

credentials$handshake(cainfo="cacert.pem")
save(credentials, file="twitter authentication.Rdata")
registerTwitterOAuth(credentials)

After following the link provided, allowing access to the Twitter app and copying the security key produced back in R, everything was all set. To check:

> registerTwitterOAuth(credentials)
[1] TRUE

Tweet Cloud

library(tm)
library(wordcloud)

tag<- searchTwitter("#AFC", n=100, cainfo="cacert.pem")
df <- do.call("rbind", lapply(tag, as.data.frame))
twitterCorpus <- Corpus(VectorSource(df$text))
twitterCorpus <- tm_map(twitterCorpus, tolower)
# remove punctuation
twitterCorpus <- tm_map(twitterCorpus, removePunctuation)
# remove numbers
twitterCorpus <- tm_map(twitterCorpus, removeNumbers)
tdm <- TermDocumentMatrix(twitterCorpus, control = list(minWordLength = 1))
m <- as.matrix(tdm)
# calculate the frequency of words
v <- sort(rowSums(m), decreasing=TRUE)
words <- names(v)
d <- data.frame(word=words, freq=v)
wordcloud(d$word, d$freq, min.freq=3)

Result

#AFC

Although it generates the cloud, at this point, it is not effective enough for analysis. One way to make it better is to get rid of common English words such as "and", "but", etc. Secondly, getting rid of Unicode and non-English languages would be an improvement also. And, there is MORE that needs to be done.
Hopefully soon.

Soon

I just learnt that the common words mentioned above are referred as stop-words. And, it is pretty easy to handle them. After removing them and adding some colors the clouds look more beautiful. I also got rid of all words with substring "http" (maybe not in the most efficient way.. I will have to look into it).

# First of all make sure the connection to twitter is authenticated.
# see file twiiterconnection.R
# 
# if that has been done previously and this is a new session, open .Rdata
# that has stored credentials
# 
#if registerTwitterOAuth(credentials) is TRUE, you are good to go

library(tm)
library(wordcloud)
library(ROAuth)
library(twitteR)

# check connection
registerTwitterOAuth(credentials)

tag<- searchTwitter("#MH370", n=700, lang= "en", cainfo="cacert.pem")
df <- do.call("rbind", lapply(tag, as.data.frame))

twitterCorpus <- Corpus(VectorSource(df$text))
tdm <- TermDocumentMatrix(
  twitterCorpus, control = list(minWordLength = 1,
                                removePunctuation = TRUE, 
                                removeNumbers = TRUE,
                                stopwords = TRUE,
                                tolower=TRUE))

m <- as.matrix(tdm)

# frequency of words
frequency <- sort(rowSums(m), decreasing=TRUE)
words <- names(frequency)

# get rid of words that contains "http"
# hopefully there is a better way to do it
words[which(grepl("http",names(frequency)))] <- ""

d <- data.frame(word=words, freq=frequency)

#save the image in png format
png("mh370.png", width=12, height=12, units="in", res=300)
wordcloud(d$word, d$freq,scale=c(10,.4),min.freq=10,
          max.words=Inf, random.order=FALSE, rot.per=.3, colors=brewer.pal(8, "Dark2"))
dev.off()

Result

#MH370

#AFC

wordcloud(d$word, d$freq,scale=c(8,.1),min.freq=7,
          max.words=Inf, random.order=FALSE, 
          rot.per=.3, vfont=c("gothic english","plain"),
          colors=brewer.pal(8, "Dark2"))

After changing the font, inverting the color and messing around in Photoshop:

I have always loved symmetry in photography. The picture above was created in Photoshop simply by coping and rotating certain parts. See SchlockySymmetry for other symmetrical photos. Recently, I stumbled upon an article in R-bloggers that used abind and biOps for some simple image manipulation. I thought it would be fun to perform some of the tricks mentioned in the article to achieve the picture above.

R Code:

#Image processing

#Read Original Image
x <- readJpeg("Boudha.jpg")

plot(x)
dev.copy(png,filename="original.png",width=2048,height=1206);
dev.off ();

#copy pixels from left to right and mirror
y <- imagedata(abind(x[,1:(ncol(x)/2),], x[,(ncol(x)/2):1,] , along=2))

plot(y)
dev.copy(png,filename="halfmirror.png",width=2048,height=1206);
dev.off ();

#copy all pixels and mirror rows
z <- imagedata(abind(y[1:(nrow(x)),,], y[(nrow(x)):1,,] , along=1))

plot(z)
dev.copy(png,filename="final.png",width=2048,height=1206);
dev.off ();

Result

original.png

halfmirror.png

final.png

The Zorganian Republic has some very strange customs. Couples only wish to have female children as only females can inherit the family's wealth, so if they have a male child they keep having more children until they have a girl. If they have a girl, they stop having children. What is the ratio of girls to boys in Zorgania?

The ratio of girls to boys in Zorgania is 1:1. This might be a little counter-intuitive at first. Here are some ways of tackling this problem.

1. Monte Carlo Simulation:

Although, Monte Carlo simulation does not necessarily show why the result is 1:1, it is appropriate because of the very counter-intuitive nature of the problem. At the very least, it helps us see that the result is indeed 1:1. Therefore, this is a good start.

The following R code estimates the probability of a child being a boy in Zorgania. 

couples <- 100000
boycount <- 0

for (i in 1:couples){
    # 0: boy
    while (sample(c(0,1),1) == 0) {
    boycount=boycount+1
  }
}
probability <- boycount/(couples+boycount)

Result:

2. Understanding the question better:

Here is another question: What is the probability of getting a tail in a fair coin toss, if all seven previous tosses resulted in heads? Since coin flips are independent events, the probability is still going to be 0.5. Similarly in this case, the child births are independent. It does not matter if the couples stop giving birth after they have a baby-girl. The expected value is unchanged. 

For example, consider five couples: C1, C2, C3, C4 and C5.

If B-> Boy and G-> Girl. Using R's sample().

For C1:  {B, G}

For C2: {G}

For C3: {B, B, G}

For C4: {B, G}

For C5: {G}

Now, ignore the couples and only consider the children. The children are {B, G, G, B, B, G, B, G, G}. The only thing happening here is simply the generation of  either a B or a G with equal probability for each generation. At this point, it is quite obvious that the part that has to do with "couple...." in the question is to mislead and confuse similar to the "previous seven tosses.." example that I mentioned in the beginning of 2. 

3.  Counting:

If we start with 512 couples (hence 512 first borns), half of them are going to have a girl as their first. Those couples will stop having children. Among, the other half couples who had a son as their first child, half of them are going to have a girl as their second child and so on. At every step there is an equal numbers of boys and girls. Therefore, the expected ratio is 1:1.

This is the naive/brute-force implementation of the Mandelbrot Set plotting. I just followed the algorithm.  

# Plotting the Mandelbrot Set 

# length of sequence for real and imaginary parts of complex numbers
length <- 1000

# sequences for real and imaginary parts
real = seq(-1.8,0.6, len=length)
imaginary = seq(-1.2,1.2, len=length)


result <- matrix(nrow = length, ncol = length)

for (i in 1:length)
{
  for (j in 1:length)
  {
    result[i,j]=inmandelbrotset(complex(real = real[i], imaginary = imaginary[j]))       
  }
}

image(result, axes=FALSE)

# function that checks if a point E mandelbrot set
inmandelbrotset <- function(c)
{
  dwell.limit <- 2048
  z <- 0  
  for (i in 1:dwell.limit)
  { 
    z <- z ** 2 + c
    if (Mod(z) > 2)
    {
      return(FALSE)
    }
  }  
  return(TRUE)
}

Adding colors:

We now have a Boolean matrix that records if a point is in the Mandelbrot Set. Since the matrix can only have two values : true or false, thus far, we have only been able to plot read and white images. The next step is to add colors such that we get more information on when a particular point escapes the radius of 2. Again, this is the naive/brute force way of doing it.

# Mandelbrot Plotting with colors 
length <- 1000
real = seq(-2.0,2.0, len=length)
imaginary = seq(-2.0,2.0, len=length)
result <- matrix(nrow = length, ncol = length)
dwell.limit <- 512

for (i in 1:length)
{
  for (j in 1:length)
  {  
    z <- 0
    c <-complex(real = real[i], imaginary = imaginary[j])
    for (k in 1:dwell.limit)
    { 
      
      z <- z ** 2 + c
      if (Mod(z) > 2)
      {      
        result[i,j]=k
        break
      }
    }     
    
  }
}
set.seed(2)
image(result,breaks=0:dwell.limit
      ,col=c(1,sample(terrain.colors
                      (dwell.limit-1,alpha = .8))),asp=1,ax=F)

and, just for the heck of it..

ASCII Mandelbrot Set using R (naive)

s <- seq(-1.7,1.2, by =.1)
a <- ""
for (i in 1:length(s))
{  
  for (j in 1:length(s))
  {
   a<-cat(a,inmandelbrotset(complex(r = s[j], i = s[i])))     
  }  
  a <- cat(a,"\n")
}

Achieved by returning a " " or "#" instead of FALSE or TRUE from function "inmandelbrotset".

A better algorithm

Utilizing R's easy to use lists in implementation:

# more efficient algorithm to plot the Mandelbrot set 

sequence <- seq(-2,2,len=1000)
dwell.limit <- 200

# matrix of points to be iterated 
complex.matrix <- t((sapply(sequence,function(x)x+1i*sequence)))
in.mandelbrot.index <- 1:length(complex.matrix)
iter=z=array(0,dim(complex.matrix)) 
           
for(i in 1:dwell.limit){ 
  # complex quadratic polynomial function for all points
  z[in.mandelbrot.index]=complex.matrix[in.mandelbrot.index]+z[in.mandelbrot.index]^2 
  # boolean matrix
  result=Mod(z[in.mandelbrot.index])<=2  
  # if result is false, store the iteration 
  iter[in.mandelbrot.index[!result]]=i
  # save all the index where points are still in the mandelbrot
  in.mandelbrot.index=in.mandelbrot.index[result]
}
set.seed(19)
image(iter,main=paste("Iterations: ", i, sep=" "), breaks=0:dwell.limit
      ,col=c(1,sample(rainbow
                      (dwell.limit-1,alpha = .8))),ax=F, asp=1)

Plotting the Julia set

A little modification to the code above (red and white Mandelbrot) produces Julia Sets. The idea here is to set a constant C and send Z to the function instead of C.

c <- complex(real=-0.1,imaginary=0.651) 
label <- toString(c)
injulia <- function(z)
{
  dwell.limit <- 128
  
  for (i in 1:dwell.limit)
  { 
    z <- z ** 2 + c
    if (Mod(z) > 2)
    {
      return(FALSE)
    }
  }  
  return(TRUE)
}

Adding colors:

This is achieved by following the same process as above.

Sierpinski Gasket using Chaos game

#### Chaos game for generation of Sierpinski Gasket
# 1. Take 3 points in a plane to form a triangle, you need not draw it.
# 2. Randomly select any point inside the triangle and consider that your current position.
# 3. Randomly select any one of the 3 vertex points.
# 4. Move half the distance from your current position to the selected vertex.
# 5. Plot the current position.
# 6. Repeat from step 3

plot.new()
iterations <- 2000
vertices <- matrix(c(0,0,0.5,1,1,0),3,2, byrow=T)
current.point <- c(0.5,0.5)
random.vertex <- sample(1:3,iterations,replace=T)
plot.result = matrix(nrow=iterations,ncol=2)
for (i in 1:iterations){
 current.point <- (current.point+vertices[random.vertex[i],])/2
 plot.result[i,] <- current.point
}
points(plot.result,pch = 46)

Adding colors:

points(plot.result,pch = 46,col=c(13,3,41)[random.vertex])

I used R for two of my classes in college: Probability and Statistics and Stochastic Probabilities. During that time, I was in reverence with the way it handled lists. However, I did not use it much after getting done with those classes. I have decided to relearn it.

In this post, I will start off with some simple programs. After ascertaining some comfort with the syntax and semantics of R, I will try some complex ones. Also, since I am in the process of learning myself, if you decide to use some of these codes, just keep in mind that there might be several better ways of doing the same thing.

1. Dice Rolls: Probability of getting a side f number of times in n dice rolls

# DiceRolls.R
# Probability of getting a side f number of times in n dice rolls

#inputs
side <- 4
f <- 2
n <- 5

# Binomial Distribution
# using R's function for binomial distribution 
# where P(X=k)
 probability <- dbinom(frequency, n, 1/6)
 probability

# Monte Carlo Method
rep <- 10000 # Number of repetitions
list <- numeric(rep) # Initialize vector
for (i in 1:rep) {
trial <- sample(1:6, n, replace=TRUE)
#1 if number of times a side appears is f
success <- if ((sum(trial==side))==f) 1 else 0
list[i] <- success 
}
mean(list)

>  probability
[1] 0.160751
> mean(list)
[1] 0.1647

2. One dimentional random walks

# Five random walks in one dimension starting at 0

# creating a color string vector
color <- c("red","green","blue","black","darkgoldenrod1")
# cleaning the plot; this is necessary because we are using par(new=TRUE)
plot.new()
for (i in 1:5) {
x0<-0
I<-1
T <-1000
x=c(x0,rnorm(T-1)) # T-1 randomly generated normal distributions
y=cumsum(x)
par(new=TRUE) # so that the previous graphs are not deleted
plot(y,type="l",col=color[i], main="Five random walks in one dimension starting at 0", xlab="", ylim=c(-60, 60))
}

3. Two dimensional random walk

I stumbled upon a presentation slideshow that discusses how R can be used more efficiently. They perform analysis on several implementations that simulate two dimensional random walk. I started out with what they call "naive implementation". According to them, this is how a person who codes in C, C++ or Java would think of doing it. I totally agree. This exercise is great because it helps understand the positives of R in comparison to other languages that I am comfortable with. Link to the aforementioned presentation slideshow. 

Naive implementation:

# naive implementation (loops are slow in R)
naive <- function(n){

# numeric creates a double-precision vector of the specified length with each element equal to 0
xpos = ypos = numeric(n)

# will be used to randomly generate direction
xdir = c(TRUE,FALSE)
nextdir = c(1,-1)

# xpos and ypos are already 0 (see comment on numeric)
for(i in 2:n){

 if(sample(xdir,1)){
  xpos[i]=xpos[i-1]+sample(nextdir,1)
  ypos[i]=ypos[i-1]
 }
 else{
  xpos[i]=xpos[i-1]
  ypos[i]=ypos[i-1]+sample(nextdir,1)
 }
}
plot(xpos,ypos,type="l",col="blue",main="Two Dimensional Random Walk", xlab="x", ylab="y")
# starting point
points(0,0,pch = 19,col="red")
}

> system.time(naive(50000))
   user  system elapsed 
   0.84    0.17    1.02 

> system.time(naive(100000))
   user  system elapsed 
   1.70    0.32    2.03 

Vector implementation:

This implementation gets rid of the for loop. It is similar to the code in one dimensional random walk code above.

# Vectorization Process
factorization = function(n){ 
 steps = sample(c(-1,1), n-1, replace = T)
 xdir=sample(c(TRUE,FALSE), n-1,replace=T)
 xpos=c(0,cumsum(ifelse(xdir,steps,0)))
 ypos=c(0,cumsum(ifelse(xdir,0,steps)))
 plot(xpos,ypos,col="green",type="l",main="Two Dimensional Random Walk", xlab="x", ylab="y")
 points(0,0,pch = 19,col="red")
}

> system.time(vectorization(50000))
   user  system elapsed 
   0.05    0.19    0.25 

> system.time(vectorization(100000))
   user  system elapsed 
   0.09    0.31    0.40 

4. Birthday Probability

Click here to go to my post on Birthday Probability using Java. Using R, I was really surprised how easy it was.

How a java programmer would have done it in R

# Birthday Probability
birthday <- function(numofpeople, numofrooms){
  count <- 0
  for (i in 1:numofrooms){
    bday <- sample(1:365,numofpeople, replace=T)
    if (any(duplicated(bday)))
      count<-count+1;
  }  
probability = count/numofrooms
probability
}

How a R programmer would have done it (I think)

# Birthday Probability
birthday <- function(numofpeople, numofrooms){ 
  list <- numeric(numofrooms)
  
  for (i in 1:numofrooms){
    bday <- sample(1:365,numofpeople, replace=T)
    list[i] <- if (any(duplicated(bday))) 1 else 0   
  }   
mean(list)
}

Result:

> source("birthday.r")
> birthday(23,5000)
[1] 0.5044

5. Mandelbrot Set Area

Click here to go to my previous blog post on the area of the Mandelbrot Set using Java and Javascript. 

monte.Carlo <- 5000
x <- runif(monte.Carlo, -2, 2)
y <- runif(monte.Carlo, -2, 2)
list <- numeric(monte.Carlo)

for (j in 1:monte.Carlo){
  list[j] <- if (inmandelbrotset(complex(real = x[j], imaginary = y[j]))) 1 else 0
}
area<-mean(list)*16

# function that checks if a point E mandelbrot set
inmandelbrotset <- function(c)
{
  dwell.limit <- 2048
  z <- 0  
  for (i in 1:dwell.limit)
  { 
    z <- z ** 2 + c
    if (Mod(z) > 2)
    {
      return(FALSE)
    }
  }  
  return(TRUE)
}

Result:

> system.time(source("mandelbrot.r"))
   user  system elapsed 
   1.91    0.00    1.90 
> area
[1] 1.5168