Drawing an Icosahedron with R

Figure 1. A blue icosahedron.

Figure 1. A blue icosahedron.

Allan Roberts

To construct an icosahedron, I began with the following premises: all edges have the same length, all faces are triangles, and five faces meet at every vertex. I then calculated the coordinates for twelve vertices equidistant from the origin. The script given below will produce an image like Figure 1. For brevity, I have simply given the vertex coordinates and shades of the faces rather than providing the underlying calculations. Do the x, y and z coordinates provided indeed describe an approximate icosahedron? If you want to, you can verify this by computing the distances between points.

R Script

#Blue planet icosahedron.
#Script written by Allan Roberts, April 2013.

X = c(1, 0.309, -0.809, -0.809,  0.309,  0.809, -0.309, -1.000, -0.309,  0.809,  0,  0);
Y = c(0,  0.951,  0.588, -0.588, -0.951,  0.588,  0.951,  0, -0.951, -0.588,  0,  0);
Z = c(0.5,  0.5,  0.5,  0.5,  0.5, -0.5, -0.5, -0.5, -0.5, -0.5,  1.118, -1.118);

Vertices = data.frame(X,Y,Z);

par(bg=gray(0));
plot.new();
par(usr=c(-3,3,-3,3));
shade = c(0.67, 0.89, 0.66, 0.30, 0.30, 0.30, 0.66, 0.30, 0, 0);

F1 = Vertices[c(1,2,12), ];
F2 = Vertices[c(2,3,12), ];
F3 = Vertices[c(3,4,12), ];
F4 = Vertices[c(4,5,12), ];
F5 = Vertices[c(5,1,12), ];
F6 = Vertices[c(1,2,6), ];
F7 = Vertices[c(2,3,7),];
F8 = Vertices[c(3,4,8),];
F9 = Vertices[c(4,5,9),];
F10 = Vertices[c(5,1,10),];

Faces = list(F1,F2,F3,F4,F5,F6,F7,F8,F9,F10);
for (i in 1:10) names(Faces[[i]]) = c(quote(X),quote(Y),quote(Z) );
for (i in 1:10) polygon(Faces[[i]]$X,Faces[[i]]$Y, col=rgb(0,0,0.2 + 0.8*shade[i]) );

Advertisements

R Bootcamp 2013 at BMSC

A beautiful Sunday in Bamfield (photo credit: Gillian Walker).

A beautiful Sunday in Bamfield (photo credit: Gillian Walker).

By Christina Suzanne

This year’s R bootcamp ran from March 3-7, 2013, immediately following the
Pacific Ecology and Evolution Conference (PEEC 2013). The workshop was
taught by Dr. Jean Richardson and Dr. Brad Anholt. The participants
ranged from R novice to R veterans; however, everyone learned something
useful and new! The five intensive days covered such things as data
input/output, managing data, scatterplots, error bars, maps,
multi-panel figures, t-tests, ANOVA, multivariate analysis, mixed effect
models, and much, much more! It also introduced various R packages, including ggplot2, vegan and lme4. Everyone left the workshop feeling much more
confident in R, and willing to tackle their own statistical problems using
R. Attending a future R bootcamp at BMSC is highly recommended to anyone who wants to start learning the R language, or increase their understanding of it!

Inside the black box, part 1: the one sample t-test

Figure 1. A figure produced using grid library view ports.

Figure 1. A figure produced using grid library view ports.

In a typical introductory statistics course, you’re given some assumptions, and provided with formulas for churning out p-values and confidence intervals. But statistics might be left feeling rather like a black box: we’re turning input into output, but it might seem difficult  to keep track of how all those assumptions and formulas work (or don’t).

One of the simplest stats tests that you’ll see in a first statistics course is the one sample t-test for testing the null hypothesis that the population mean for a random variable with a normal distribution is equal to zero. To implement this test in R all you need is the t.test( ) function. But what are the inner workings? The R script given below generates a pseudorandom sample and then puts it through the standard R function for a t.test, as well as a user defined version that should give you just about the same p-value. But wait… inside the t-test “box” other mysteries are to be found. For example: Where did the t-distribution come from in the first place? How should you interpret that p-value? And what’s up with that null hypothesis anyways? Sometimes asking questions leads to more questions… .

You can run the script provided below by cutting and pasting it onto the R command line, and pressing <enter>.

R Script

data = rnorm(20,0,1);

T.TEST = function(data){
#a one sample t test against the null hypothesis that the population mean is zero.

n = length(data);
df = n-1;
t = sqrt(n)*mean(data)/sd(data);
p.value = 2 * (pt(-abs(t),df));
output=list(t,p.value);
names(output) = c(quote(t),quote(p.value));
return(output);
}

t.test(data);
T.TEST(data);

References

R Core Team (2012). R: A language and environment for statistical computing. R
Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org/

 

A Seastar Rearrangement

Figure 1. If you look closely you may recognize that this is a scrambled image of the Madreporite logo.

Figure 1. If you look closely you may recognize the scrambled pieces of the Madreporite logo.

By Allan Roberts

One of the things that makes the stats and graphics application R so flexible is that there are lots of extension packages, or libraries, available. Some of these come with R by default, and others are available from the R Project website. The image in Figure 1 was made by using a function in the “png” library to read a .png image of the Madreporite logo into R. I used the “grid” library to create a view port for each tile. Then I used a modified version of the script posted in the Sliding Tile Game with R post to scramble the tiles.

References

Urbanek, S., 2011. png: Read and write PNG images. R package
version 0.1-4. http://CRAN.R-project.org/package=png

R Core Team (2012). R: A language and environment for statistical computing. R
Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org/

Plotting Data on an Angle

Figure 1. A line graph at a  30 degree angle.

Figure 1. A line graph on angle.

You can try to make your data more dynamic looking by taking the y-axis off the perpendicular. Of course, that might be bad judgement. This blog post does not advocate rotating the frame of reference for your data plots. It merely shows you how.

The example provided below uses the “grid” package. If you don’t have this installed already, you can do so by making sure you are connected to the internet, and then entering “install.packages(“grid”)” on the R command line. I figured out how to use the grid library commands in this example by referring to Murrell (2006).

R Script
#Example written by Allan Roberts, April 2013

library(grid);
pushViewport(viewport());
pushViewport(viewport(height=0.5,width=0.5,angle=30));
time = seq(0,1,0.05);
x = rnorm(21,0.5,0.1);
grid.lines(time,x);
grid.xaxis(gp=gpar(col=rgb(1,0,0)));
grid.yaxis(gp=gpar(col=rgb(0,0,1)));

Reference

Murrell, Paul. 2011. R Graphics, 2nd Edition. CRC Press.

Buffon’s Needle

By Allan Roberts

Figure 1. A simulation of Buffon's Needle Experiment.

Figure 1. A simulation of Buffon’s Needle Experiment. An estimate for Pi has been computed from the proportion of needles crossing the horizontal lines (n = 100,000; not all lie within the margins of the figure).

Today is π Day. Which is sort of actually a real holiday (Lamb 2012). (But not, presumably, a rational holiday.) And what better way to celebrate π Day than with a probabilistic estimation of π? The R script provided below will simulate a version of Buffon’s Needle experiment, and use the proportion of needles crossing the horizontal lines to estimate π. (For details on Buffon’s Needle see the Wikipedia article cited below.)

π, what is it good for? Well, it gets used in the formulas for the area of a circle, and the volume of a sphere. It’s also a constant in the formula for the probability density function for the normal distribution. Thus every time you make an assumption of normality, you are making a tacit reference to π.

R Script

BuffonsNeedle <- function(n=100, sd=2){
#Implemented in R by Allan Roberts, 2013.

plot.new();
par( usr = c(-2,12,-2,12), mar=c(4,4,4,4) );
arrows(rep(-4,16),-2:12, rep(14,16),-2:12, length=0);

X <- rnorm(n,5,sd);
Y <- rnorm(n,5,sd);
Angle <- runif(n,0,2*pi); #Note: the script itself uses pi. Homework: write a pi-less version.
X2 <- X + cos(Angle);
Y2 <- Y + sin(Angle);
CrossesLine <- ( floor(Y) != floor(Y2) );

arrows(X[!CrossesLine],Y[!CrossesLine],X2[!CrossesLine],Y2[!CrossesLine],length=0,col=rgb(0,1,0),lwd=1);

arrows(X[CrossesLine],Y[CrossesLine],X2[CrossesLine],Y2[CrossesLine],length=0,col=rgb(0,0.5,0,0.5),lwd=1);

p <- sum(CrossesLine)/n;
return(c(p,2/p));
}

BuffonsNeedle( );

References

Buffon’s needle. Wikipedia. http://en.wikipedia.org/wiki/Buffon%27s_needle

Lamb, Evelyn. July 21, 2012. Pi Approximation Day Celebrated July 22, How Much Pi Do You Need? Huffington Post.

Pi Day. Wikipedia.

Logistic Growth Model with R

Figure 1. Example.

Figure 1. This graph was produced by the R script provided below.

Models like the discrete logistic growth model are famous for producing complex behaviour from simple equations. You can cut and paste the R script provided below onto the R command line, to produce a graph like the one given Figure 1. Try varying the parameter r to see what happens; for example, to set the parameter r to 1.5, enter the following:

> LogisticGrowthExample( r = 1.5 );

References

Logistic Equation. Wolfram Mathworld: mathworld.wolfram.com/LogisticEquation.html
Logistic Map. Wikipedia: en.wikipedia.org/wiki/Logistic_map

R Script

#NOTE: You will need to replace the right and left quotation marks with
#straight quotes before running this script in R.

LogisticGrowthExample <- function(r = 3.7, N0 = 0.1, steps = 100){
#Implemention by Allan Roberts, March 2013.

N <- numeric(steps);
N[1] <- N0; #the initial population

Next <- function(N) {N <-  N + (r-1)*N*(1-N)};
for (i in 2:(steps)) N[i] <- Next(N[i-1]);

plot(N~seq(0,(steps-1)),type=”l”,xlab=quote(time), ylim=c(0,1.5), las = 1 );
title(main=”Discrete Logistic Growth Model”, font.main=1, cex.main=2);
text(80,0.2,paste(“r =”,as.character(r) ), adj=c(0,0.5));
text(80,0.1, expression(“N”[0]), adj=c(0,0.5 ) );
text(84.14,0.1, paste(” =”,as.character(N[1])), adj=c(0,0.5) );
}

LogisticGrowthExample( );