Friday, August 1, 2014

UPDATE: 100 Prisoners, 100 lines of code, a game theory example in R

I  was looking for some Game-Theory post in R-bloggers and found this very cool post,  obviously I decided to read the post, to look at the code and implement a version of it.

The idea of the post was very cool for me that I thought it would be a nice idea to get more "juice" from it.

All the credit for this post, is for its original author.

The code is available here.

As Matt Asher describes in his original post:

# ***************
# ***************

In this game, the warden places 100 numbers in 100 boxes, at random with equal probability that any number will be in any box. Each convict is assigned a number. One by one they enter the room with the boxes, and try to find their corresponding number. They can open up to 50 different boxes. Once they either find their number or fail, they move on to a different room and all of the boxes are returned to exactly how they were before the prisoner entered the room.

The prisoners can communicate with each other before the game begins, but as soon as it starts they have no way to signal to each other. The warden is requiring that all 100 prisoners find their numbers, otherwise she will force them to listen to hundreds of hours of non-stop, loud rock musician interviews. Can they avoid this fate?

The first thing you might notice is that if every prisoner opens 50 boxes at random, they will have a 0.5 probability of finding their number.

There's actually a strategy that can guarantee all the prisoners to get their numbers.

For this post, we will call that strategy, model, and the strategy with no strategy at all, we will call it random.model.

The model strategy consist in that at the beginning, each prisoner will open the box based on their number, for example, the prisoner 42 will open the box at the position 42, after open it, say it contains the number 12, and given that 42 is not equal to 12, he will have to open the box at position 12. This will happen until he had opened 50 boxes or found his number.

The random.model strategy consist in no strategy at all, only opening boxes at random until he had opened 50 boxes or found his number.

Every time a prisoner found his number we will call this event a success.

Our interest is to find the strategy that maximizes the success rate. We will run each simulation 1000 times.

# ***************
# ***************

Let's take a look at the performance of each model individually with this histogram:

We can see that the random.model behaves as we were expecting assuming that by chance we have a 0.5 chance of getting our number by selecting 50 random boxes, On the other hand, the other strategy has a very different distribution, there's a gap between 51 and 99 with 0 of success rate!, any suggestion to explain how is this possible?

Another way to plot the distribution of the success rates obtained from each model is by looking at a density plot:

This histogram is telling the same story than the previous plots, but in the same plot:

And as usual, I always left the best plot for the grand finale. 

What we can say about this plot?, before going there, I will continue using Matt's explanation:

The theoretical success rate of the model strategy is about 31%

One way to look at the distribution of numbers in boxes is to see it as a permutation of the numbers from 1 to 100. Each permutation can be partitioned into what are called cycles. A cycle works like this: pick any number in your permutation. Let's say it's 23. Then you look at the number the 23rd place (ie the number in the 23rd box, counting from the left). If that number is 16, you look at the number in the 16th place. If that number is 87, go open box number 87 and follow that number. Eventually, the box you open up will have the number that brings you back to where you started, completing the cycle. Different permutations have different cycles.

The key for the prisoner is that by starting with the box that is the same place from the left as his number, and by following the numbers in the boxes, the prisoner guarantees that if he is in a cycle of length less than 50, he will eventually open the box with his number in it, which would complete the cycle he began. One way to envision cycles of different lengths is to think about the extreme cases. If a particular permutation shifted every single number over one to the left (and wrapped number 1 onto the end), you would have a single cycle of length 100. Box 1 would contain number 2, box 2 number 3 and so on. On the other hand, if a permutation flipped every pair of consecutive numbers, you would have 50 cycles, each of length 2: box 1 would have number 2, box 2 would have number 1. Of course if your permutation doesn't change anything you have 100 cycles of length 1.

As you can see from the histogram of the model strategy, you can never have between 50 and 100 winning prisoners. Anytime you have a single cycle of length greater than 50, for example 55, then all 55 prisoners who start on that cycle will fail to find their number. If no cycles are longer than 50, everyone wins.

Now, looking at the boxplot, OK we would say, using the model strategy we have like a 30% percent of chance that every one of the prisoners will get their number, but on the other hand, the random strategy guarantees that at least every prisoner will have 50% chance of getting his number, and therefore win.

So, the model strategy is like an all of nothing strategy, you can win the lottery or die trying. This is observable looking at the large standard deviation in the boxplot, on the other hand, the random.model strategy is very consistent with low dispersion values. This means that if the prisoners decides to use the model strategy, they will have to face tree possible outcomes, or everybody wins, or nobody does, or only a few does. On the other hand, with the random.model strategy, they will expect that at least 50 prisoners will win, thinking as a prisoner as a logical individual, I think this is the best strategy he can use.

Write your opinions.

Thanks to Matt Asher for the idea


Monday, July 7, 2014

What are the names of the school principals in Mexico?, If your name is Maria, probably this post will interest you. Trends and cool plots from the national education census of Mexico in 2013

I will start this post with a disclaimer:

The main intention of the post is to show how is the distribution of the school principal names in Mexico, for example, to show basic trends regarding about what is the most common nation-wide first name and so on, also to show trends delimited by state and regions.

These trends in data would answer questions such:

1. Are the most common first names distributed equally among the states?
2. Does the states sharing the same region also share the same "naming" behavior?

Additionally, this post includes cool wordclouds.

Finally and the last part of my disclaimer is that, I am really concerned about the privacy of the persons involved. I am not by any sense promoting the exploitation of this personal data, if you decide to download the dataset, I would really ask you to study it and to generate information that is beneficial, do not join the Dark side.



The database is located here
The R code can be downloaded here
Additional data can be downloaded here

All the results were computed exploring 202,118 schools across the 32 states of Mexico from the 2013 census


Here is the wordcloud of names (by name, I am referring to first name only), it can be concluded that MARIA is by far the most common first name of a school principal in all Mexican schools, followed by JOSE and then by JUAN

The following wordcloud includes every word in the responsible_name column (this includes, first name, last names). Now the plot shows that besides the common first name of MARIA, also the last names of HERNANDEZ, MARTINEZ and GARCIA are very common.


Looking at this barplot, the name MARIA is by far the most common name of the Mexican school's principals, with a frequency ~25,000. The next most popular name is JOSE with a frequency of ~7,500

Looking at the same data, just adjusted to represent the % of each name inside the pool of first names we have that MARIA occupy ~11% of the names pool.


 With this heatmap, my intention is to show the distribution of the top 20 most common first names across all the Mexican states

It can be concluded that there is a small cluster of states which keep the most number of principals named MARIA(but no so fast!, some states, for example Mexico and Distrito Federal are very populated, so I will reduce this effect in the following plot). In summary the message of this plot is the distribution of frequency of the top 20 most frequent first-names across the country.


For me, a young data-science-padawan, this is my favorite analysis: "hunting down the trends".

The setup of the experiment is very simple, map the top 1,000 most frequent nation-wide names across each state to create a 32 x 1000 matrix (32 states and 1,000 most nation-wide frequent names).

With this matrix, normalize the values by diving each row by the sum of it  (this will minimize the effect of the populated states vs the non populated while maintaining the proportion of the name frequencies per state). Then I just computed a distance matrix and plotted it as a heatmap.

What I can conclude with this plot is that, there are clusters of states that seems to maintain a geographical preference to be clustered within its region, this would be concluded that it is likely that states sharing the same regions would be more likely to share the "naming" trends due to some cultural factors (like the cluster that includes Chihuahua, Sonora and Sinaloa). But this effect is not present in all the clusters.

All images can be downloaded in PDF format here, just don't do evil with them!

Plot 1 here
Plot 2 here
Plot 3 here
Plot 4 here
Plot 5 here
Plot 6 here


Wednesday, June 18, 2014

[SOLVED] Problems compiling Rcpp dependent R packages in Crunchbang Linux 11 (Debian 7.5 (wheezy) 64-bit)

I had a very weird issue when I tried to compile (install) certain R packages like "wordcloud", "RSNNS" or "GOSemSim" (a Bioconductor package). The installation always ended with a compiling error at the end.

At first, I tried to solve the problem by finding anyone that had faced the same issue, so I Duck-Duck-go-it, Google-it and at least for me, I did not find anyone.

At the end, and which is the core of my post, I solved the problem by removing or commenting all my customizations in the /usr/lib/R/etc/ file (for example, loading libraries automatically).

In summary, using the default worked for me.

I wrote this post to make a record of the issue, so maybe it would help someone to solve the same in the future.


Tuesday, February 25, 2014

Convert Ensembl, Unigene, Uniprot and RefSeq IDs to Symbol IDs in R using Bioconductor

Hello, I have programmed a function that converts different sources of IDs to Symbol IDs.

The input ID types allowed are (at the moment):  Ensembl, Unigene, Uniprot and RefSeq.

The code is available clicking here

NOTE: The function depends on the Bioconductor package "" available here

For example, lets show 10 Ensembl IDs:

> id[1:10]
 [1] "ENSG00000121410" "ENSG00000175899" "ENSG00000256069" "ENSG00000171428"
 [5] "ENSG00000156006" "ENSG00000196136" "ENSG00000114771" "ENSG00000127837"
 [9] "ENSG00000129673" "ENSG00000090861"

And their Symbol IDs:

> res[1:10]
 [1] "A1BG"     "A2M"      "A2MP1"    "NAT1"     "NAT2"     "SERPINA3"
 [7] "AADAC"    "AAMP"     "AANAT"    "AARS"    

This is a running example of the function to convert Unigene IDs to Symbol IDs (For all the other IDs types, just replace "unigene" to "ensembl" or "refseq" or "uniprot"):

unigene <- toTable(org.Hs.egUNIGENE)
# extract 100 random unigene entries
id  <- unigene[sample(1:length(unigene[,2]),100),2]
id.type  <- "unigene"
res <- get.symbolIDs(id,id.type)


Monday, February 10, 2014

Upgrade and update R 2.15 to R 3.0 in Debian Wheezy

Following the instructions from CRAN, you need to add the R backports in your source list.


First, open a Terminal and open the sources.list file:

$ gksudo gedit /etc/apt/sources.list

Then, add these lines at the bottom of the file (Note, I use the Revolution Analytics Dallas, TX server, but this can be easily changed taking a look here for the mirrors):

deb wheezy-cran3/
#deb-src wheezy-cran3/


There's a folder where R uses to store the packages we download, just rename it to the current version of R. For example, mine was "2.15" and then I just renamed it to "3.0" and was inside this path:


Remember that some packages also needs to install some files in folders that belongs to the root, so, I would recommend to open R in sudo mode (only if you're sure about what you're doing :P) just by executing R this way: "sudo R" and then, in the R console type :

update.packages(checkBuilt=TRUE, ask=FALSE)


The Debian backports archives on CRAN are signed with the key ID 381BA480, to add them, in a Terminal prompt type:

gpg --keyserver --recv-key 381BA480
gpg -a --export 381BA480 > jranke_cran.asc
sudo  apt-key add jranke_cran.asc


Save the file and you can either enter to Synaptic, update the packages list and then just upgrade the packages or in a terminal type:

sudo apt-get update
sudo apt-get upgrade

And that's all.

Saturday, October 19, 2013

Optimizing a multivariable function parameters using a random method, genetic algorithm and simulated annealing in R

Say that you are implementing a non-linear regression analysis, which is shortly described by wikipedia as:

"In statistics, nonlinear regression is a form of regression analysis in which observational data are modeled by a function which is a nonlinear combination of the model parameters and depends on one or more independent variables."

For the training set, we have the following:

And the function to optimize the parameters is:

Which leads us to the following equality:

In other words, we want to optimize the value of theta in order to minimize the sum of the error among y and predicted.y:

Given theta (each parameter a0,..a3 has a range from 0 to 15):

And the error function:

finally, the goal function:

In other words, the goal function searches for the value of theta that minimizes the error.


This is the scatter plot of the training set:

Here is the implementation in R, you can download the file clicking here

Here is a result plot using the genetic algorithm:


Saturday, August 17, 2013

Lieutenant Dan You Got a New Interface..

After some days of thinking, I realize that this blog deserved a little bit more attention, so I decided to change the interface and I'am happy about how it looks now.

Hope to see you again in my personal blog.

keep on programming!