Tag Archives: R

Space Invaders

I burned through all of my extra lives in a matter of minutes, and my two least-favorite words appeared on the screen: GAME OVER (Ernest Cline, Ready Player One)

Inspired by the book I read this summer and by this previous post, I decided to draw these aliens:

Invader1 Invader3
Invader4Invader2

Do not miss to check this indispensable document to choose your favorite colors:

require("ggplot2")
require("reshape")
mars1=matrix(c(0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,1,0,0,0,0,1,0,0,0,
0,0,0,0,1,0,0,1,0,0,0,0,
0,0,0,1,1,1,1,1,1,0,0,0,
0,0,1,1,0,1,1,0,1,1,0,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,1,0,1,0,0,0,0,1,0,1,0,
0,1,0,0,1,0,0,1,0,0,1,0,
0,0,0,0,0,0,0,0,0,0,0,0), nrow=12, byrow = TRUE)
mars2=matrix(c(0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,1,1,0,0,0,0,0,
0,0,0,0,1,1,1,1,0,0,0,0,
0,0,0,1,1,1,1,1,1,0,0,0,
0,0,1,1,0,1,1,0,1,1,0,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,0,0,0,1,0,0,1,0,0,0,0,
0,0,0,1,0,1,1,0,1,0,0,0,
0,0,1,0,1,0,0,1,0,1,0,0,
0,1,0,1,0,0,0,0,1,0,1,0,
0,0,0,0,0,0,0,0,0,0,0,0), nrow=12, byrow = TRUE)
mars3=matrix(c(0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,1,0,0,1,0,0,0,0,
0,0,0,1,0,0,0,0,1,0,0,0,
0,0,0,1,1,1,1,1,1,0,0,0,
0,0,1,1,0,1,1,0,1,1,0,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,1,0,1,0,0,0,0,1,0,1,0,
0,1,0,0,1,1,1,1,0,0,1,0,
0,0,0,0,1,0,0,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0), nrow=12, byrow = TRUE)
mars4=matrix(c(0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,1,1,0,0,0,0,0,
0,0,0,0,1,1,1,1,0,0,0,0,
0,0,0,1,1,1,1,1,1,0,0,0,
0,0,1,1,0,1,1,0,1,1,0,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,1,1,1,1,1,1,1,1,1,1,0,
0,1,0,0,1,0,0,1,0,0,1,0,
0,0,0,1,0,0,0,0,1,0,0,0,
0,0,0,0,1,0,0,1,0,0,0,0,
0,0,0,1,0,0,0,0,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0), nrow=12, byrow = TRUE)
opt=theme(legend.position="none",
panel.background = element_blank(),
panel.grid = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
axis.text = element_blank())
p1=ggplot(melt(mars1), aes(x=X2, y=X1))+geom_tile(aes(fill=jitter(value, amount=.1)), colour="gray65", lwd=.025)+
scale_fill_gradientn(colours = c("chartreuse", "navy"))+scale_y_reverse()+opt
p2=ggplot(melt(mars2), aes(x=X2, y=X1))+geom_tile(aes(fill=jitter(value, amount=.1)), colour="gray65", lwd=.025)+
scale_fill_gradientn(colours = c("olivedrab1", "magenta4"))+scale_y_reverse()+opt
p3=ggplot(melt(mars3), aes(x=X2, y=X1))+geom_tile(aes(fill=jitter(value, amount=.1)), colour="gray65", lwd=.025)+
scale_fill_gradientn(colours = c("violetred4", "yellow"))+scale_y_reverse()+opt
p4=ggplot(melt(mars4), aes(x=X2, y=X1))+geom_tile(aes(fill=jitter(value, amount=.1)), colour="gray65", lwd=.025)+
scale_fill_gradientn(colours = c("tomato4", "lawngreen"))+scale_y_reverse()+opt

Princess Jasmine’s Trick

I’m history! No, I’m mythology! Nah, I don’t care what I am; I’m free hee! (Genie, when he is released from the magical oil lamp by Aladdin)

Princess JasmineA long time ago, in a kingdom far away, lived a beautiful princess named Jasmine. There also lived a very rich and evil wizard named Jafar, who was in love with the princess. In order to married with Jasmine, Jafar bought  her father’s will with treasures, but the princess was harder to convince. One day Jafar told the princess: Request me whatever you want and if  I am able to bring it to you, you will become my wife. The princess, tired of the insistence of Jafar, answered: I only want a gold chain, but I want you to give it to me as follows: the first day I should have just one link of the chain. The second day I should have two links. The third day, three … and so on. When you give me all the links of the chain I will marry you. Jafar, intrigued, asked: But how many links should have the chain?  And Jasmine replied: I want you to give me the longest chain that allows you to pay me breaking only 30 links. Jafar began to laugh out loud as he walked away and said to the princess: Tomorrow I’ll bring you such chain!. But as he went to his palace, his happiness turned into anger: he realized that there was not enough gold in the world to build the chain that asked Jasmine.

This is my own version of one of my favorite anti-common-sense mathematical curiosities. To explain it, let me start with an example. Imagine a simple chain with 7 links. If you open the 3rd link, the you split the chain into 3 pieces: a single link (the one you opened), a piece of 2 links and another one of 4 links. You could pay to Jasmine during seven days combining these 3 pieces:

  • Day 1: Give her the single link
  • Day 2: Give her the 2-links piece and take the single link, leaving her with 2 links
  • Day 3: Give her the single link again, leaving her with 3 links
  • Day 4: Give her the 4-links piece and take all pieces she has, leaving her with 4 links
  • Day 5: Give her the single link again, leaving her with 5 links
  • Day 6: Give her the 2-links piece and take 2-links piece, leaving her with 6 links
  • Day 7: Give her the single link piece, leaving her with all links

Is easy to see that having a chain with 63 links, you could pay Jasmine breaking only 3 links (positions 5th, 14th and 31st). It easy to prove that the length of the biggest chain you can manage breaking only n links is (2n+1-1)*(n+1)+n

Next plot represents the minimum number of breaks to pay Jasmine daily for a given chain’s length. I call it the Jasmine’s Staircase:

Generated With R #rstats

Some curiosities around chains:

  • Jasmine asked Jafar a chain of 66.571.993.087 links
  • Supposing one link weights 4 grams, the chain of Jasmine would weight around 266 tons. It is supposed to be around 171 tons of gold in the world
  • If you spend 1 second to climb the first step of the staircase, you will spend 302 years to climb the step number 100

Jafar was right. Jasmine was clever:

library(sqldf)
library(ggplot2)
library(extrafont)
max.breaks=5
CalculateLength = function(n) {n+sum(sapply(0:n, function(x) 2^x*(n+1)))}
results=data.frame(breaks=1:max.breaks, length=sapply(1:max.breaks, CalculateLength))
links=data.frame(links=2:CalculateLength(max.breaks))
results=sqldf("SELECT links.links, min(results.breaks) as minbreaks FROM links, results WHERE links.links <= results.length GROUP BY 1")
opts=theme(
panel.background = element_rect(fill="mistyrose"),
panel.border = element_rect(colour="black", fill=NA),
axis.line = element_line(size = 0.5, colour = "black"),
axis.ticks = element_line(colour="black"),
panel.grid = element_line(colour="white", linetype = 2),
axis.text.y = element_text(colour="black"),
axis.text.x = element_text(colour="black"),
text = element_text(size=20, family="Humor Sans"),
plot.title = element_text(size = 40)
)
ggplot(results, aes(links,minbreaks))+
geom_area(fill="violet", alpha=.4)+
geom_step(color="violetred", lwd=1.5)+
labs(x="Chain's Length", y="Minimum Number of Breaks", title="Princess Jasmine's Staircase")+
scale_x_continuous(expand = c(0, 0), breaks = sapply(1:max.breaks, CalculateLength))+
opts

Looking For Life

 

Machines take me by surprise with great frequency (Alan Turing)

Imagine a 8×8 grid in which cells can be alive (black colored) or empty (light gray colored): Generated with R #RstatsAs with the One-dimensional Cellular Automata, the next state of a cell is a function of the states of the cell’s nearest neighbors (the only neighbors that we consider are the eight cells that form the immediate perimeter of a cell). To evolve this community of cells over time, rules are quite simple:

  • If a live cell has less than two neighbors, then it dies (loneliness)
  • If a live cell has more than three neighbors, then it dies (overcrowding)
  • If an dead cell has three live neighbors, then it comes to life (reproduction)
  • Otherwise, a cell stays as is (stasis)

These are the rules of the famous Conway’s Game Of Life, a cellular automaton invented in the late 1960s by John Conway trying to refine the description of a cellular automaton to the simplest one that could support universal computation. In 1970 Martin Gardner described Conway’s work in his Scientific American column. Gardner’s article inspired many people around the world to experiment with Conway’s, which eventually led to the final pieces of how the Game Of Life could support universal computation in what was surely a global collaborative effort.

In this experiment I will try to find interesting objects that can be found in the Game Of Life: static (remain the same over the time), periodic (change but repeating their initial configuration in some iterations) or moving objects (travel through the grid before disappearing). Why are interesting? Because these are the kind of objects required to perform computations.

The plan is simple: I will generate some random grids and will evolve them over time a significant number of times. After doing this, I will check those grids having still some alive cells inside. Will I find there what I am looking for?

I generated 81 grids in which live cells are randomly located using binomial random variables with probabilities equal to i/80 with i from 0 (all cells empty) to 80 (all cells alive). This is a quick way to try a set of populations with a wide range of live cells. I measure % of alive cells after each iteration. I will analyze those grids which still have live cells after all iterations. This is what happens after 150 iterations:

Generated with R #Rstats

I find some interesting objects. Since I keep them in my script, I can list them with ls(pattern= "Conway", all.names = TRUE). Two of them are specially interesting because are not static. Are those which produce non-constant lines in the previous plot.

First one is a periodic object which reproduces itself after every 3 iterations:

Generated with R #Rstats

Second one is a bit more complex. After 8 iterations appears rotated:

Generated with R #Rstats

What kind of calculations can be done with these objects? I don’t now yet but let’s give time to time. Do you want to look for Life?

library(ggplot2)
library(scales)
library(gridExtra)
SumNeighbors = function (m) #Summarizes number of alive neighbors for each cell
{
  rbind(m[-1,],0)+rbind(0,m[-nrow(m),])+cbind(m[,-1],0)+cbind(0,m[,-ncol(m)])+
    cbind(rbind(m[-1,],0)[,-1],0)+cbind(0, rbind(0,m[-nrow(m),])[,-nrow(m)])+
    cbind(0,rbind(m[-1,],0)[,-nrow(m)])+cbind(rbind(0,m[-nrow(m),])[,-1],0)
}
NextNeighborhood = function (m) #Calculates next state for each cell according to Conway's rules
{
  (1-(SumNeighbors(m)<2 | SumNeighbors(m)>3)*1)-(SumNeighbors(m)==2 & m==0)*1
}
splits=80 #Number of different populations to simulate. Each population s initialized randomly
          #according a binomial with probability i/splits with i from 0 to splits
iter=150
results=data.frame()
rm(list = ls(pattern = "conway")) #Remove previous solutions (don't worry!)
for (i in 0:splits)
{
  z=matrix(rbinom(size=1, n=8^2, prob=i/splits), nrow=8); z0=z
  results=rbind(results, c(i/splits, 0, sum(z)/(nrow(z)*ncol(z))))
  for(j in 1:iter)
  {z=NextNeighborhood(z); results=rbind(results, c(i/splits, j, sum(z)/(nrow(z)*ncol(z))))}
  #Save interesting solutions
  if (sum(z)/(nrow(z)*ncol(z))>0) assign(paste("Conway",format(i/splits, nsmall = 4), sep=""), z)
}
colnames(results)=c("start", "iter", "sparsity")
#Plot reults of simulation
opt1=theme(panel.background = element_rect(fill="gray85"),
          panel.grid.minor = element_blank(),
          panel.grid.major.x = element_blank(),
          panel.grid.major.y = element_line(color="white", size=.5, linetype=2),
          plot.title = element_text(size = 45, color="black"),
          axis.title = element_text(size = 24, color="black"),
          axis.text = element_text(size=20, color="black"),
          axis.ticks = element_blank(),
          axis.line = element_line(colour = "black", size=1))
ggplot(results, aes(iter, sparsity, group = start))+
  geom_path(size=.8, alpha = 0.5, colour="black")+
  scale_x_continuous("Iteration", expand = c(0, 0), limits=c(0, iter), breaks = seq(0,iter,10))+
  scale_y_continuous("Alive Cells", labels = percent, expand = c(0, 0), limits=c(0, 1), breaks = seq(0, 1,.1))+
  labs(title = "Conway's Game Of Life Simulation")+opt1
#List of all interesting solutions
ls(pattern= "Conway", all.names = TRUE)
#Example to plot the evolution of an interesting solution (in this case "Conway0.5500")
require(reshape)
opt=theme(legend.position="none",
            panel.background = element_blank(),
            panel.grid = element_blank(),
            axis.ticks=element_blank(),
            axis.title=element_blank(),
            axis.text =element_blank())
p1=ggplot(melt(Conway0.5500), aes(x=X1, y=X2))+geom_tile(aes(fill=value), colour="white", lwd=2)+
  scale_fill_gradientn(colours = c("gray85", "black"))+opt
p2=ggplot(melt(NextNeighborhood(Conway0.5500)), aes(x=X1, y=X2))+geom_tile(aes(fill=value), colour="white", lwd=2)+
  scale_fill_gradientn(colours = c("gray85", "black"))+opt
p3=ggplot(melt(NextNeighborhood(NextNeighborhood(Conway0.5500))), aes(x=X1, y=X2))+geom_tile(aes(fill=value), colour="white", lwd=2)+
  scale_fill_gradientn(colours = c("gray85", "black"))+opt
p4=ggplot(melt(NextNeighborhood(NextNeighborhood(NextNeighborhood(Conway0.5500)))), aes(x=X1, y=X2))+geom_tile(aes(fill=value), colour="white", lwd=2)+
  scale_fill_gradientn(colours = c("gray85", "black"))+opt
#Arrange four plots in a 2x2 grid
grid.arrange(p1, p2, p3, p4, ncol=2)

The Andrica’s Conjecture

Things should be as simple as possible, but not simpler (Albert Einstein)

Following with conjectures about primes, it is time for Andrica’s conjecture. The great mathematician Leonhard Euler (1707-1783) pointed: “Mathematicians have tried with no success to find some kind of order in the sequence of prime numbers and today we have reasons to believe that this is a mystery that human mind will never understand”.

In 1985, the Romanian mathematician Dorin Andrica published his conjecture, still unproved, which makes reference to gap between consecutive prime numbers. In concrete, his conjecture establishes that difference between square roots between two consecutive prime numbers is always less than 1. The highest difference encountered until now is 0.67087, located between p4=7 and p5=11.

Following you can find the plot of these differences for first 400 prime numbers:

andricaIt is very interesting how dots form hyperbolic patterns. Does not seem similar in some sense to the Ulam spiral? Primes: how challenging you are!

Two more comments:

  • It is better to find primes using matlab package than doing with schoolmath one. Reason is simple: for schoolmath package, 133 is prime!
  • Why did Andrica formulated his conjecture as √pn+1-√pn < 1 instead of √pn+1-√pn < 3/4? In terms of statistical error, the second formulation is more accurate. Maybe the charisma of number 1 is hard to avoid.

This is the code. I learned how to insert mathematical expressions inside a ggplot chart:

library(matlab)
library(ggplot2)
ubound=2800
primes=primes(ubound)
andrica=data.frame(X=seq(1:(length(primes)-1)), Y=diff(sqrt(primes)))
opt=theme(panel.background = element_rect(fill="gray92"),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.major.y = element_line(color="white", size=1.5),
plot.title = element_text(size = 45),
axis.title = element_text(size = 28, color="gray35"),
axis.text = element_text(size=16),
axis.ticks = element_blank(),
axis.line = element_line(colour = "white"))
ggplot(andrica, aes(X, Y, colour=Y))+geom_point(size=5, alpha=.75)+
scale_colour_continuous(guide = FALSE)+
scale_x_continuous("n", limits=c(0, length(primes)-1), breaks = seq(0,length(primes)-1,50))+
scale_y_continuous(expression(A[n]==sqrt(p[n+1])-sqrt(p[n])), limits=c(0, .75), breaks = seq(0,.75,.05))+
labs(title = "The Andrica's Conjecture")+
opt

Summer Summary

The universe is full of magical things patiently waiting for our wits to grow sharper (Eden Phillpots)

I launched this blog 7 months ago and published 30 posts during this time. These are some of my figures until now:

My favourite post? I don’t really know, but I am very proud of this one and this one. I have received more positive critics than negative ones and the future sounds good: I have lots of experiments in my head to try with R.

Thanks a lot.

map

Batman’s Choice

A hero can be anyone, even a man doing something as simple and reassuring as putting a coat on a young boy’s shoulders to let him know the world hadn’t ended (Batman in The Dark Knight Rises)

Joker has captured Batman and keeps him into a dark and cold dungeon of Gotham City. Showing his sadistic character, Joker proposes the following game to Batman:

This is a six shooter revolver with two bullets in the cylinder. Bullets are inside two consecutive chambers. I will spin the cylinder and will fire the gun aiming to my head. If I survive you will have to do the same but you decide if you want to spin the cylinder before firing or not. If you still keep you head over your shoulders after firing, you will be free.

Joker fires and nothing happens. He survives and passes the revolver to Batman. What should Batman do? Spinning or not? What would you do?

From my point of view, answer is quite anti-intutive because the best option is not spinning the cylinder again. Spinning case is clear: probability of losing the head is 2/6=33% but what about not spinning? Doing the next shoot directly eliminates two possibilities: the previous shot of Joker and the second bullet according to direction of cylinder rotation (remember two bullets are consecutive and Joker is unfortunately  still alive). It means there is only 1 chance to dead between 4, so probability of losing the head in this scenario is 1/4=25% which is significantly lower than the first one.

Here you can find the resulting graph of simulating the game up to 500 times:

batmans_choice2

Will it be the end of Batman? Not sure.

This is the code of this experiment:

library(ggplot2)
library(extrafont)
niter=500
results=data.frame()
for (i in 1:niter)
{
bullet1=sample(1:6,1)
Joker=sample((1:6)[-c(bullet1, bullet1%%6+1)],1)
#Option 1: Shooting
Batman1=Joker%%6+1
dead1=(Batman1 %in% c(bullet1, bullet1%%6+1))*1
#Option 2: Rolling and Shooting
Batman2=sample(1:6,1)
dead2=(Batman2 %in% c(bullet1, bullet1%%6+1))*1
results=rbind(results, c(i, dead1, dead2))
}
colnames(results)=c("iter", "dead1", "dead2")
results$csum1=cumsum(results$dead1)/as.numeric(rownames(results))
results$csum2=cumsum(results$dead2)/as.numeric(rownames(results))
theme_xkcd=theme(
panel.background = element_rect(fill="darkolivegreen1"),
panel.border = element_rect(colour="black", fill=NA),
axis.line = element_line(size = 0.5, colour = "black"),
axis.ticks = element_line(colour="black"),
panel.grid = element_line(colour="white", linetype = 2),
axis.text.y = element_text(colour="black"),
axis.text.x = element_text(colour="black"),
text = element_text(size=18, family="Humor Sans"),
plot.title = element_text(size = 50)
)
p=ggplot(data=results, aes(x=iter, y=csum1))+
geom_abline(intercept = 1/4, slope = 0, size = 0.4, linetype=2, colour = "black", alpha=0.8)+
geom_abline(intercept = 1/3, slope = 0, size = 0.4, linetype=2, colour = "black", alpha=0.8)+
geom_line(aes(y=csum2), colour="green4", size=1.5, fill=NA)+
geom_line(colour="green4", size=1.5, fill=NA)+
coord_cartesian(ylim=c(.0, 1), xlim=c(1, niter))+
scale_y_continuous(breaks = c(0,round(1/4, digits = 2),round(1/3, digits = 2),1))+
geom_text(data=results[niter*.75, ], family="Humor Sans", colour="green4", y=0.38, label="Rotating Cylinder and Shooting ...", size=4, adjust=1)+
geom_text(data=results[niter*.75, ], family="Humor Sans", colour="green4", y=0.20, label="Shooting without Rotating Cylinder ...", size=4, adjust=1)+
labs(x="Number Of Trials", y="Prob. of Losing The Head", title="Batman's Choice")+
theme_xkcd
ggsave("batmans_choice.jpg", plot=p, width=8, height=5)

The Zebra Of Riemann

Mathematics is the art of giving the same name to different things (Henri Poincare)

Many surveys among experts point that demonstration of the Riemann Hypothesis is the most important pending mathematical issue in this world. This hypothesis is related to Riemann zeta function, which is supossed to be zero only for those complex whose real part is equal to 1/2 (this is the conjecture of Riemann itself). Confirming the conjecture would imply deep consequences in prime numbers teory and also in our knowledge of their properties. Next plot represents the argument of zeta function over complex with real part between -75 and 5 (x axis) and imaginary part between -40 and 40 (y axis). An explanation of this kind of graph can be found here. Does not it remind you of something?
zebra
Solving the hypothesis of Riemann is one of the seven Clay Mathematics Institute Millennium Prize Problems. Maybe zebras keep the secret to solve it in their skins.

This is the code to plot the graph:

require(pracma)
z=outer(seq(-75, 5, by =.1),1i*seq(-40, 40, by =.1),'+')
z=apply(z, c(1,2), function(x) Arg(zeta(x)))
opt=theme(legend.position="none",
panel.background = element_blank(),
panel.grid = element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text =element_blank())
z=data.frame(expand.grid(x=seq(ncol(z)), y=seq(nrow(z))), z=as.vector(z))
require(ggplot2)
ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=c("black","white", "gray80", "gray40")) + opt

Four Simple Turtle Graphs To Play With Kids

Technology is just a tool: in terms of getting the kids working together and motivating them, the teacher is the most important (Bill Gates)

Some days ago I read this article in R-bloggers and I discovered the TurtleGraphics package. I knew about turtle graphics long time ago and I was thinking of writing a post about them, doing some draw on my own using this technique. In fact, some of my previous drawings could have been done using this package.

Turtle graphics are a good tool to teach kids some very important mathematical and computational concepts such as length, coordinates, location, angles, iterations … Small changes in some of this parameters (especially in angles) can produce very different drawings. Those who want to know more about turtle graphics can read the Wikipedia article.

Following you can find four simple graphics preceded by the code to generate them. Please, send me your creations if you want:

library(TurtleGraphics)
turtle_init()
turtle_col("gray25")
for (i in 1:150) {
  turtle_forward(dist=1+0.5*i)
  turtle_right(angle=89.5)}
turtle_hide()

turtle1

library(TurtleGraphics)
turtle_init()
turtle_col("gray25")
turtle_right(angle=234)
for (i in 1:100) {
  turtle_forward(dist=0.9*i)
  turtle_right(angle=144.3)}
turtle_hide()

turtle2

library(TurtleGraphics)
turtle_init()
turtle_col("gray25")
turtle_setpos(48,36)
d=50
for (i in 1:300) {
  turtle_forward(dist=d)
  if (i%%4==0) {
    turtle_right(angle=75)
    d=d*.95}
else turtle_right(angle=90)}
turtle_hide()

turtle3

library(TurtleGraphics)
turtle_init()
turtle_col("gray25")
turtle_setpos(50,35)
turtle_right(angle=30)
d=25
turtle_setpos(50-d/2,50-d/2*tan(pi/6))
for (i in 1:100) {
  turtle_forward(dist=d)
  d=d+.5
  turtle_right(angle=120+1)}
turtle_hide()

turtle4

How To Approximate Pi With A Short Pencil And A Big Paper

Experiment, be curious: though interfering friends may frown, get furious at each attempt to hold you down (Tony Bennett, Experiment)

Instructions:

  1. Take a pencil and measure it
  2. Take a piece of paper and draw parallel lines on it (you can use the pencil, of course); separation between lines should double the length of the pencil
  3. Toss the pencil over the paper 100 times (or more)
  4. Make note of how many times do the pencil cross some of the lines
  5. Calculate ratio between tosses and crosses: this is your approximation of Pi

Some time ago, I published a post about one of the most amazing places where PI was discovered. This is another example of the ubiquity of this mathematical constant. This experiment is based on Buffon’s needle problem, another amazing experiment of 18th century. Next plot represents ratio of tosses to crosses depending on the length of pencil. When the pencil is half the length of the separation between lines, the previous ratio is approximately PI:
Buffon
If you get very bored some afternoon you can replicate this experiment with your children. Use a short pencil. If not, you will need an extremely big piece of paper. Meanwhile here you have the code:

trials=100000
results=sapply(seq(.1, 2, by = .05), function(x)
{
  r=x #Length of pencil in relation to separation between lines
  Needles=t(sapply(1:trials, function(y) c(100*runif(1),2*pi*runif(1))))
  Needles=cbind(Needles,Needles[,1]+r*cos(Needles[,2]))
  Needles=data.frame(x1=Needles[,1], x2=Needles[,3], Cross=abs(trunc(Needles[,1])-trunc(Needles[,3])))
  c(r, trials/(trials-nrow(Needles[Needles$Cross==0,])))
})
results=data.frame(t(results))
colnames(results)=c(c("ratio", "inv.perc.crosses"))
require(ggplot2)
require(extrafont)
require(grid)
opts=theme(
  panel.background = element_rect(fill="darkolivegreen1"),
  panel.border = element_rect(colour="black", fill=NA),
  axis.line = element_line(size = 0.5, colour = "black"),
  axis.ticks = element_line(colour="black"),
  panel.grid.major = element_line(colour="white", linetype = 1),
  panel.grid.minor = element_blank(),
  axis.text.y = element_text(colour="black", size=20),
  axis.text.x = element_text(colour="black", size=20),
  text = element_text(size=25, family="xkcd"),
  legend.key = element_blank(),
  legend.position = c(.2,.75),
  legend.background = element_blank(),
  plot.title = element_text(size = 50))
c=ggplot(results, aes(ratio, inv.perc.crosses))
c +geom_abline(intercept = pi, slope = 0, size = 0.4, linetype=2, colour = "black", alpha=0.8)+
  geom_line(color="green4", size=1.5)+
  geom_point(color="gray92", size=8, pch=16)+
  geom_point(color="green4", size=6, pch=16)+
  geom_text(aes(0.55, 5), hjust=0, family="xkcd", label="PI is around here!", size=10)+
  ggtitle("Hot to approximate PI with \na short pencil and a big paper")+
  xlab("Length of pencil divided by separation between lines") +
  ylab("Number of tosses divided by number of crosses")+
  geom_segment(aes(x = .625, y = 4.6, xend = .52, yend = 3.45), size=1.5, arrow = arrow(length = unit(0.5, "cm")))+
  scale_y_continuous(breaks=c(pi, 4, 8, 12, 16), labels=c("3.141593","4","8","12","16"))+
  scale_x_continuous(breaks=seq(.1, 2, by = .1))+opts

The Goldbach’s Comet

Every even integer greater than 2 can be expressed as the sum of two primes (Christian Goldbach, 1742)

The point cloud known as Goldbach’s Comet represents the amount of different ways (y axis) an even number (x axis) can be writen as sum of two prime numbers. In this plot, x axis is between 2 and 50.000 (stars are just ornaments):

Goldbach

Mathematicians are still waiting for a proof of this conjecture. This is the code for drawing this plot:

require(schoolmath)
require(utils)
library(plyr)
n=50000
data(primlist)
primes=as.data.frame(t(combn(primlist[primlist>2 & primlist<n-2], 2)))
primes$V3=primes$V1+primes$V2
primes2=count(primes, "V3")
primes2=primes2[primes2$V3<=n,]
stars=cbind(runif(50, min=-n*0.05, max=n), runif(200, min=-n*0.001, max=max(primes2$freq)))
plot.new()
par(mai = rep(0, 4), bg = "gray12")
plot(NA,type="n", xlim=c(-n*0.05,n), xaxs="i", ylim=c(-n*0.001,max(primes2$freq)))
points(stars, col = "blue4", cex=.7, pch=16)
points(stars, col = "blue", cex=.3, pch=16)
points(stars, col = "gray75", cex=.1, pch=16)
apply(primes2, 1, function(x) points(x=x[1],y=x[2], col = if (runif(1)&gt;x[1]/n) {"white"} else {sample(colours(),1)}, cex=.1, pch=16))