Just for Fun

3月 142022
 

On this Pi Day, let's explore the "πth roots of unity." (Pi Day is celebrated in the US on 3/14 to celebrate π ≈ 3.14159....) It's okay if you've never heard of the πth roots of unity. This article starts by reviewing the better-known nth roots of unity. It then shows all the math you'll need to understand the πth roots of unity and how to visualize them!

The nth roots of unity

Recall that the nth root of unity for a positive integer, n, is a complex number, z, such that the nth power of z equals 1. In symbols, z satisfies the equation zn = 1.

For any positive integer, n, there are n roots of unity. The nth roots of unity are the complex numbers
\(\exp(2\pi i k /n) = \cos(2 \pi k/n) + i \sin(2 \pi k/n)\)
for k=0, 1, 2, ..., n. Geometrically, these points are the equally spaced points on the unit circle. You get them by tracing the image of the point 1 + 0*i as you rotate the circle by 2π/n radians. Some simple examples are:

  • n=2. The square roots of unity are the numbers 1 and -1. Geometrically, these points are represented in the Argand plane by the ordered pairs (1, 0) and (-1, 0). They are the orbit of the point (1,0) under rotations by π radians.
  • n=3. The cube roots of unity are the numbers 1, -1/2 + i sqrt(3)/2, and -1/2 - i sqrt(3)/2. These points are the orbit of the point (1,0) under rotations by 2π/3 radians.
  • n=4. The 4th roots of unity are the numbers 1, i, -1, and -i. The coordinates of these points on the unit circle are (1,0), (0, 1), (-1, 0), and (0, -1). These points are the orbit of the point (1,0) under rotations by π/2 radians.

The πth roots of unity

So, what is a πth root of unity? It is complex number, z, such that zπ = 1. Geometrically, it is one of the points on the orbit of the point (1,0) by rotations of 2*pi/pi = 2 radians. These complex numbers have the formula \(z = \exp(2 i n)\) for any integer, n. You can verify that these numbers are the πth roots of unity by showing that \(z^\pi = 1\), as follows:
\(\exp(2 i n)^\pi = \exp(2\pi i n) = \cos(2 \pi n) + i \sin(2 \pi n) = 1\) for all integers n.

Notice that the "circle map" x → exp(i x) sends the even integers to the πth roots of unity. Thus, you can also refer to the points as the image of the even integers under the mapping.

Geometrically, the πth roots are generated by rotating the unit circle by 2 radians or by -2 radians. For simplicity, the subsequent sections visualize the πth roots for only the positive integers, which correspond to counterclockwise rotations by 2 radians.

Visualize the πth roots of unity

Let's generate some positive integers and plot the complex points exp(2 i n) in the Argand plane. To get started, let's plot only the images of the first 25 even integers:

%let maxN = 355;
data Roots;
do n = 1 to &maxN;
   s = 2*n;
   xn = cos(s);   /* image of even integers under mapping to unit circle */
   yn = sin(s);
   output;
end;
run;
 
/* define a helpful macro */
%macro PlotCircle();
ellipseparm semimajor=1 semiminor=1 / xorigin=0 yorigin=0 outline lineattrs=(color=cxCCCCCC);
%mend;
 
ods graphics / width=400px height=400px;
title "Wrap Positive Even Integers onto a Circle";
title2 "First 25 Even Integers";
proc sgplot data=Roots(obs=25) aspect=1 noautolegend;
   %PlotCircle;
   scatter x=xn y=yn / datalabel=s;
   refline 0 / axis=x;   refline 0 / axis=y;
   xaxis display=none;   yaxis display=none;
run;

Geometrically, you can think of wrapping the real number line onto the circle like thread onto a spool. The points show the images of the even integers. You can see that 2 is mapped into the second quadrant, 4 is mapped to the third quadrant, and 6 is mapped to the fourth quadrant just below the point (1,0) because 6 radians is slightly less than 2π. This pattern repeats for the images of 8, 10, and 12. Each set of three points slightly lags the previous set of three points.

The image of 22 is very close to (-1, 0). This is because 22/7 is a close approximation of π, so 22 ≈ 7π. Similarly, the image of 44 ≈ 14π is close to the point (1, 0).

The first 22 images are nearly equally spaced on the circle, but not quite. You can see that the images of 46, 48, and 50 almost overlay the images of 2, 4, and 6, but are slightly separated. This pattern continues for the next set of 22 even integers, and the next, and the next. No image ever hits a previous image because we are rotating by 2 radians and there are 2π radians in a circle. (These numbers are not commensurable.) The following graph shows the image of a few hundred even integers onto the circle:

This graph visualizes the πth roots of unity. The images of all positive and all negative integers will never overlap or repeat. By the irrational rotation theorem (rescaled), the πth roots of unity are dense in the circle.

The Beatty sequence on the circle

In a previous article, I showed that every positive integer belongs to one of two disjoint sets. Each positive integer belongs to the Beatty sequence for π, or it belongs to the complementary sequence. The same is true for the even positive integers: some are elements of the Beatty sequence for π, and the rest are not. Let's see what happens if we color each point on the circle according to whether it is the image of the Beatty sequence:

Did you expect this result? Wow! The graph shows that all of the integers that belong to the Beatty sequence are mapped into two arcs on the unit circle. The red arc in the second quadrant is for radian angles in the interval (π-1, π). The red arc in the fourth quadrant is for radian angles in the interval (2π-1, 2π). These values can be deduced by using the definition of the Beatty sequence, and the fact that a Beatty number can be represented as floor(k π) = k π - f for some fractional part f in the interval (0, 1).

The graph shows that (1/π)th of the πth roots of unity belong to the Beatty sequence for π. And these numbers are all clumped together into two arcs on the circle.

Technically, the Beatty sequence is defined only for positive integers. Thus, I've illustrated this result for positive integers. However, you can extend the Beatty sequence to all non-zero integers in a natural way. If you do that, then the roots of unity that belong to the "negative Beatty sequence" are visualized by reflecting the circle across the horizontal axis. If you overlay the two circles, the extended red and blue markers interlace, and the result is not as striking.

Summary

There are infinitely many πth roots of unity. They are complex numbers of the form \(\exp(2 i n)\) for any integer, n. Geometrically, you can visualize these points as the rotations by ±2 radians of the point (1,0) in the Argand plane. The πth roots of unity are dense in the unit circle. If you restrict your attention to only positive rotations (the images of the positive integers), you can classify each root according to whether it is the image of an even integer in the Beatty sequence for π or in the complementary sequence. The points that are in the Beatty sequence are mapped to the arcs for angles (π-1, π) and (2π-1, 2π). These arcs are (1/π)th of the arclength of the unit circle.

You can download the SAS program that creates these images.

Further reading

The post The pi_th roots of unity appeared first on The DO Loop.

2月 092022
 

For some reason, SAS programmers like to express their love by writing SAS programs. Since Valentine's Day is next week, I thought I would add another SAS graphic to the collection of ways to use SAS to express your love.

Last week, I showed how to use vector operation and a root-finding algorithm to visualize the path of a ball rolling on an elliptical billiards table. Every time the ball hits the boundary of the table, it rebounds and continues rolling in a new direction. At the end of the article, I commented that "the same program can compute the path of a ball on any table" for which certain conditions hold. With Valentine's Day coming soon, I wondered whether I could visualize the dynamics of a ball bouncing across a heart-shaped table. The graph to the right is the result of my efforts.

In the figure, the pink area is the interior of the heart-shaped billiard table. The red lines trace the path of a ball that starts in the middle of the heart with an initial motion at 45 degrees to the horizontal. The figure shows the ball's path for the first 250 bounces.

Although the caption of the figure is "Love Makes My Heart Bounce," the figure might also resonate with people who feel like they bounce around and never find love. It might appeal to those who have recently ended a relationship and feel like they are on the rebound. And during the global pandemic, maybe you have stayed close to home, never venturing very far from those you love? This heart is for you, too!

Details of the computation

Here are a few details about the computation:

  • The shape of the heart is defined by \(H(x,y)=0\), where \(H(x,y) = (x^2 + y^2 - 1)^3 - x^2 y^3\).
  • The gradient of H vanishes at the four points on the curve, namely (0, ±1) and (±1, 0). If the virtual ball "hits" the table's virtual boundary at those locations, the ball will "stick" to the boundary because the program uses the gradient of the function to compute the normal vector to the curve. The program could be modified to compute the normal vector at (±1, 0). However, the curve has a cusp singularity at (0, ±1), so it is impossible to define a normal vector at those points. On a real billiard table, you would want to round off the singularity at those points.
  • Because the heart shape is not convex, I had to modify the routine that finds the intersection of a ray and the curve H=0. The previous article on this topic assumes that a ray intersects the curve at only one point. Along the top of the heart, a ray might intersect the curve at three points. I modified the program to find the nearest intersection point.

You can download the program that I used to make the graph. Feel free to change the colors and caption and spread the love.

The post Billiards on a heart-shaped table appeared first on The DO Loop.

2月 072022
 

I recently showed how to find the intersection between a line and a circle. While working on the problem, I was reminded of a fun mathematical game. Suppose you make a billiard table in the shape of a circle or an ellipse. What is the path for a ball at position p on the table if it rolls in a specified direction, as determined by a vector, v? A possible path is shown in the graph at the right for a circular table. The figure shows the first 50 bounces for a ball that starts at position p=(0.5 -0.5) and travels in the direction of v=(1, 1).

This article shows how to use SAS to compute images like this for circular and elliptical tables.

Vectors and geometry

One of the advantages of the SAS/IML matrix language is its natural syntax: You can often translate an algorithm from its mathematical formulation to a computer program by using only a few statements. This section presents the mathematical formulation for finding the trajectories of a ball on a circular or elliptical billiard table, then shows how to translate the math into SAS/IML statements.

The fundamental mathematical assumption is that all collisions are perfectly elastic and that the angle of incidence equals the angle of reflection. For a curved surface, the angles are measured by using the normal and tangent vectors to the curve. The geometry is shown in the following figure.

A ball moves toward the boundary of the table with velocity vector v. The boundary is represented as the level set H = 0 for some bivariate function, H(x). For example, a circular table uses H(x,y) = x2 + y2 - 1. The previous article showed how to find the point at which the ball hits the curve.

You can compute the unit normal vector, N, at that point by using the gradient of H. Project the velocity vector onto the unit normal to obtain vN. The linear space orthogonal to N is the tangent space. The vector vT = v - vN is the projection of v onto the tangent space. (Note that v = vT + vN.) The reflection of v is defined to be the vector that has the same tangent components but the opposite normal components, which is vrefl = vT - vN. A little algebra shows that vrefl = v - 2vN, so, in practice, the projection onto the tangent space is not needed.

The following SAS/IML program implements functions that compute these vectors. To test the functions, assume a ball is initially at the point p=(0.5, -0.5) and moves with unit velocity v=(1, 1)/sqrt(2), which is 45 degrees from the horizontal.

proc iml;
/* Define the multivariate function H and its gradient.
   Example: H(x)=0 defines a circle and grad(H) is the gradient */
start H(x);               
   return( x[ ,1]##2 + x[ ,2]##2 - 1 );  /* = x##2 + y##2 - 1 */
finish;
start gradH(x);
   return( 2*x );                        /* { dH/dx  dH/dy } */
finish;
 
/* given x such that H(x)=0 and a unit vector, v, find the 
   normal and tangent components of v such that v = v_N + v_T, 
   where v_N is normal to H=0 and v_T is tangent to H=0 */
/* project v onto the normal space at x */
start NormalH(x, v);
   grad = gradH(x);      /* if H(x)=0, grad(H) is normal to H=0 at x */
   N = grad / norm(grad);/* unit normal vector to curve at x */
   v_N = (v*N`) * N;     /* projection of v onto normal vetor */
   return( v_N ); 
finish;
/* project v onto the tangent space at x */
start TangentH(x, v);     
   v_N = NormalH(x, v);
   v_T = v - v_N;
   return( v_T ); 
finish;
/* reflect v across tangent space at x */
start Reflect(x, v);
   v_N = NormalH(x, v);
   v_refl = v - 2*v_N;
   return( v_refl );
finish;
 
/* Test the definitions for the following points and direction vector */
p = {0.5 -0.5};           /* initial position of ball */
v = {1 1} / sqrt(2);      /* unit velocity vector: sqrt(2) = 0.7071068 */
q = {1 0};                /* ball hits circle at this point */
N = NormalH(q, v);        /* normal component to curve at q */
T = TangentH(q, v);       /* tangential component to curve at q */
v_refl = Reflect(q, v);   /* reflection of v */
print v, N, T, v_refl;

From the initial position and velocity, it is clear that the ball will hit the circle at q=(1,0) and reflect off at 135 degrees, which is the direction vector vrefl=(-1, 1)/sqrt(2). For this geometry, the unit normal vector at q is N=(1,0) and the unit tangent vector is T=(0,1), so it is easy to mentally verify that the program is giving the correct values for the vectors in this problem.

Bouncing around a billiard table

Let's play a game of billiards on a table that has a nonstandard shape. To begin the game, you must specify the initial position (p) and velocity (v) of the ball, as well as the shape of the table (a function H). The algorithm is simple:

  1. Given p and v, find the point q where the ball hits the boundary of the table. This step is described in a previous article..
  2. Reflect v across the normal vector at q to obtain a new direction, vrefl.
  3. Repeat this process some number of times. Optionally, plot the path to visualize how the ball bounces across the table.

The following SAS/IML program implements this algorithm by using the F_Line and FindIntersection functions from the previous article:

/* To find the intersection of H=0 with the line through G_p in the direction of G_v, see
   https://blogs.sas.com/content/iml/2022/02/02/line-search-sas.html */
/* Restriction of H to the line through G_p in the direction of G_v */
start F_Line(t) global(G_p, G_v);
  return( H( G_p + t*G_v ) );
finish;
/* return the intersection of H=0 and the line */
start FindIntersection(p, v, step=20) global(G_p, G_v);
   G_p = p; G_v = v;
   t = froot("F_Line", 1E-3 || step);  /* find F(t) = 0 for t in [1E-3, step] */
   q = p + t*v;                        /* H(q) = 0 */
   return( q );
finish;
/* given a point, p, and a direction vector, v, find the point
   q on the line through p in the direction of v such that H(q)=0.
   Reflect the velocity vector to get a new direction vector, v_refl */
start HitTable(q, v_refl, p, v);
   q = FindIntersection(p, v);
   v_refl = Reflect(q, v);
finish;
 
/* Example: A periodic path with four segments */
p = {0.5 -0.5};     /* point p */
v = {1 1}/sqrt(2);  /* direction vector v */
 
NSeg = 6;           /* follow for this many bounces */
x = j(NSeg, ncol(p));
x[1,] = p;
do i = 2 to NSeg;
   run HitTable( q, v, x[i-1,], v ); /* overwrite v on each iteration */
   x[i,] = q;
end;
print x[L="Positions"];

These initial conditions are highly symmetric and therefore easy to visualize. From the initial position and velocity, the ball hits the boundary of the circular table at (1,0), (0,1), (-1,0), and (0,-1) before returning again to (1,0). The ball will follow this periodic trajectory forever. It will always hit the same four points on the boundary of the circle.

For a general initial condition, you will want to draw a graph that visualizes the orbit of the ball as it bounces across the table. You can download the SAS program that computes the graphs in this article. For ease of use, I encapsulated the computations and visualization into a module called Billiards. You can call the Billiards subroutine with an initial condition and the number of bounces and the module will return a graph of the trajectory of the ball. For example, the following statements generate the first 50 bounces of a ball on a circular table:

/* Example 2. Either periodic orbit or bounces inside an annulus */
p = {0.5 -0.5};     /* point p */
v = {1 1.4};        /* vector v */
v = v / norm(v);    /* standardize to unit vector */
print v;
 
title "50 Bounces on a Round Table";
title2 "p=(0.5, -0.5); v=(0.58, 0.81)";
run Billiards(p, v, 50);  /* NSeg=50 bounces */

The graph is shown at the top of this article. For this example, the ball bounces around the table inside an annular region. For a circular table, you can prove mathematically that every trajectory is either a periodic orbit or is dense in an annular region.

Billiards on an elliptical table

If you study the previous program, you will see that the program does not use any special knowledge of the shape of the table. Given any smooth functions H and grad(H), the program should be able to handle tables of other shapes. That is the advantage of writing a program that uses vectors in a coordinate-free manner.

Let's change the functions for H and grad(H) to use an elliptical table with parameters a and b. That is, the shape of the table is determined by H(x)=0, where \(H(\mathbf{x}) = x_1^2/a^2 + x_2^2/b^2 - 1\) is the standard form of an ellipse. For the next example, let's use a=sqrt(5) and b=2.

/* define the elliptical region and gradient vector */
start H(x); 
   a = sqrt(5); b = 2;
   return( (x[ ,1]/a)##2 + (x[ ,2]/b)##2 - 1 ); /* ellipse(a,b) */
finish;
start gradH(x);
   a = sqrt(5); b = 2;
   return( 2*x[ ,1]/a##2 || 2*x[ ,2]/b##2 ); 
finish;

An ellipse is determined by its two foci, which are located at (+/-sqrt(a2 - b2), 0) when a > b. For our example, the foci are at (-1, 0) and (1, 0).

You can show that the path of a ball on an elliptical table is one of three shapes:

  • If the ball does not pass between the foci, its path is inside an annular region.
  • If the ball passes between the foci, its path is inside a hyperbolic region.
  • If the ball passes over one focal point, it also passes over the second focal point. The path converges to the major axis that contains the foci.

The following examples show each case:

/* If the ball does not pass between the foci, its path is inside an annular region. */
p = {1.5 0};        /* point p */
v = {0   1};        /* vector v */
title "250 Bounces on an Elliptical Table";
title2 "Ball Does Not Pass Between Foci";
run Billiards(p, v, 250);
/* If the ball passes between the foci, its path is inside a hyperbolic region. */
p = {0  0};         /* point p */
v = {1  1}/sqrt(2); /* vector v */
title "250 Bounces on an Elliptical Table";
title2 "Ball Passes Between Foci";
run Billiards(p, v, 250);
/* If the ball passes over one focal point, it also passes over the second focal point. 
   The path converges to the major axis that contains the foci. */
p = {1 -1};         /* point p */
v = {0  1};         /* vector v */
title "10 Bounces on an Elliptical Table";
title2 "Ball Passes Through a Focal Point";
run Billiards(p, v, 10);

Summary

In summary, you can use analytical geometry to figure out the trajectory of a ball on nonstandard billiard table that has a curved boundary. By using vectors and a root-finding algorithm in the SAS/IML language, you can write a program that displays the path of a ball on a circular table. The advantage of a coordinate-free vector formulation is that the same program can compute the path of a ball on a table that has a more complicated boundary, such as an ellipse. In fact, the same program can compute the path of a ball on any table for which H=0 defines a smooth convex curve on which the gradient of H does not vanish.

Incidentally, people have built elliptical billiard tables and have put a hole at one of the foci. To sink a ball, you can shoot it at the hole or you can shoot it at the other focal point. Either way, the ball rolls into the hole.

You can download the SAS program that I used to create the graphs in this article.

The post Billiards on an elliptical table appeared first on The DO Loop.

12月 152021
 

Suppose you are creating a craft project for the Christmas holidays, and you want to choose a palette of Christmas colors to give it a cheery holiday appearance. You could use one of the many online collections of color palettes to choose a Christmas-themed palette. However, I didn't want to merely take a palette that someone else designed. As a statistician, I wanted to create my own Christmas palette by using some sort of mathematical technique that creates a "statistically perfect" set of Christmas colors from other palettes!

Artists use their talents to create aesthetically pleasing and harmonious palettes of colors. But why should artists have all the fun? In this article, I show how a statistician might create a palette of Christmas colors! Starting from existing palettes of colors, the technique uses principal component analysis to select a linear subspace of colors. From within that subspace, you can create your own palette of colors. An example is the green-gold-red palette that is shown below. It was made by math, but I think it is pretty!

The same technique could also be used to create autumn palettes, spring palettes, and so forth, but don't take this too seriously: this article is intended to be geeky fun for the holidays. I'll let others debate whether this mathematical technique is useful for creating art.

Outline of the method

Here is the general idea behind creating a statistical palette of colors:

  1. Collect 6-10 existing palettes that you like. This process should produce 40-60 colors that are related to your theme (such as "Christmas"). I have written about how to read and visualize palettes in SAS.
  2. Each color can be represented a triplet of values in RGB space. You can use principal component (PC) analysis to project the colors into a plane that captures most of the variation in the colors. The graph at the right shows the principal component scores for 41 colors in seven Christmas palettes.
  3. Each palette is represented as a path through this linear subspace. So, you can choose a new sequence of points in this subspace and then REVERSE the projection to obtain a new palette of RGB colors.
  4. There are infinitely many paths you could choose, but one choice is to use vertices of a regular polygon to create a "statistically perfect" set of Christmas colors.

Paths through principal component space

For each palette, if you connect the sequence of colors in the palette, you obtain a path through the plane of the principal component scores. If you try to display the paths for all seven palettes, you will get a jumbled mess of lines that crisscross the graph. To make it easier to visualize the paths for each palette, the following graph limits the display to two palettes. The "Jingle Bell" palette has six colors. It starts with two green colors, moves to golden colors, and ends with two red colors. The "Unwrapping My Gifts" palette has five colors. It starts with two green colors, moves to a silver color, and ends with a red and dark burgundy color. The other palettes trace similar paths through the plane of scores.

Create your own path

If each palette corresponds to a path, then each path corresponds to a palette. In other words, if you choose a sequence five or six coordinates in the space of principal component scores, the path defines a new Christmas palette.

The ancient Greeks thought that a circle was the most perfect shape. Regular polygons, which approximate the circle, are also seen as perfect in their symmetry and beauty. Therefore, it makes sense to use a regular hexagon to create a new Christmas palette. The following graph demonstrates how you can choose a new Christmas palette. The circle represents a path through the space of colors. The gray stars represent six points on the circle, which are also vertices on a regular hexagon. I will define a new palette by using math to figure out what colors to associate with the coordinates of the six stars. I've ordered the six points starting from the point that I know will be green (because it is in the lower-left corner) and ending with the point that I know will be red (because it is in the upper-left corner).

Inverting the principal component transformation

Until now, we've only considered the projection of the RGB colors onto the plane of the first two PCs. Of course, this projection is not invertible: There are many colors that share the same two PC scores. However, in general, the transformation from RGB space to the space of all principal component scores is invertible. You can represent the transformation as a matrix by using the spectral decomposition of the correlation matrix of the colors. To embed the plane of the first two PCs, you need to set a value for the coordinate of the third PC. The obvious statistical choice is to use PC3=0, which represents the average score in the third PC. In other words, you can associate the first two PC scores with the plane that has ordered triplets (PC1, PC2, 0).

With that choice, it is straightforward to apply the inverse PC transformation. The inverse transformation maps PC scores into colors in RGB space. You can use the SAS IML language to carry out this transformation. If you visualize a dense set of points on the circle, you get the following spectrum of colors:

Any smaller palette will contain a subset of these colors. For example, if you choose six points on the circle that are vertices of a regular histogram, you get the palette that is shown at the top of this article. The new palette is based on the original seven palettes but is distinct from them. The new palette incorporates the properties of the original palettes but is defined by defining a new path through the PC scores.

Summary

In summary, this article shows how to create a new palette of colors that is based on statistical properties of a set of colors. You can perform a principal component analysis in RGB space to obtain a plane of PC scores. Every path through that plane corresponds to a new palette of colors. For my palette, I chose a sequence of vertices on a regular hexagon.

This method is not guaranteed to give an aesthetically pleasing palette of colors, but I think it does in this example. And if you don't like my palette, you can choose any other path through the set of PC scores to create a new palette! If you would like to experiment with this technique, you can download the SAS program that I used to generate the figures in this article.

The post A statistical palette of Christmas colors appeared first on The DO Loop.

8月 042021
 

In a previous article, I discussed a beautiful painting called "Phantom’s Shadow, 2018" by the Nigerian-born artist, Odili Donald Odita. I noted that if you overlay a 4 x 4 grid on the painting, then each cell contains a four-bladed pinwheel shape. The cells display rotations and reflections of the pinwheel.

The figure at the right overlays the grid on Odita's painting. The cell in the upper right (marked R0) displays one pinwheel. When I first saw Odita's artwork, I thought that all 16 cells were a rotation or reflection of this one pinwheel. However, it turns out that Odita's work contains a mathematical twist. To generate all the pinwheels in the image, you need to use two different colorings of the pinwheel. The first coloring is (in counter-clockwise order) teal, orange, blue, and salmon. The second coloring is teal, salmon, blue, and orange. You cannot get the second colored pinwheel by rotating and reflecting the first.

Rotations and reflections of a pinwheel

The "fundamental pinwheel" (R0) is displayed to the left. As I showed in the previous article, you can generate this pinwheel by starting with a four-sided polygon and rotating it around the origin.

For simplicity, you can start with the coordinates of the pinwheel. The following SAS data set contains the coordinates of the pinwheel vertices. As shown in the previous article, you can use the POLYGON statement in PROC SGPLOT to visualize the pinwheel.

data Shape;
input ID x y @@;
datalines;
0  0 0  0  0.667  0.333  0  1  1  0  0  1  0  0  0  
1  0 0  1 -0.333  0.667  1 -1  1  1 -1  0  1  0  0 
2  0 0  2 -0.667 -0.333  2 -1 -1  2  0 -1  2  0  0 
3  0 0  3  0.333 -0.667  3  1 -1  3  1  0  3  0  0 
;

If you ignore the colors, the pinwheel shape has rotational and reflectional symmetries. Let's understand how this pinwheel looks when it is rotated or reflected under the symmetries of the square. The symmetries make up a finite group of order eight, which is known as the dihedral group of the square or D4. The group has a representation in terms of eight 2 x 2 orthogonal matrices. The following SAS/IML program defines a function that will transform a planar set of points according to one of the matrices in the dihedral group. The program reads the pinwheel coordinates into a data matrix and transforms the pinwheel according to each element of D4. The results are plotted by using PROC SGPANEL:

proc iml;
/* actions of the D4 dihedral group */
start D4Action(v, act);
   /* the subroup of rotations by 90 degrees */
   if      act=0 then M = { 1  0,  0  1};
   else if act=1 then M = { 0 -1,  1  0};
   else if act=2 then M = {-1  0,  0 -1};
   else if act=3 then M = { 0  1, -1  0};
   /* the subgroup of reflections across horz, vert, or diagonals */
   else if act=4 then M = {-1  0,  0  1};
   else if act=6 then M = { 0 -1, -1  0};
   else if act=7 then M = { 1  0,  0 -1};
   else if act=5 then M = { 0  1,  1  0};
   return( v*M` );  /* = (M*z`)` */
finish;
/* read in the pinwheel shape */
use Shape; read all var {x y} into P; 
           read all var "ID"; close;
/* write out the transformation of the pinwheel under the D4 actions */
OpNames = {"Identity" "R1" "R2" "R3" "S0" "S1" "S2" "S3"};
Name = OpNames[1];
Q = {. . .};
create Panel from Name Q[c={'Name' 'ID' 'x' 'y'}];
do i = 0 to 7;
   R = D4Action(P, i);
   Q = ID || R;
   Name = j(nrow(Q), 1, OpNames[i+1]);
   append from Name Q;
end;
close;
QUIT;
 
ods graphics / width=720px height=480px;
/* for convenience, define macros for the COLAXIS and ROWAXIS options */
%macro colOpts; colaxis offsetmin=0 offsetmax=0 display=(nolabel noticks novalues); %mend;
%macro rowOpts;  rowaxis offsetmin=0 offsetmax=0 display=(nolabel noticks novalues); %mend;
 
title "Dihedral Group (D4) Actions on Pinwheel Figure";
title2 "Colors=(Teal, Orange, Blue, Salmon)";
proc sgpanel data=Panel noautolegend;
   styleattrs wallcolor=&gray datacolors=(&teal &orange &blue &salmon);
   panelby Name / columns=4 onepanel;
   polygon x=x y=y ID=ID / group=ID fill;
   %colOpts; %rowOpts;
run;

The labels for each cell indicate the element of D4 that transforms the original pinwheel. The first row represents rotations by 0, 90, 180, and 270 degrees. The second row represents a vertical reflection (S0), a reflection across a diagonal line (S1), a reflection across a different diagonal (S2), and a horizontal reflection (S3). As mentioned earlier, I assumed that these eight images would be sufficient to reproduce the cells in Odita's artwork, but I was wrong. Let's look at the 16 cells in Odita's image and label each cell according to the elements of the D4 dihedral group:

I was surprised to discover that only half of the cells are transformations of the upper-left pinwheel. Furthermore, the four cells in the upper-left corner are arranged in the same order as the four cells in the lower-right corner. Out of the eight possible transformations of the pinwheel under the D4 group, Odita only uses four. The remaining cells either have the colors in a different order or have the "vane" of the pinwheel pointing in a different direction.

A second generator

Mathematically, Odita's pattern has a second pinwheel "generator." You could choose any of the unlabeled cells to be the generator of the remaining cells, but I will choose the cell on the second row and third column, which is displayed to the left.

For this pinwheel, the vanes are pointing in the same direction as the first pinwheel, but the colors (in counter-clockwise order) are teal, salmon, blue, and orange. Notice that the salmon and orange colors have switched their order.

In terms of the SAS data set, we can create this new colored pinwheel by switching the ID variable for the observations that have ID=1 and ID=3. You can then create a panel of images for this new pinwheel under the dihedral transformations, as follows:

/* switch the colors of the 2nd and 4th vanes */
data Panel2;
set Panel;
if      ID=1 then ID=3;
else if ID=3 then ID=1;
run;
 
title "Dihedral Group (D4) Actions on Pinwheel Figure";
title2 "Colors=(Teal, Salmon, Blue, Orange)";
proc sgpanel data=Panel2 noautolegend;
   styleattrs wallcolor=&gray datacolors=(&teal &orange &blue &salmon);
   panelby Name / columns=4 onepanel;
   polygon x=x y=y ID=ID / group=ID fill;
   %colOpts; %rowOpts;
run;

Let's use a prime to distinguish the elements of this second panel. That is, the first row represents the elements R`0, R`1, R`2, and R`3. The second row represents the elements S`0, S`1, S`2, and S`3. We can use these labels to identify the remaining cells in Odita's work, as follows:

The remaining eight cells display six of the possible patterns for the dihedral actions on the second generator. The S`0 and S`2 actions are repeated. The S`1 and S`3 actions do not appear.

Using SAS to generate mathematical art

Now that we have identified all the cells in Odita's painting, we could actually use PROC SGPANEL in SAS to create a mathematical reproduction of his work. Or, we can use SAS to create new Odita-inspired images.

I will do the latter. The following statements stack the panels of the dihedral group for the two fundamental pinwheels:

data Pinwheels;
length Name $20;
set Panel Panel2;
Group = floor((_N_-1)/20);
run;
 
ods graphics / width=640px height=640px;
title "Dihedral Shadow";
proc sgpanel data=Pinwheels noautolegend;
   styleattrs wallcolor=&gray datacolors=(&teal &orange &blue &salmon);
   panelby Group / columns=4 onepanel noheader;
   polygon x=x y=y ID=ID / group=ID fill;
   %colOpts; %rowOpts;
run;

In homage to Odita, I name this mathematical image "Dihedral's Shadow." The first two rows show the actions (rotations and reflections) of the D4 dihedral group on a pinwheel shape whose colors are in the sequential order "teal, orange, blue, and salmon." The last two rows show the dihedral actions on a pinwheel shape whose colors are in the order "teal, salmon, blue, and orange."

Notice that this mathematical image does not share some of the aesthetic properties of Odita's work. In Odita's image, the colors of the polygons are chosen so that the four neighbors (up-down and left-right) of a polygon are different colors from the polygon itself. My creation does not have that property. For example, if you think of the image as an 8 x 8 grid and traverse the first row, you see two teal polygons followed by two salmon polygons followed by two blue polygons. In Odita's work, no row has two adjacent polygons that have the same color.

I don't like the three orange polygons in the middle of my image. Odita arranged his pinwheel images to avoid neighboring polygons that have the same color. I suppose I could try arranging the pinwheels in a different order.

Summary

There is much more that could be said about the mathematical structures in Odita's work, "Phantom’s Shadow, 2018." In this article, I show how you can analyze components of the painting as reflections and rotations of two four-color pinwheel-shaped figures. Mathematically, these images can be generated by using orthogonal matrices. The images can be assembled into a lattice by using PROC SGPANEL in SAS.

This article shows how to create a purely mathematical image: the dihedral actions on two pinwheel shapes. It arranges the images without consideration of how a pinwheel looks adjacent to its neighbors. In a triumph of art over mathematics, I think Odita's art is more pleasing than my imitation.

Would you like to use SAS to create mathematical art? You can download the SAS program that I used to create the images in this article. I challenge my readers to use the Pinwheels data set to create their own mathematical works. You can use one of the following ideas, or explore your own ideas:

  • The design of experiments uses arrays to describe factors that are varied in the experiment. Use these pinwheel images to create a visual representation of an experimental design. The factors are rotation (4 levels), reflection (4 levels), and color-sequence (2 levels).
  • Rearrange the cells in "Dihedral's Shadow" to reduce the number of adjacent polygons that have the same color. You can do this manually, but extra points if you use math or SAS/OR software to minimize the number of adjacent polygons that have the same color!
  • Use a random number generation to create a random arrangement of the pinwheels.

I have created a thread on the SAS Support Communities that you can use to post your own creations. Have fun!

The post Use SAS to create mathematical art appeared first on The DO Loop.

8月 022021
 

Art evokes an emotional response in the viewer, but sometimes art also evokes a cerebral response. When I see patterns and symmetries in art, I think about a related mathematical object or process. Recently, a Twitter user tweeted about a painting called "Phantom’s Shadow, 2018" by the Nigerian-born artist, Odili Donald Odita. A modified version of the artwork is shown to the right.

The artwork is beautiful, but it also contains a lot of math. The image shows 64 rotations, reflections, and translations of a polygon in four colors. As I will soon explain, the image can also be viewed as 4 x 4 grid where each cell contains a four-bladed pinwheel. The grid displays rotations and reflections of a pinwheel shape.

When I saw this artwork, it inspired me to look closely at its mathematical structure and to create my own mathematical version of the artwork in SAS. This article shows how to use rotations to create a pinwheel from a polygon. A second article discusses how to use rotations and reflections to create a mathematical interpretation of Odita's painting.

Create a polygon

Look closely at the upper left corner of "Phantom's Shadow." You will see the following pinwheel-shaped figure:

If the center of the pinwheel is the origin, then this pinwheel is based on a sequence of 90-degree rotations of the teal-colored polygon about the origin. Each rotation is a different color (teal, orange, blue, and salmon) on a gray background.

You can assign coordinates to the vertices of the teal polygon. If three vertices are on the corners of the unit square, the fourth vertex appears to be at (2/3, 1/3). You can put the vertices into a SAS data set and graph the polygon by using the POLYGON statement in PROC SGPLOT:

%let alpha = 0.333;          /* Odita's work uses a vertex as (1-alpha, alpha) */
data Poly1;
ID = 1;
input x y @@;
if x=. then do;              /* just for fun, support other locations of vertex */
   x = 1 - α y = α
end;
datalines;
0 0   . .   1 1   0 1   0 0
;
 
%let blue   = CX0f5098;
%let orange = CXeba411;
%let teal   = CX288c95;
%let salmon = CXd5856e;
%let gray   = CX929386;
 
ods graphics / width=480px height=480px;
title "Base Polygon";
proc sgplot data=Poly1 aspect=1;
   styleattrs wallcolor=&gray;
   polygon x=x y=y ID=ID / fill fillattrs=(color=&teal) outline;
   xaxis offsetmin=0.005 offsetmax=0;
   yaxis offsetmin=0 offsetmax=0.005;
run;

The complete artwork repeats this polygon 64 times in a grid by using different colors and different orientations. But the colors and orientations are not random—although that would also look cool! Instead, there is an additional structure. The polygon is part of a pinwheel, and the pinwheel shape is rotated and reflected 16 times in a 4 x 4 grid. The next section shows how to create the pinwheel.

Create a pinwheel

It is possible to perform 2-D rotations and reflections in the DATA step, but the rotations and reflections of a planar figure are most simply expressed in terms of 2 x 2 orthogonal matrices.

To form the pinwheel from the basic polygon, you can use a group of four matrices. All rotations are in the counter-clockwise direction about the origin:

  • R0 is the identity matrix
  • R1 is the matrix that rotates a point 90 degrees (about the origin)
  • R2 is the matrix that rotates a point 180 degrees
  • R3 is the matrix that rotates a point 270 degrees

You do not need advanced mathematics to understand rotations by 90 degrees. However, these four rotation matrices are an example of an algebraic structure called a finite group. These matrices form the cyclic group of order 4 (C4), where the group operation is matrix multiplication.

Regardless of their name, you can use these matrices to transform the polygon into a pinwheel. The following SAS/IML program reads in the base polygon and transforms it according to each of the four matrices. You can then plot the four images, each in a different color.

/* rotate polygon about the origin to form a pinwheel */
proc iml;
/* actions of the C4 cyclic group: rotations by 90 degrees */
start C4Action(v, act);
   if      act=0 then M = { 1  0,  0  1};
   else if act=1 then M = { 0 -1,  1  0};
   else if act=2 then M = {-1  0,  0 -1};
   else if act=3 then M = { 0  1, -1  0};
   return( v*M` );  /* = (M*z`)` */
finish;
 
use Poly1;  read all var {x y} into P;  close;   /* read the polygon */
 
/* write the pinwheel */
Q = {. . .};
create Shape from Q[c={'ID' 'x' 'y'}];
do i = 0 to 3;
   R = C4Action(P, i);          /* rotated image of polygon     */
   Q = j(nrow(R), 1, i) || R;   /* save the image to a data set */
   append from Q;
end;
close;
QUIT;
 
title "Pinwheel Figure";
title2 "Colors=(Teal, Orange, Blue, Salmon)";
proc sgplot data=Shape aspect=1;
   styleattrs wallcolor=&gray datacolors=(&teal &orange &blue &salmon);
   polygon x=x y=y ID=ID / group=ID fill;
   xaxis offsetmin=0 offsetmax=0;
   yaxis offsetmin=0 offsetmax=0;
run;

The resulting image is the four-blade pinwheel, which is the basis of Odita's artwork.

A grid of pinwheels

The rest of "Phantom’s Shadow, 2018" can be described in terms of rotations and reflections of the pinwheel shape...with a twist! The next article introduces the dihedral group of order 4 (D4). This is a finite group of order eight and is the symmetry group of the square. It turns out that you can use the elements of D4 to generate the 16 pinwheels in Odita's artwork. You can then use PROC SGPANEL in SAS to generate a graph that looks like "Phantom's Shadow," or you can create your own Odita-inspired mathematical artwork! Stay tuned!

The post The art of rotations and reflections appeared first on The DO Loop.

11月 302020
 

Most games of skill are transitive. If Player A wins against Player B and Player B wins against Player C, then you expect Player A to win against Player C, should they play. Because of this, you can rank the players: A > B > C

Interestingly, not all games are transitive. The game "Rock, Paper, Scissors," is an intransitive game. In this game, Rock beats Scissors, Scissors beats Paper, and Paper beats Rock. You cannot rank the game shapes. You cannot declare that one shape wins more often than another.

In "Rock, Paper, Scissors," both players display their shapes simultaneously. If not, the second player could win 100% of the time by knowing the shape that is shown by the first player.

Intransitive dice

Recently, I saw an intransitive game that involves special six-sided dice. In this game, there are four dice: A, B, C, and D. Two players each choose a die and roll it. The higher number wins. If the players choose their die sequentially, the second player is at a huge advantage: No matter which die the first player chooses, the second player can choose a die such that his probability of winning the game is 2/3.

Clearly, the numbers on the faces of the dice are important since this advantage does not exist for the usual six-sided dice. One set of intransitive dice have faces as follows:
A = {0, 0, 4, 4, 4, 4}
B = {3, 3, 3, 3, 3, 3}
C = {2, 2, 2, 2, 6, 6}
D = {1, 1, 1, 5, 5, 5}

For these dice:

  • A wins against B with probability 2/3.
  • B wins against C with probability 2/3.
  • C wins against D with probability 2/3.
  • D wins against A with probability 2/3.

The probability of winning for pairs that involve die B are easy to compute because that die has a 3 on every face. By inspection, die A will win against B whenever A shows a 4, which occurs 2/3 of the time. Similarly, die B wins against die C whenever C shows a 2, which is 2/3 of the time.

The remaining probabilities can be computed by using a tree diagram for conditional probability. A tree diagram for the dice C and D are shown below. If C shows a 6 (which happens 1/3 of the time), then C wins. If C shows a 2 (which happens 2/3 of the time), then the game depends on the roll of D. If D rolls a 1 (which happens 1/2 of the time), then C wins. But if D rolls a 5 (which happens 1/2 of the time), then C loses. Since the rolls are independent, you multiply the individual probabilities to obtain the total probability of each pair of events. For C versus D, C wins with probability 2/3.

You can construct a similar diagram for the D versus A game. It is also interesting to compute the probability of A versus C and B versus D for these dice.

Simulate many games

Because conditional probability can be tricky, I always like to check my computations by running a simulation. The following SAS/IML program uses the SAMPLE function to simulate 10,000 rolls of each die. The estimated probability for winning is the number of times that one die showed the greater number, divided by 10,000, which is the total number of simulated games. For each of the four games considered, the estimated probability is approximately 2/3.

proc iml;
/* define the dice */
A = {0, 0, 4, 4, 4, 4};
B = {3, 3, 3, 3, 3, 3};
C = {2, 2, 2, 2, 6, 6};
D = {1, 1, 1, 5, 5, 5}; 
 
call randseed(1234);
N = 10000;
/* simulate N rolls for each die */
sA = sample(A, N);
sB = sample(B, N);
sC = sample(C, N);
sD = sample(D, N);
 
/* estimate the probability of 'A vs B', 'B vs C', etc */
AvsB = sum(sA > sB) / N;
BvsC = sum(sB > sC) / N;
CvsD = sum(sC > sD) / N;
DvsA = sum(sD > sA) / N;
print AvsB BvsC CvsD DvsA;

Further reading

There are other examples of intransitive dice. The ones in this article are called "Efron's Dice" because they were discovered by Brad Efron in 1970. Yes, the same Brad Efron who also invented bootstrap resampling methods! The following articles provide more examples and more explanation.

  • "Nontransitive dice", Wikipedia article with lots of examples and a section devoted to Efron's dice, as used in this article.
  • Grimes, J., (2010) "Curious Dice," Plus Magazine. A wonderful exposition.

The post Intransitive dice appeared first on The DO Loop.

8月 262020
 

In the paper "Tips and Techniques for Using the Random-Number Generators in SAS" (Sarle and Wicklin, 2018), I discussed an example that uses the new STREAMREWIND subroutine in Base SAS 9.4M5. As its name implies, the STREAMREWIND subroutine rewinds a random number stream, essentially resetting the stream to the beginning. I struggled to create a compelling example for the STREAMREWIND routine because using the subroutine "results in dependent streams of numbers" and because "it is usually not necessary in simulation studies" (p. 12). Regarding an application, I asserted that the subroutine "is convenient for testing."

But recently I was thinking about two-factor authentication and realized that I could use the STREAMREWIND subroutine to emulate generating a random token that changes every 30 seconds. I think it is a cool example, and it gives me the opportunity to revisit some of the newer features of random-number generation in SAS, including new generators and random number keys.

A brief overview of two-factor authentication

I am not an expert on two-factor authentication (TFA), but I use it to access my work computer, my bank accounts, and other sensitive accounts. The main idea behind TFA is that before you can access a secure account, you must authenticate yourself in two ways:

  • Provide a valid username and password.
  • Provide information that depends on a physical device that you own and that you have previously registered.

Most people use a smartphone as the physical device, but it can also be a PC or laptop. If you do an internet search for "two factor authentication tokens," you can find many images like the one on the right. This is the display from a software program that runs on a PC, laptop, or phone. The "Credential ID" field is a long string that is unique to each device. (For simplicity, I've replaced the long string with "12345.") The "Security Code" field displays a pseudorandom number that changes every 30 seconds. The Security Code depends on the device and on the time of day (within a 30-second interval). In the image, you can see a small clock and the number 28, which indicates that the Security Code will be valid for another 28 seconds before a new number is generated.

After you provide a valid username and password, the account challenges you to type in the current Security Code for your registered device. When you submit the Security Code, the remote server checks whether the code is valid for your device and for the current time of day. If so, you can access your account.

Two-factor random number streams

I love the fact that the Security Code is pseudorandom and yet verifiable. And it occurred to me that I can use the main idea of TFA to demonstrate some of the newer features in the SAS random number generators (RNGs).

Long-time SAS programmers know that each stream is determined by a random number seed. But a newer feature is that you can also set a "key" for a random number stream. For several of the new RNGs, streams that have the same seed but different keys are independent. You can use this fact to emulate the TFA app:

  • The Credential ID (which is unique to each device) is the "seed" for an RNG.
  • The time of day is the "key" for an RNG. Because the Security Code must be valid for 30 seconds, round the time to the nearest 30-second boundary.
  • Usually each call to the RAND function advances the state of the RNG so that the next call to RAND produces a new pseudorandom number. For this application, we want to get the same number for any call within a 30-second period. One way to do this is to reset the random number stream before each call so that RAND always returns the FIRST number in the stream for the (seed, time) combination.

Using a key to change a random-number stream

Before worrying about using the time of day as the key value, let's look at a simpler program that returns the first pseudorandom number from independent streams that have the same seed but different key values. I will use PROC FCMP to write a function that can be called from the SAS DATA step. Within the DATA step, I set the seed value and use the "Threefry 2" (TF2) RNG. I then call the Rnd6Int function for six different key values.

proc fcmp outlib=work.TFAFunc.Access;
   /* this function sets the key of a random-numbers stream and 
      returns the first 6-digit pseudorandom number in that stream */
   function Rnd6Int(Key);
      call stream(Key);               /* set the Key for the stream */
      call streamrewind(Key);         /* rewind stream with this Key */
      x = rand("Integer", 0, 999999); /* first 6-digit random number in stream */
      return( x );
   endsub;
quit;
 
options cmplib=(work.TFAFunc);       /* DATA step looks here for unresolved functions */
data Test;
DeviceID = 12345;                    /* ID for some device */
call streaminit('TF2', DeviceID);    /* set RNG and seed (once per data step) */
do Key = 1 to 6;
   SecCode = Rnd6Int(Key);           /* get random number from seed and key values */
   /* Call the function again. Should produce the same value b/c of STREAMREWIND */
   SecCodeDup = Rnd6Int(Key);  
   output;
end;
keep DeviceID Key SecCode:;
format SecCode SecCodeDup Z6.;
run;
 
proc print data=Test noobs; run;

Each key generates a different pseudorandom six-digit integer. Notice that the program calls the Rnd6Int function twice for each seed value. The function returns the same number each time because the random number stream for the (seed, key) combination gets reset by the STREAMREWIND call during each call. Without the STREAMREWIND call, the function would return a different value for each call.

Using a time value as a key

With a slight modification, the program in the previous section can be made to emulate the program/app that generates a new TFA token every 30 seconds. However, so that we don't have to wait so long, the following program sets the time interval (the DT macro) to 10 seconds instead of 30. Instead of talking about a 30-second interval or a 10-second interval, I will use the term "DT-second interval," where DT can be any time interval.

The program below gets the "key" by looking at the current datetime value and rounding it to the nearest DT-second interval. This value (the RefTime variable) is sent to the Rnd6Int function to generate a pseudorandom Security Code. To demonstrate that the program generates a new Security Code every DT seconds, I call the Rnd6Int function 10 times, waiting 3 seconds between each call. The results are printed below:

%let DT = 10;                  /* change the Security Code every DT seconds */
 
/* The following DATA step takes 30 seconds to run because it
   performs 10 iterations and waits 3 secs between iterations */
data TFA_Device;
keep DeviceID Time SecCode;
DeviceID = 12345;
call streaminit('TF2', DeviceID);   /* set the RNG and seed */
do i = 1 to 10;
   t = datetime();                  /* get the current time */
   /* round to the nearest DT seconds and save the "reference time" */
   RefTime = round(t, &DT); 
   SecCode = Rnd6Int(RefTime);      /* get a random Security Code */
   Time = timepart(t);              /* output only the time */
   call sleep(3, 1);                /* delay 3 seconds; unit=1 sec */
   output;
end;
format Time TIME10. SecCode Z6.;
run;
 
proc print data=TFA_Device noobs; 
   var DeviceId Time SecCode;
run;

The output shows that the program generated three different Security Codes. Each code is constant for a DT-second period (here, DT=10) and then changes to a new value. For example, when the seconds are in the interval [05, 15), the Security Code has the same value. The Security Code is also constant when the seconds are in the interval [15, 25) and so forth. A program like this emulates the behavior of an app that generates a new pseudorandom Security Code every DT seconds.

Different seeds for different devices

For TFA, every device has a unique Device ID. Because the Device ID is used to set the random number seed, the pseudorandom numbers that are generated on one device will be different than the numbers generated on another device. The following program uses the Device ID as the seed value for the RNG and the time of day for the key value. I wrapped a macro around the program and called it for three hypothetical values of the Device ID.

%macro GenerateCode(ID, DT);
data GenCode;
   keep DeviceID Time SecCode;
   format DeviceID 10. Time TIME10. SecCode Z6.;
   DeviceID = &ID;
   call streaminit('TF2', DeviceID); /* set the seed from the device */
   t = datetime();                   /* look at the current time */
   /* round to the nearest DT seconds and save the "reference time" */
   RefTime = round(t, &DT);          /* round to nearest DT seconds */
   SecCode = Rnd6Int(RefTime);       /* get a random Security Code */
   Time = timepart(t);               /* output only the time */
run;
 
proc print data=GenCode noobs; run;
%mend;
 
/* each device has a unique ID */
%GenerateCode(12345, 30);
%GenerateCode(24680, 30);
%GenerateCode(97531, 30);

As expected, the program produces different Security Codes for different Device IDs, even though the time (key) value is the same.

Summary

In summary, you can use features of the SAS random number generators in SAS 9.4M5 to emulate the behavior of a TFA token generator. The SAS program in this article uses the Device ID as the "seed" and the time of day as a "key" to select an independent stream. (Round the time into a certain time interval.) For this application, you don't want the RAND function to advance the state of the RNG, so you can use the STREAMREWIND call to rewind the stream before each call. In this way, you can generate a pseudorandom Security Code that depends on the device and is valid for a certain length of time.

The post Rewinding random number streams: An application appeared first on The DO Loop.

8月 122019
 
What is this math good for, anyway?
     –Every student, everywhere

I am a professional applied mathematician, yet many of the mathematical and statistical techniques that I use every day are not from advanced university courses but are based on simple ideas taught in high school or even in grade school. I've written many blog posts in which the solution to an interesting problem requires little more than high-school math. Even when the solution requires advanced techniques, high-school math often provides the basis for solving the problem.

In celebration of the upcoming school year, here are 12 articles that show connections between advanced topics and mathematical ideas that are introduced in high-school or earlier. If you have a school-age child, read some of the articles to prepare yourself for the inevitable mid-year whine, "But, Mom/Dad, when will I ever use this stuff?!"

Grade-school math

Obviously, most adults use basic arithmetic, fractions, decimals, and percents, but here are some less obvious uses of grade-school math:

High-school math

Algebra, linear transformations, geometry, and trigonometry are the main topics in high school mathematics. These topics are the bread-and-butter of applied mathematics:

  • Linear transformations: Anytime you create a graph on the computer screen, a linear transformation transforms the data from physical units (weight, cost, time) into pixel values. Although modern graphical software performs that linear transformation for you, a situation in which you have to manually apply a linear transformation is when you want to display data in two units (pounds and kilograms, dollars and Euros, etc). You can use a simple linear transformation to align the tick labels on the two axes.
  • Intersections: In high school, you learn to compute the intersection between two lines. You can extend the problem to find the intersection of two line segments. I needed that result to solve a problem in probability theory.
  • Solve a system of equations: In Algebra II, you learn to solve a system of equations. Solving linear systems is the foundation of linear algebra. Solving nonlinear systems is among the most important skills in applied mathematics.
  • Find the roots of a nonlinear equation: Numerically finding the roots of a function is taught in pre-calculus. It is the basis for all "inverse problems" in which you want to find inputs to a function that produce a specified output. In statistics, a common "inverse problem" is finding the quantile of a cumulative probability distribution.
  • Binomial coefficients: In algebra, many teachers use the mnemonic FOIL (First, Outer, Inner, Last) to teach students how to compute the quadratic expansion of (a + b)2. Later, students learn the binomial expansion of an arbitrary power, (a + b)n. The coefficients in this expansion are called the binomial coefficients and appear in Pascal's triangle as well as in many discrete probability distributions such as the negative binomial and hypergeometric distributions.
  • Pythagorean triples: In trigonometry, a huge number of homework problems involve using right triangles with side lengths that are proportional to (3, 4, 5) and (5, 12, 13). These are two examples of Pythagorean triples: right triangles whose side lengths are all integers. It turns out that you can use linear transformations to generate all primitive triples from the single triple (3, 4, 5). A high school student can understand this process, although the process is most naturally expressed in terms of matrix multiplication, which is not always taught in high school.

High-school statistics

Many high schools offer a unit on probability and statistics, and some students take AP Statistics.

Einstein famously said, "everything should be made as simple as possible, but not simpler." It is surprising to me how often an advanced technique can be simplified and explained by using elementary math. I don't claim that "everything I needed to know about math I learned in kindergarten," but I often return to elementary techniques when I describe how to solve non-elementary problems.

What about you? What are some elementary math or statistics concepts that you use regularly in your professional life? Are there fundamental topics that you learned in high school that are deeper and more useful than you realized at the time? Leave a comment.

The post The math you learned in school: Yes, it’s useful! appeared first on The DO Loop.

8月 122019
 
What is this math good for, anyway?
     –Every student, everywhere

I am a professional applied mathematician, yet many of the mathematical and statistical techniques that I use every day are not from advanced university courses but are based on simple ideas taught in high school or even in grade school. I've written many blog posts in which the solution to an interesting problem requires little more than high-school math. Even when the solution requires advanced techniques, high-school math often provides the basis for solving the problem.

In celebration of the upcoming school year, here are 12 articles that show connections between advanced topics and mathematical ideas that are introduced in high-school or earlier. If you have a school-age child, read some of the articles to prepare yourself for the inevitable mid-year whine, "But, Mom/Dad, when will I ever use this stuff?!"

Grade-school math

Obviously, most adults use basic arithmetic, fractions, decimals, and percents, but here are some less obvious uses of grade-school math:

High-school math

Algebra, linear transformations, geometry, and trigonometry are the main topics in high school mathematics. These topics are the bread-and-butter of applied mathematics:

  • Linear transformations: Anytime you create a graph on the computer screen, a linear transformation transforms the data from physical units (weight, cost, time) into pixel values. Although modern graphical software performs that linear transformation for you, a situation in which you have to manually apply a linear transformation is when you want to display data in two units (pounds and kilograms, dollars and Euros, etc). You can use a simple linear transformation to align the tick labels on the two axes.
  • Intersections: In high school, you learn to compute the intersection between two lines. You can extend the problem to find the intersection of two line segments. I needed that result to solve a problem in probability theory.
  • Solve a system of equations: In Algebra II, you learn to solve a system of equations. Solving linear systems is the foundation of linear algebra. Solving nonlinear systems is among the most important skills in applied mathematics.
  • Find the roots of a nonlinear equation: Numerically finding the roots of a function is taught in pre-calculus. It is the basis for all "inverse problems" in which you want to find inputs to a function that produce a specified output. In statistics, a common "inverse problem" is finding the quantile of a cumulative probability distribution.
  • Binomial coefficients: In algebra, many teachers use the mnemonic FOIL (First, Outer, Inner, Last) to teach students how to compute the quadratic expansion of (a + b)2. Later, students learn the binomial expansion of an arbitrary power, (a + b)n. The coefficients in this expansion are called the binomial coefficients and appear in Pascal's triangle as well as in many discrete probability distributions such as the negative binomial and hypergeometric distributions.
  • Pythagorean triples: In trigonometry, a huge number of homework problems involve using right triangles with side lengths that are proportional to (3, 4, 5) and (5, 12, 13). These are two examples of Pythagorean triples: right triangles whose side lengths are all integers. It turns out that you can use linear transformations to generate all primitive triples from the single triple (3, 4, 5). A high school student can understand this process, although the process is most naturally expressed in terms of matrix multiplication, which is not always taught in high school.

High-school statistics

Many high schools offer a unit on probability and statistics, and some students take AP Statistics.

Einstein famously said, "everything should be made as simple as possible, but not simpler." It is surprising to me how often an advanced technique can be simplified and explained by using elementary math. I don't claim that "everything I needed to know about math I learned in kindergarten," but I often return to elementary techniques when I describe how to solve non-elementary problems.

What about you? What are some elementary math or statistics concepts that you use regularly in your professional life? Are there fundamental topics that you learned in high school that are deeper and more useful than you realized at the time? Leave a comment.

The post The math you learned in school: Yes, it’s useful! appeared first on The DO Loop.