Just for Fun

8月 122013
 

Editor's Note: My 8th grade son, David, created a poster that he submitted to the 2013 ASA Poster Competition. The competition encourages students to display "two or more related graphics that summarize a set of data, look at the data from different points of view, and answer specific questions about the data." In this guest blog, David explains his experiment and analysis. All prize-winning posters are featured in the August 2013 issue of Amstat News.
-- Rick

Hi! My name is David Wicklin. This blog post describes a research project and statistical graphs that I created for the 2013 ASA Poster Competition.

My poster began one day when my sister saw a small box on a store shelf that contained two spiked plastic balls. The balls were about the size of tennis balls. The box said that you should put these "dryer balls" in your dryer with a load of wet laundry. The box claimed that the balls would "reduce drying time by up to 25%." I was skeptical. Could putting plastic balls in the dryer really reduce drying time? I was looking for a project for my school science fair, so I decided to conduct an experiment to find out.

My family's dryer has a "dampness sensor" that beeps when a load of laundry is dry. I decided to use this feature and a stop watch to determine the drying time for a load of laundry. For half the loads, I'd toss in two dryer balls. For others I would not.

Experimental Method

The first issue was how to design the experiment. One option was a controlled experiment in which I would select clothes of different materials (jeans, sweatshirts, tee shirts, and so forth) and different quantities (small, medium, and large loads). I would wash and dry each load twice: once with dryer balls and once without. That would make comparing drying times easy because the only varying factor would be the dryer balls; the clothes would be exactly the same.

However, from a practical point of view, a controlled experiment was not ideal. It would require wasting time and energy to wash and dry items twice. Plus, my family had real laundry that needed washing and drying! My dad suggested that we wash the regular family clothes and compare the mean drying times. If I washed enough laundry, I should be able to account for the random differences between the loads, such as the materials being dried. I would weigh the clothes to keep track of the size of the load.

I got my mom to agree to let me dry our family's laundry over the next few months. As you might suspect, she was very much in favor of this activity! I dried a total of 40 loads of laundry. For each load, I weighed the wet clothes by using a spring scale. For 20 randomly selected loads, I added two dryer balls to the dryer. The other 20 loads had no balls, just the clothes. Each load was dried using the same dryer setting. All types of clothing were used to randomize the samples. For each load, I recorded the weight and the drying time in a notebook.

Analyzing and graphing the data

I created two different graphs. The first is a side-by-side box and whisker diagram. This shows that the average drying times are not separated by much. For loads without dryer balls, the average drying time was 28.4 minutes. The median drying time was 27.5 minutes, and the first and third quartiles were 22.25 and 34 minutes. For loads with dryer balls added, the average drying time was 28.3 minutes. The median drying time was 28.75 minutes, and the first and third quartiles were 23 and 33 minutes.

The box and whisker diagram shows that it does not matter if you add dryer balls or not. If dryer balls really "reduce drying time by up to 25%," I would expect the red box plot to be shifted down by 5–8 minutes as compared to the blue box. It isn't. The two boxes are essentially the same.

But maybe dryer balls only "work" for big loads (or small loads). To find out, I created a scatter plot that displays the drying time and the weight of the wet clothes. I plotted blue squares for drying times that did not use dryer balls. I plotted red circles for times that included dryer balls.

The scatter plot shows that dryer balls do not affect drying time whether it is a heavy load or a light one. For light loads, it takes about 22 minutes to dry the laundry. For heavy loads it takes about 35 minutes. The red and blue markers are mixed together. If dryer balls actually reduce drying time, I would expect the red markers to be mostly below the blue markers.

The scatter plot enables me to estimate how long it will take to dry clothes, which is pretty cool! If I know that a load of laundry weighs 4 kg, I can predict that the load will take about 27 minutes to dry. I would be very confident to predict that the load will dry between 24 and 34 minutes.

The drying time obviously increases with the weight of the wet clothes. However, it is not clear if the drying time increases in a straight line or maybe a curve.

Conclusions and Reflections

A limitation of my experiment is that I only used one washer and dryer. Our family's dryer is fairly new and our washer is "high efficiency," which means that it spins out a lot of the water. I have not done the experiment with an older, conventional, washer.

Based on my results, I conclude that the mean drying time does not depend on whether dryer balls are used. In particular, dryer balls do not “reduce drying time by up to 25%.” The observed difference was about 1%, but that could be due to randomness. The time required to dry clothes depends on the weight of the wet clothes, but using dryer balls does not make any difference in the time it takes to dry the clothes, whether the load is small or large.

I am very happy to report that this poster was awarded second place in the nation for the 7–9 grade level! I am excited to be chosen for this honor. I have to thank my dad because without him I never would have been motivated enough to make this poster or record and plot 40 loads of laundry.

I got presented the award at a very special ceremony that featured the 2013 president of the ASA, Marie Davidian, a professor in the Statistics Department at NC State University. Also at the ceremony was the 2012 president, Bob Rodriguez, a senior director at SAS. They gave me a plaque and I told them all about my project. They were really nice and very friendly. It really meant a lot to me that they would take time out of their busy day to show interest in my poster and to encourage me to pursue science and math in the future.

tags: Data Analysis, Just for Fun
12月 142012
 

In my previous post, I described how to implement an iterated function system (IFS) in the SAS/IML language to draw fractals. I used the famous Barnsley fern example to illustrate the technique. At the end of the article I issued a challenge: can you construct an IFS whose fractal attractor looks like a Christmas tree?

I started my attempt by "straightening" the fern example, but it ended up looking like a "fishbone fern" instead of an evergreen tree. I decided to abandon the four-function IFS and adopt a seven-function system. In addition to "left" and "right" branches, I introduced transformations for "left diagonal" and "right diagonal" branches. I also played with the relative probabilities and the translation parameters to get the right "fill" for the tree.

Eventually I constructed the following program, which implements an IFS that creates the "SAS Christmas tree" that is shown at the left. I'm not 100% satisfied with it, but as Linus says in A Charlie Brown Christmas, it's not "such a bad little tree."

proc iml;
/* For an explanation of how to construct an iterated function system in SAS, see
   http://blogs.sas.com/content/iml/2012/12/12/iterated-function-systems-and-barnsleys-fern-in-sas/
*/
/* Each row is a 2x2 linear transformation */
/* Christmas tree */
L = {0.03  0     0    0.1,
     0.85  0.00  0.00 0.85,
     0.8   0.00  0.00 0.8,
     0.2  -0.08  0.15 0.22,
    -0.2   0.08  0.15 0.22,
     0.25 -0.1   0.12 0.25,
    -0.2   0.1   0.12 0.2};
/* ... and each row is a translation vector */
B = {0 0,
     0 1.5,
     0 1.5,
     0 0.85,
     0 0.85,
     0 0.3,
     0 0.4 };
prob = { 0.02 0.6 0.1 0.07 0.07 0.07 0.07};
L = L`; B = B`; /* For convenience, transpose the L and B matrices */
 
/* Iterate the discrete stochastic map */
N = 1e5;          /* number of iterations */
x = j(2,N); k = j(N,1);
x[,1] = {0, 2};   /* initial point */
call randgen(k, "Table", prob); /* values 1-7 */
 
do i = 2 to N;
   x[,i] = shape(L[,k[i]], 2)*x[,i-1] + B[,k[i]]; /* iterate */
end;
 
/* Plot the iteration history */
y = x`;
create IFS from y[c={"x" "y"}]; append from y; close IFS;
quit;
 
/* basic IFS Christmas Tree */
ods graphics / width=200px height=400px;
proc sgplot data=IFS;
  title "SAS Christmas Tree";
  scatter x=x y=y / markerattrs=(size=1 color=ForestGreen);
  yaxis display=none;
  xaxis display=none;
run;

As I looked at the tree, however, I decided that it seemed rather bare. It seemed to "need a little love," so I decided to add some brightly colored, ornamental, decorations. I took a random sample of points on the attractor and drew balls at those points. I also added a star on top, all courtesy of PROC SGPLOT. You can download the complete program that creates the final image, which is shown at left. I used a few programming tricks, such as a data attribute map to assign customized, bright, colors for the ornaments.

So now it's your turn! In the comments, submit or link to a SAS program that creates a Christmas tree, a snowflake, or some other seasonal image. Be creative! My colleague, Robert Allison, has already posted a few whimsical SAS graphs for the holidays.

As for me, I'm finished blogging until the New Year. However, I'm not exactly taking a vacation from SAS programming because I'll be working to finish my forthcoming book, Simulating Data with SAS.

Thank you, readers, for an amazing year of blogging, interacting, and sharing ideas in 2012. Thanks to you, this blog has become the most widely-read blog at SAS!

tags: Just for Fun
12月 142012
 

In my previous post, I described how to implement an iterated function system (IFS) in the SAS/IML language to draw fractals. I used the famous Barnsley fern example to illustrate the technique. At the end of the article I issued a challenge: can you construct an IFS whose fractal attractor looks like a Christmas tree?

I started my attempt by "straightening" the fern example, but it ended up looking like a "fishbone fern" instead of an evergreen tree. I decided to abandon the four-function IFS and adopt a seven-function system. In addition to "left" and "right" branches, I introduced transformations for "left diagonal" and "right diagonal" branches. I also played with the relative probabilities and the translation parameters to get the right "fill" for the tree.

Eventually I constructed the following program, which implements an IFS that creates the "SAS Christmas tree" that is shown at the left. I'm not 100% satisfied with it, but as Linus says in A Charlie Brown Christmas, it's not "such a bad little tree."

proc iml;
/* For an explanation of how to construct an iterated function system in SAS, see
   http://blogs.sas.com/content/iml/2012/12/12/iterated-function-systems-and-barnsleys-fern-in-sas/
*/
/* Each row is a 2x2 linear transformation */
/* Christmas tree */
L = {0.03  0     0    0.1,
     0.85  0.00  0.00 0.85,
     0.8   0.00  0.00 0.8,
     0.2  -0.08  0.15 0.22,
    -0.2   0.08  0.15 0.22,
     0.25 -0.1   0.12 0.25,
    -0.2   0.1   0.12 0.2};
/* ... and each row is a translation vector */
B = {0 0,
     0 1.5,
     0 1.5,
     0 0.85,
     0 0.85,
     0 0.3,
     0 0.4 };
prob = { 0.02 0.6 0.1 0.07 0.07 0.07 0.07};
L = L`; B = B`; /* For convenience, transpose the L and B matrices */
 
/* Iterate the discrete stochastic map */
N = 1e5;          /* number of iterations */
x = j(2,N); k = j(N,1);
x[,1] = {0, 2};   /* initial point */
call randgen(k, "Table", prob); /* values 1-7 */
 
do i = 2 to N;
   x[,i] = shape(L[,k[i]], 2)*x[,i-1] + B[,k[i]]; /* iterate */
end;
 
/* Plot the iteration history */
y = x`;
create IFS from y[c={"x" "y"}]; append from y; close IFS;
quit;
 
/* basic IFS Christmas Tree */
ods graphics / width=200px height=400px;
proc sgplot data=IFS;
  title "SAS Christmas Tree";
  scatter x=x y=y / markerattrs=(size=1 color=ForestGreen);
  yaxis display=none;
  xaxis display=none;
run;

As I looked at the tree, however, I decided that it seemed rather bare. It seemed to "need a little love," so I decided to add some brightly colored, ornamental, decorations. I took a random sample of points on the attractor and drew balls at those points. I also added a star on top, all courtesy of PROC SGPLOT. You can download the complete program that creates the final image, which is shown at left. I used a few programming tricks, such as a data attribute map to assign customized, bright, colors for the ornaments.

So now it's your turn! In the comments, submit or link to a SAS program that creates a Christmas tree, a snowflake, or some other seasonal image. Be creative! My colleague, Robert Allison, has already posted a few whimsical SAS graphs for the holidays.

As for me, I'm finished blogging until the New Year. However, I'm not exactly taking a vacation from SAS programming because I'll be working to finish my forthcoming book, Simulating Data with SAS.

Thank you, readers, for an amazing year of blogging, interacting, and sharing ideas in 2012. Thanks to you, this blog has become the most widely-read blog at SAS!

tags: Just for Fun
12月 122012
 

Fractals. If you grew up in the 1980s or '90s and were interested in math and computers, chances are you played with computer generation of fractals. Who knows how many hours of computer time was spent computing Mandelbrot sets and Julia sets to ever-increasing resolutions?

When I was a kid, I realized that there were two kinds of fractals: the deterministic ones, such as the Mandelbrot set, and the probabilistic ones, such as the Barnsley fern. I spent most of my time generating the deterministic fractals, because they were simpler to understand for someone who knew little about probability or linear transformations.

Yet, I was always fascinated by how to create objects by using Michael Barnsley's iterated function systems (IFS). You can construct an IFS by defining some number of affine planar transformations, selecting one at random with some probability, and applying it to an initial point. By iteratively selecting and applying these affine transformations, the trajectory of the point traces out some fractal shape, such as the fern shown at the left. Now that I am older and know a little bit about probability and linear transformations, let's investigate how to construct Barnsley's fern in the SAS/IML matrix language. I invite DATA step programmers to create a similar program by using the SAS DATA step.

Construct Barnsley's Fern in SAS

Before I began programming the IFS, I reminded myself how to construct an iterated function system. Take a few minutes to read about the underlying mathematics. You are then ready to construct the IFS for Barnsley's fern:

  1. Create four affine transformations in the plane. Each affine transformation is of the form T(x) = A*x + b, where A is a 2x2 matrix, b is a 2x1 vector, and x is a 2x1 vector that represents a point in the plane.
  2. Specify the probability of selecting each of the four transformations.
  3. Iterate the discrete stochastic map by selecting an initial point and repeatedly doing the following a large number of times:
    • Choose one of the transformations at random, using the given (unequal) probabilities.
    • Apply the chosen transformation to the current point to obtain a new (image) point.
  4. Plot the trajectory (or orbit) of the initial point.
Under certain conditions on the transformations (they need to be contracting), the orbit of the initial point is attracted to a fractal set.

The following SAS/IML program implements an iterated function system for creating Barnsley's fern. The SGPLOT procedure is used to plot the trajectory. The main programming details are that the RANDGEN subroutine is called to generate random numbers from the "Table" distribution, and the SHAPE function is used to convert a row of coefficients into a 2x2 matrix. The SAS/IML matrix-vector operations are used to apply each affine transformation.

/*
http://en.wikipedia.org/wiki/Barnsley_fern
The Barnsley fern is created with four affine transformations:
T1(x,y) = [ 0 0   ] * [x]  with prob=0.01
          [ 0 0.16]   [y] 
 
T2(x,y) = [ 0.85 0.04] * [x] + [0  ] with prob=0.85
          [-0.04 0.85]   [y]   [1.6]
 
T3(x,y) = [ 0.2 -0.26] * [x] + [0  ] with prob=0.07
          [ 0.23 0.22]   [y]   [1.6]
 
T4(x,y) = [-0.15 0.28] * [x] + [0   ] with prob=0.07
          [ 0.26 0.24]   [y]   [0.44]
*/
proc iml;
/* 1. Each row is a 2x2 linear transformation */
L = {0    0     0    0.16,  /* L1 */
     0.85 0.04 -0.04 0.85,  /* L2 */
     0.2 -0.26  0.23 0.22,  /* L3 */
    -0.15 0.28  0.26 0.24}; /* L4 */
 
/* ... and each row is a translation vector */
B = {0 0,     /* B1 */
     0 1.6,   /* B2 */
     0 1.6,   /* B3 */
     0 0.44}; /* B4 */
/* For convenience, transpose the L and B matrices */
L = L`; B = B`;
 
/* 2. The transformations are of the form T(x) = A*x + b 
      where A and b are chosen randomly from the columns of 
      L and B. Specify probability of each transformation. */
prob = {0.01 0.85 0.07 0.07};
 
/* 3. iterate the discrete stochastic map */
N = 1e5;          /* number of iterations */
x = j(2,N);
x[,1] = {0, 2};   /* initial point */
 
do i = 2 to N;
   call randgen(k, "Table", prob);           /* k is a number 1,2,3 or 4 */
   x[,i] = shape(L[,k], 2)*x[,i-1] + B[,k];  /* apply transformation to x[i-1] */
end;
 
/* 4. plot the iteration history */
y = x`;
create IFS from y[c={"x" "y"}]; append from y; close IFS;
quit;
 
ods graphics / width=200px height=400px;
proc sgplot data=IFS;
  title "Barnsley Fern";
  scatter x=x y=y / markerattrs=(size=1);
  yaxis display=none;
  xaxis display=none;
run;

The algorithm creates the image that appears at the beginning of this article.

The holiday challenge: Construct an IFS for a Christmas tree

The cool thing about iterated function systems is that the affine transformations and the probabilities are actually parameters that you can vary. If you choose a different set of affine transformations and/or a different set of probabilities, you get a different fractal attractor. You can find Web sites that show other ferns, leaves, and branches that result from using different parameter values.

Last year during the holiday season I published some attempts to construct a Christmas tree by using the SAS DATA step. It was fun to see contributions from many SAS programmers, including a menorah! This year I issue the following challenge: can you find parameter values for an IFS that create a Christmas tree image?

I will publish my program that creates a fractal Christmas tree in a few days, so start working on your ideas! You can submit your program as a comment to my solution. Not into Christmas trees? How about a Koch snowflake instead?

tags: Just for Fun
10月 102012
 

Last week I wrote a SAS/IML program that computes the odds of winning the game of craps. I noted that the program remains valid even if the dice are not fair. For convenience, here is a SAS/IML function that computes the probability of winning at craps, given the probability vector for each of the two six-sided dice. The program is explained in the previous article.

proc iml;
start CrapsWin(dieProb);
   probs = dieProb` * dieProb; /* joint prob = product of marginals */
   events = repeat(1:6, 6) + T(1:6);
   P = j(1, 12, 0);            /* probability of each event */
   do i = 2 to ncol(P);        /* P(sum=2), ..., P(sum=12) */
      idx = loc( events=i );
      P[i] = sum( probs[idx] );
   end;
   point = (4:6) || (8:10);
   /*      Pr(7 or 11)  + Prob(Making Point) */
   Pwin  = P[7] + P[11] + sum(P[point]##2 / (P[point] + P[7]));
   return( Pwin );
finish;

I got interested in this problem when I stumbled upon a Web page by Jung-Chao Wang for a course at Western Michigan University. Wang discusses the probability of winning at craps when the game is played with unfair dice. In particular, the Web page shows how the probabilities of winning at craps change if the dice are "shaved." A "shaved die" is one that is not perfectly cubical, so that the faces with less area have a smaller probability of appearing than the faces with larger areas.

For example, you can increase the probability of rolling a 1 or 6 by sanding or trimming those two faces. The area of the 1 and 6 faces do not change, but the other four faces now have less area. Of course, you could also shave the non-1 and non-6 faces, which has the opposite effect. Wang calls this "negative shaving." (By "you," I really mean "an unscrupulous cheat." This discussion is purely a theoretical exercise in probability; I do not advocate cheating.)

For each pairs of faces (such as 1-6) you can consider the odds of winning at craps when you play with two identical dice that have been negatively or positively shaved. What is interesting to me is that the effect of shaving depends (strongly) on the faces that you shave. If you shave the 1-6 face, you get a dramatically different behavior than if you shave the 2-5 or 3-4 faces.

The results are shown in the preceding graph (click to enlarge). If you shave the 1 and 6 faces, you dramatically decrease the probability of the shooter winning. If you shave the 3 and 4 faces, you increase the probability, but the effect is not as dramatic. Presumably "the house" would want to shave the 1 and 6 sides, whereas a crooked shooter would want to try to replace the house dice with dice that are shaved on the 3 and 4 faces. However, to increase the probability of the shooter winning by, say, 5%, the 3-4 faces would need to be shaved by a lot. In contrast, to decrease the probability of the shooter winning by 5%, much less tampering is required of the 1-6 faces. Thus the "house" that cheats is less likely to be detected than the shooter that cheats!

Here's the program that creates the graph:

/* change probability of landing on 1-6, 2-5, or 3-4 sides */
/* Idea from J. C. Wang:
   http://www.stat.wmich.edu/wang/667/egs/craps2.html */
shaved = T( do( -1/12, 1/6, 0.01 ) );
WinProb = j(nrow(shaved), 3); /* results: winning probability */
do side = 1 to 3;
   do i = 1 to nrow(shaved);
      dieProb = j(1, 6, 1/6 - shaved[i] );
      dieProb[side||7-side] = 1/6 + 2*shaved[i];
      WinProb[i, side] = CrapsWin(dieProb);
   end;
end;
 
/* graph the data with PROC SGPLOT */
A = shaved || WinProb;
create Craps from A[c={"Shaved" "Side1" "Side2" "Side3"}];
append from A;  close Craps;
quit;
 
proc sgplot data=Craps;
title "Probability of Winning at Craps with Unfair Dice";
series x=Shaved y=Side1 / curvelabel="Sides 1 & 6";
series x=Shaved y=Side2 / curvelabel="Sides 2 & 5";
series x=Shaved y=Side3 / curvelabel="Sides 3 & 4";
yaxis grid values=(0.35 to 0.7 by 0.05) label="Probability of Winning";
xaxis grid label="Change in Probability";
run;

I really like this graph because it reveals how important various pairs of faces are to making a point.

For example, a pair of dice that rolls many 1s and 6s are bad for the shooter, because the most common sums will be 2, 7, and 12. You lose if 2 or 12 appear on the first roll, although you win with a 7. However, if you happen to establish a point, you are in trouble. Rolling a 2 or 12 doesn't help you make your point, so the most likely scenario is that you will roll a 7 and lose.

In contrast, a pair of dice that rolls many 3s and 4s is good for the shooter. The most common sums will be 6, 7, and 8. You win with a 7. If you establish a point, it is likely to be a 6 or 8, which are the easiest points to make.

The ideas discussed in this article also apply to other kinds of "crooked dice." The SAS/IML program opens up a variety of scenarios to simulate. For example, if only one face appears more than the others, which face would the shooter want it to be? Which face would the house want it to be? You could also modify the SAS/IML function to take two probability vectors, one for each die. Have fun simulating these scenarios, but remember: no cheating in real life!

tags: Just for Fun, Statistical Programming
10月 042012
 

Gambling games that use dice, such as the game of "craps," are often used to demonstrate the laws of probability. For two dice, the possible rolls and probability of each roll are usually represented by a matrix. Consequently, the SAS/IML language makes it easy to compute the probabilities of various events.

For two fair dice, the probability of any side appearing is 1/6. Because the rolls are independent, the joint probability is computed as the product of the marginal probabilities. You can also create a matrix that summarizes the various ways that each particular number can occur:

proc iml;
dieProb = j(1, 6, 1/6); /* each die is fair: Prob(any face)=1/6 */
/* joint prob is indep = product of marginals */
probs = dieProb` * dieProb; /* 6x6 matrix; all elements = 1/36 */
events = repeat(1:6, 6) + T(1:6);
print events[c=("1":"6") r=("1":"6")];

You can use the event and probs matrices to compute the probabilities of each possible sum of two dice. There are several ways to do this, but I like to use the LOC function to find the elements of the event matrix that correspond to each roll, and then add the associated probabilities:

P = j(1, 12, 0); /* probability of each event */
do i = 2 to ncol(P);
   idx = loc( events=i );
   P[i] = sum( probs[idx] );
end;
print P[format=6.4 c=("P1":"P12")];

The preceding table shows the probabilities of each roll. (Note that the probability of rolling a 1 with 2 dice is zero; the P[1] element is not used for any computations.) You can use the table to compute the probability of winning at craps. If you roll a 7 or 11 on the first roll, you win. If you roll a 2, 3, or 12, you lose. If you roll a 4, 5, 6, 8, 9, or 10 on the first roll (called "the point"), you continue rolling. You win if you can successfully re-roll your "point" before you roll a 7. It is an exercise in conditional probability that you can compute the probability of making a "point" in craps by summing a ratio of probabilities, as follows:

win = {7 11};
lose = {2 3 12};
point = (4:6) || (8:10);
 
Pwin1 = sum( P[win] );   /* Prob of winning on first roll */
PLose1 = sum( P[lose] ); /* Prob of losing on first roll */
/* Prob winning = Pr(7 or 11)  + Prob(Making Point) */
Pwin  = P[7] + P[11] + sum(P[point]##2 / (P[point] + P[7]));
print Pwin1 PLose1 Pwin;

According to the computation, the probability of winning at craps is almost—but not quite!—50%.

This exercise shows that having a matrix and vector language is useful for computing with probabilities. More importantly, however, this post sets the foundations for looking at the interesting case where the two dice are not fair. If you change the definition of dieProb (on the first line of the program) so that it is no longer a constant vector, all of the computations are still valid!

Next week I will post an article that shows how the odds change if some sides of the dice are more likely to appear than others. Feel free to conduct your own investigation of unfair dice before then.

tags: Just for Fun, Statistical Programming
10月 022012
 

John D. Cook posted a story about Hardy, Ramanujan, and Euler and discusses a conjecture in number theory from 1937. Cook says,

Euler discovered 635,318,657 = 158^4 + 59^4 = 134^4 + 133^4 and that this was the smallest [integer] known to be the sum of two fourth powers in two ways. It seems odd now to think of such questions being unresolved. Today we’d ask Hardy “What do you mean 635318657 is the smallest known example? Why didn’t you write a little program to find out whether it really is the smallest?”

Cook goes on to encourage his readers to write a program that settles the conjecture, which is obviously no longer an open problem. Here's a SAS/IML solution. First, introduce a function that I use a lot:

proc iml;
start ExpandGrid( x, y ); 
/* Return ordered pairs on a uniform grid of points */
   Nx = nrow(x);    Ny = nrow(y);
   x = repeat(x, Ny);
   y = shape( repeat(y, 1, Nx), 0, 1 );
   return ( x || y );
finish;

The ExpandGrid function generates a uniform grid of all ordered pairs of it's arguments. To solve the "conjecture," generate all ordered pairs that contain the integers 1–158 and check to see if the sum of fourth powers contain duplicate values. The number 158 is the greatest integer that is less than the fourth root of 635,318,657.

/* settle conjecture of Euler/Hardy: Is 635,318,657 the smallest
   integer that can be written as the sum of 4th powers in two ways? */
k = floor(635318657##0.25);        /* only need to test up to 4th root */
g = ExpandGrid(t(1:k), t(1:k));    /* all ordered pairs */
g = g[ loc(g[,1] <= g[,2]), ];     /* omit symmetric pairs */
f = g[,1]##4 + g[,2]##4;           /* sum of 4th powers */
print (nrow(f)) (ncol(unique(f))); /* exactly one sum is the same */

The output proves the conjecture: there are 12,561 ordered pairs of integers to consider, and there are 12,560 unique sums of fourth powers. That means that there is exactly one sum of fourth powers that can be written in two different ways, and Euler has already told us the particular values.

If you didn't know Euler's results, you could use the same computation to find Euler's values: sort the numbers by the sum of fourth powers and use the DIF function to take differences between consecutive elements. Print out the rows for which the difference is zero (that is, that the sum of fourth powers are equal):

/* find the actual pairs of integers that have the same sum */
g = g || f;                    /* concatenate fourth powers onto pairs */
call sort(g, 3);               /* sort by sum of 4th powers */
idx = loc(dif(g[,3])=0);       /* dif=0 <==> consecutive values equal */
print (g[idx-1:idx, ])[c={N1 N2 "Sum"}];

It is fun to think of what someone like Euler might have accomplished if had the tools that we do today!

tags: Just for Fun, Numerical Analysis
8月 292012
 

Magic squares are cool. Algorithms that create magic squares are even cooler.

You probably remember magic squares from your childhood: they are n x n matrices that contain the numbers 1,2,...,n2 and for which the row sum, column sum, and the sum of both diagonals are the same value. There are many algorithms to generate magic squares. I remember learning as a child how to construct a magic square for any odd number, n, by using the "Siamese method." I also remember being fascinated by Ben Franklin's construction of "semi-magic squares" in which he used the sum of "bent diagonals" instead of straight diagonals.

Last week I had email discussions with several readers who had the idea to generate magic squares in SAS. The readers noted that my recent blog post on "How to return multiple values from a SAS/IML function," which features a module that computes row sums and column sums of a matrix, could be extended to check whether a matrix is a magic square.

When I first came to SAS, I was disappointed because there is no function in the SAS/IML language that generates magic squares. Prior to arriving at SAS, I had used MATLAB for matrix computations, and the MATLAB documentation is full of calls to the magic function, which generates magic squares of any size. An internet search for "magic squares in SAS" did not yield any useful results. However, when I searched for "MATLAB magic squares," I was thrilled to discover that Cleve Moler—cofounder of The MathWorks and numerical analyst extraordinaire—has written an e-book (C. Moler, 2011, Experiments with MATLAB) that includes a chapter on the construction of magic squares!

The following program implements Moler's algorithm in the SAS/IML language. For the description of the algorithm, see Chapter 10 (p. 7–10) of Moler's e-book.

proc iml;
 
start Magic(n);
/* Moler's algorithm for constructing a magic square for any n.
   See Moler, C. (2011) "Magic Square," Chapter 10 in 
   Experiments with MATLAB, p. 7-10
   E-Book: http://www.mathworks.com/moler/exm/chapters/magic.pdf
*/
 
   /* define helper functions */
   start ModPos(a, q); /* always return positive value for mod(a,q) */
      return( a - q*floor(a/q) );
   finish;
 
   /* Magic square when n is odd */
   start MagicOdd(n);
      I = repeat( T(1:n), 1, n );
      J = repeat( 1:n, n );  /* or T(I) */
      A = ModPos(I+J-(n+3)/2, n); /* modify formula to match MATLAB output */
      B = mod(I+2*J-2, n);
      return( n*A + B + 1 );
   finish;
 
   /* Magic square when n is a multiple of 4 */
   start Magic4(n);
      M = shape(1:n##2, n, n);
      I = repeat( T(1:n), 1, n );
      J = repeat( 1:n, n );  /* or T(I) */
      idx = loc(floor(mod(I,4)/2) = floor(mod(J,4)/2));
      M[idx] = n##2 + 1 - M[idx];
      return( M );
   finish;
 
   /* Magic square when n>4 is even but not a multiple of 4 */
   start Magic2(n);
      s = n/2;           /* s is odd */
      A = MagicOdd(s);
      M = (A         || (A + 2*s##2)) // /* 2x2 blocks of magic squares */
         ((A+3*s##2) || (A +   s##2));
      /* col sums are correct; permute within cols to adjust row sums */
      i = 1:s;
      r = (n-2)/4;
      h = do(n-r+2,n,1); /* might be empty */
      j = (1:r) || h;
      M[i || i+s, j] = M[i+s || i, j]; /* swap rows for cols in j */
      i = r+1;
      j = 1 || i;
      M[i || i+s, j] = M[i+s || i, j]; /* swap rows for cols in j */
      return( M );
   finish;
 
   /* main algorithm */
   if n<=1 then return(1);
   if n=2 then return( {. ., . .} ); /* no magic square exists */
   if mod(n,2)=1 then return( MagicOdd(n) ); /* n is odd */
   if mod(n,4)=0 then return( Magic4(n) );   /* n is divisible by 4 */
   return( Magic2(n) );        /* n is divisible by 2, but not by 4 */
finish;

The implementation uses two interesting features. First, I define a helper function that always returns a positive value for the expression mod(a,q), because the MOD function in SAS can sometimes return a negative value. Secondly, the computation for the Magic2 function constructs a block matrix of size 2n from a magic square of size n. Ian Wakeling told me that this is an application of a general technique to construct a large magic square from smaller components.

To test the implementation, you can construct and print a few magic squares, as follows:

 
do n = 3 to 6;
   M = Magic(n);
   L = cat("Magic Square, n=", char(n,1)); /* form label */
   print M[label=L];
end;

At long last, I have my wish: a function in PROC IML that generates magic squares!

Do you have a favorite memory regarding magic squares? Do you have a favorite algorithm that constructs a magic square? Post a comment!

tags: Just for Fun, Matrix Computations
6月 182012
 

To celebrate special occasions like Father's Day, I like to relax with a cup of coffee and read the newspaper. When I looked at the weather page, I was astonished by the seeming uniformity of temperatures across the contiguous US. The weather map in my newspaper was almost entirely yellow and orange, which told me that the midday temperatures were predicted to have extremely small variation. It looked like the range of temperatures might vary by less than 30 degrees Fahrenheit (17 degrees Celsius), with most temperatures in the 70's and 80's.

To me this seemed like an extremely small variation, so that afternoon I downloaded the actual midday temperatures from Weather.com, which are shown below:

The Weather.com map uses a different color scale than my local paper. Nevertheless I quickly estimated that most of the temperatures were between 62 F and 92 F, with just a few hot spots in southern Texas, Arizona, and Nevada. The statistician in me wondered, what is the record for the "most uniform" temperatures at midday in the contiguous US? I'll bet that someone can automate the downloading of temperatures for standard US locations and answer the question. I'd love to see a graph of "temperature variation vs. date" for, say, the last ten or twenty years.

I don't know of a source for temperature data, but I figured that I might as well do SOME statistics to celebrate Father's Day weekend! I enlisted my youngest child to assist me in typing in the temperatures. (Yes, it was a father/daughter bonding moment.) I ran the following univariate analysis of the temperatures:

data Temps06162012;
input Temperature @@;
datalines;
66 72 72 73 86 81 68 67 70 64
64 67 63 73 71 75 81 75 94 96 77 96 93
67 71 72 70 74 69 86 80 83 78 86
63 69 71 71 79 78 82 83 85 84 89 90 88 89 92 94
79 66 79 74 72 75 88 83 88 82 91 88 89 90 86 86 83 84
79 81 85 86 80 84 84 77 85 74 79 83 84 85
62 63 81 61 78 78 81 80 74 78 81 83 84 88 87
;
 
proc univariate data=Temps06162012;
   var Temperature;
   histogram Temperature;
run;

The histogram and the accompanying tabular output told me that the mean temperature was a comfortable 79 degrees F, and that the standard deviations of temperatures was a paltry 8.56 degrees F. It was, indeed, a beautiful weekend for celebrating Father's Day.

Can anyone find a date for which the temperatures exhibited less variation? I'd love to see the analysis.

tags: Data Analysis, Just for Fun
6月 132012
 

A collegue who works with time series sent me the following code snippet. He said that the calculation was overflowing and wanted to know if this was a bug in SAS:

data A(drop=m);
call streaminit(12345);
m = 2;
x = 0;
do i = 1 to 5000;
   x = m*x + rand("Normal");
   output;
end;
run;
NOTE: A number has become too large at line 6 column 11.
      The number is >1.80E+308 or <-1.80E+308.

This behavior is expected, and is definitely not a bug. The behavior follows from a basic analysis of the computation within the DO loop. Take a moment to figure out why the magnitude of x becomes large in this program.

A simple discrete dynamical system

Before I became a statistician, I studied dynamical systems, which is the study of the long-term qualitative behavior of a differential equation or of an iterative map. A lot can be learned about the behavior of iterates of x by just ignoring the random term. (You can justify ignoring the random term by noting that rand("Normal") has zero mean.)

The iterative mapping x → mx is linear, so it has simple dynamical properties:

  • If m < 1, all initial values for x approach 0 under iteration because the sequence x, mx, m2x, ..., mnx approaches zero for all values of x.
  • If m > 1, all nonzero values of x diverge under iteration.

A discrete dynamical system with a stochastic component

Adding a random term with mean zero does not change the mathematics of this problem very much. However, now the dynamics are described in terms of "expected behavior," because the system is no longer deterministic. For the iterative mapping x → mx + ε, ε ∼ N(0,1), you can expect the following long-term behavior for the iterates:

  • the iterates to bounce around a neighborhood of zero when m < 1
  • the iterates diverge for m > 1

The following graphs show typical behavior:

data A;
call streaminit(12345);
phi = 0.99; x = 0;
do i = 1 to 5000;
   x = phi*x + rand("Normal"); output;
end;
phi = 1.01; x = 0;
do i = 1 to 5000;
   x = phi*x + rand("Normal"); output;
end;
run;
 
proc sgpanel data=A;
panelby phi / uniscale=column rows=2;
refline 0 / axis=y;
scatter x=i y=x;
run;

You might wonder about the case m = 1. There is a theorem that says that a random walk x → x + ε, ε ∼ N(0,1), will eventually get arbitrarily far from it's initial position. Therefore the case m = 1 would eventually diverge in finite precision arithmetic...but you might have to wait a REALLY long time!

tags: Just for Fun, Numerical Analysis, Statistical Thinking