SAS博客列表

SAS programmers love to make special graphs for Valentine's Day. In fact, there is a long history of heart-shaped graphs and love-inspired programs written in SAS! Last year, I added to the collection by showing how a ball bounces on a heart-shaped billiards table. This year, I create a similar image, but one that is stochastic instead of deterministic. This year's image is a random walk in a heart-shaped region.

In the figure to the right, the red lines trace the path of a random walk that starts in the middle of the heart. For each step, the new location is determined by a random bivariate normal vector. That is, from the position (x,y), the next location is (x+dx, y+dy) where (dx, dy) is a random bivariate normal vector. If the next location would be outside of the heart-shaped region, then shorten the step and land on the boundary. The figure shows the path of a constrained random walk for 2000 steps.

Perhaps this image is a visual metaphor for randomly exploring love? Some readers will be content to view this image and reflect on the love in their lives. Others will want to learn more about how the image is created in SAS. If you are in the second group, read on!

An unconstrained random walk

Before I show how to constrain the random walk to stay inside the heart, let's create and visualize an unconstrained random walk. There are many kinds of random walks, but this one takes a step according to a random bivariate normal vector. Starting from a point (x,y), you generate a vector (dx,dy) ~ MVN(0, I(2)) and move to the new point (x+dx, y+dy). You repeat this process many time, then draw the path that results. The following SAS DATA step generates the random walk by using 0.2 as the standard deviation of each variate:

/* random 2-D walk where each step is a vector v, which is bivariate normal
   v ~ MVN(0, I(2)) */
data RandomWalk;
N = 2000;
call streaminit(12345);
x = 0; y = 0;
do i = 1 to N;
   dx = rand("Normal", 0, 0.2);
   dy = rand("Normal", 0, 0.2);
   x = x + dx;
   y = y + dy;
   output;
end;
keep x y;
run;
ods graphics / reset width=480px height=480px;
title "Random Walk with Bivariate Normal Steps";
proc sgplot data=RandomWalk;
   series x=x y=y / lineattrs=(color=lightred);
   xaxis grid;
   yaxis grid;
run;

This graph displays a random walk with 2,000 steps, starting from the origin. Most of the time, the two consecutive points are close to each other. However, sometimes the next point is far away from the previous point.

A constrained random walk

Now imagine a variation on the previous random walk. The initial point is inside a compact closed region. Instead of always moving from an initial point (x,y) to a final point (x+dx, y+dy), you compute whether the final point is inside the region. If so, move to the point. If not, truncate the step so that you land on the boundary of the region. Repeat this many times to obtain a constrained random walk that is always inside the region.

For Valentine's Day, I will choose the region to be heart-shaped. It is easiest to perform these computations in the SAS IML language. The following program defines a few helper functions and generates 2,000 points for the constrained random walk. The points are written to a SAS data set and overlaid on a graph of the heart-shaped region:

/* Instead of an unconstrained walk, at each step check whether the 
   new location is outside of a region. If so, truncate the step to 
   land on the boundary of the region. */
proc iml;
/* define the heart-shaped region */
start H(xy);
  x = xy[,1]; y = xy[,2];
  return ( (x**2 + y**2 - 1)**3 - x**2 * y**3 );
finish;
 
/* return 1 if (x,y) is inside the region */
start InRegion(xy, tol=1e-14);
  return (H(xy) <= tol);
finish;
 
/* given a point, x, and a vector, v, this function returns the function
   f(t) = H(x + t*v), t in [0,1]. If f has a root, then the line segment from 
   x to x+v intersects the boundary of the reqion. 
   This function is used by FROOT to find points on the boundary of the region. */
start OnBoundary(t) global(G_x, G_v);
  return ( H( G_x + t*G_v ) );
finish;
 
/* Start at x. Try to step to x+v.
   If x+v is in the region, take the step.
   Otherwise, see if we can take part of the step to land on the boundary. */
start StepInRegion(x, v) global(G_x, G_v);
   if InRegion(x + v) then 
      return (x + v);
   if InRegion(x) then do;
      G_x = x; G_v = v;
      /* does the line from x to x+v cross the boundary? */
      t = froot("OnBoundary", {0 1});  /* find intersection of line and region */
      if t > 0 then 
         return (x + t*v);       /* step onto the boundary */
      else
         return (x);             /* cannot step in this direction */
   end;
   /* something is wrong: x is not in the region */
   return ( {. .} );
finish;
 
N = 2000;
call randseed(12345);
vel = randfun(N//2, "Normal", 0, 0.2);
 
x = j(1, 2, 0);
create walk from x [c={x y}];
do i = 1 to N;
   x = StepInRegion(x, vel[i,]);
   append from x;
end;
close;
QUIT;
 
title "Random Walk Inside a Heart-Shaped Region";
proc sgplot data=Walk;
   series x=x y=y / lineattrs=(color=lightred);
run;

The trajectory of the random walk is shown for 2,000 random steps inside a heart-shaped region. You can carry out this computation for any compact region, provided that you can represent the region as the level set of some continuous function H(x,y)=0 that divides the plane into an inside region (where H(x,y) < 0) and an outside region (where H(x,y) > 0).

To visualize the heart-shaped region, you can overlay the trajectory on the region, as done at the top of this article. The images below visualize how the random walk evolves and gradually fills the region. For details about how these images were created, see the SAS program that creates these images.

Summary

It is easy to simulate an unconstrained random walk. This article shows how to use SAS to generate a random walk in which the points are contained in a compact region. During the simulation, you can monitor whether the next step will occur outside the region. If so, you can decrease the step so that the point is on the boundary of the region. This is illustrated by using a heart-shaped region, but the same algorithm can constrain the random walk to any compact region of the plane.

The post A random walk inside a heart appeared first on The DO Loop.

One of the benefits of social media is the opportunity to learn new things. Recently, I saw a post on Twitter that intrigued me. The tweet said that the expected volume of a random tetrahedron in the unit cube (in 3-D) is E[Volume] = 0.0138427757.... This number seems surprisingly small! After all, the cube has unit volume, and the volume of the largest tetrahedron in the cube is 1/3. So, the result says that the average volume of a random tetrahedron is 24 times smaller than the largest tetrahedron!

The result is attributed to Zinani, A. (2003) "The Expected Volume of a Tetrahedron whose Vertices are Chosen at Random in the Interior of a Cube," which is behind a paywall. The abstract says that the analytical result is exactly
      E[Volume] = 3977 / 216000 - π² / 2160
and that this number is "in convincing agreement with a simulation."

This article shows how to run a Monte Carlo simulation in SAS to simulate random tetrahedra and compute their volumes. We use the simulation to confirm Zinani's result and to visualize the distribution of the volume.

The volume of a tetrahedron

In multivariable calculus, students are introduced to the scalar triple product of three vectors. Given three vectors u, v, and w, the scalar triple product is the quantity (u x v)⋅w, where (u x v) is the cross product of the vectors, and z⋅w = z^T*w is the scalar dot product. The scalar triple product is the (signed) volume of the parallelepiped formed by the vectors.

Recall that the volume of a tetrahedron with vertices (A,B,C,D) is 1/6 of the volume of the parallelepiped formed by the vectors between vertices. Consequently, the volume of a tetrahedron is
      V = 1/6 |(AB x AC)⋅AD|
where AB is the vector from A to B, AC is the vector from A to C, and AD is the vector from A to D.

Compute the volume of a tetrahedron in SAS

Since the volume formula includes vectors and vector operations, it is natural to use the SAS/IML language to compute the volume. The following functions implement the definition of the vector cross-product. To test the function, the program computes the volume of the largest tetrahedron that can be imbedded in the unit cube. One such tetrahedron is shown to the right. The tetrahedron has vertices A(0,0,0), B(1,1,0), C(0,1,1), and D(1,0,1). The volume of this tetrahedron is 1/3, as shown by the following program.

/* volume of tetrahedron in 3D */
proc iml;
/* vector cross product (3D vectors only)
   https://blogs.sas.com/content/iml/2013/05/20/matri-multiplication.html
*/
start CrossProd(u, v);
   i231 = {2 3 1};   i312 = {3 1 2};
   return( u[i231]#v[i312] - v[i231]#u[i312] );
finish;
 
start VolumeTetra(A,B,C,D);
   z = CrossProd(B-A, C-A);    /* z = AB x AC */
   w = colvec(D - A);          /* w = AD */
   volume = abs(z`*w) / 6;     /* V = 1/6 | z`*w | */
   return( volume );
finish;
 
/* test the formula: compute the volume of the largest tetrahedron inside the unit cube */
A = {0 0 0};
B = {1 1 0};
C = {0 1 1}; 
D = {1 0 1};
maxVol = VolumeTetra(A,B,C,D);
print maxVol;

A Monte Carlo simulation of random tetrahedron

It is not hard to write a SAS/IML program that generates four random vertices for a tetrahedron inside the unit cube in 3-D. To generate the four vertices, you can generate 12 uniform variates independently in the interval [0,1]. Assign the first three numbers to be the (x,y,z) Euclidean coordinates of the first vertex, the next three numbers to be the coordinates of the second vertex, and so forth. From the vertices, you can call the VolumeTetra function to produce the volume.

The following program generates 100,000 random tetrahedra and computes their volumes. The program then calls the MEAN function computes the average volume, which is the Monte Carlo estimate of the expected value.

call randseed(1234);
N = 1E5;
X = randfun(N//12, "Uniform");   /* Z x 12 matrix of uniform variates */
Vol = j(N, 1, .);                /* allocate array for volumes */
do i = 1 to N;
   A = X[i, 1:3];   B = X[i, 4:6];  /* form four vertices from the 12 random numbers */
   C = X[i, 7:9];   D = X[i, 10:12];
   Vol[i] = VolumeTetra(A,B,C,D);
end;
 
avgVol = mean(Vol);                                       /* Monte Carlo estimate of mean */
expectedValue = 3977/216000 - constant('pi')**2 / 2160;   /* Zinani's exact result for E[Volume] */
print avgVol expectedValue (expectedValue-avgVol)[L="Diff"];

The output shows that the Monte Carlo estimate is very close to the expected value. This confirms Zinani's result: the average volume of a random tetrahedron is smaller than you might think.

Visualize the distribution of volumes

The expected value for the distribution of tetrahedral volumes is small. Since we know that at least one tetrahedron has a relatively large value (1/3), we might conjecture that most of the tetrahedra have very small volumes and only a few have larger volumes. A way to test that conjecture is to visualize the distribution by drawing a histogram of the volumes for the 100,000 random tetrahedra. The following SAS/IML statements write the expected value to a macro variable and write the volumes of the random tetrahedra to a data set. PROC UNIVARIATE then computes descriptive statistics for the volumes and plots a histogram of the distribution of volumes:

/* E[Volume] = 3977 / 216000 - pi**2 / 2160 = 0.01384277574...  */
call symputx("EV",expectedValue);       /* write E[Volume] to a macro variable */
create Vol var "Vol";  append;  close;  /* write volumes to a data set */
QUIT;
 
title "Volume of a Random Tetrahedron";
proc univariate data=Vol;
   histogram Vol / endpoints href=&EV odstitle=title;
run;

The graph indicates that a random tetrahedron is likely to have a volume that is near 0. The vertical line is the expected value of the distribution. Notice that only a small proportion of the tetrahedra have volumes that exceed 0.1. In this collection of 100,000 random tetrahedra, the largest volume was 0.1777, which is about 50% of the volume of the largest possible tetrahedron. The distribution has a long but thin tail.

If you look carefully at the table of descriptive statistics (the "Moments" table), you will see that the mean and the standard deviation are nearly the same. This property is characteristic of an exponential distribution. An exponential distribution has a skewness of 2 and an excess kurtosis of 6. However, the sample skewness and sample kurtosis are biased statistics, and the estimates shown here (1.8 and 5.1, respectively) are typical for a large sample from an exponential distribution. This suggests that the tetrahedral volume is distributed as V ~ Exp(λ=E[Volume]). Let's use PROC UNIVARIATE again and overlay the density of the Exp(λ=0.0139) distribution. Furthermore, you can use a Q-Q plot to visualize how well the exponential model fits the simulated volumes.

proc univariate data=Vol;
   histogram Vol / exponential(theta=0 sigma=&EV) endpoints;
   qqplot Vol / exponential(theta=0 sigma=&EV);
run;

The fit is excellent! Not only does the density curve overlay perfectly on the histogram, but the quantiles of the data (not shown) are in excellent agreement with the expected quantiles of the exponential distribution. If you look at the Q-Q plot, you again see excellent agreement between the data and the exponential model. The only shortcoming is that the simulated volumes have fewer extreme values than predicted by the model.

The model enables you to compute the probability of certain events. For example, the probability that a random tetrahedron has volume greater than 0.1 is 1-cdf("Expo", 0.1, &EV), which is 0.00073.

Summary

A result by Zinani (2003) gives the expected volume of a tetrahedron whose vertices are chosen uniformly at random in the unit cube. Zinani's result is E[Volume] = 0.0138427757.... Although this number seems surprisingly small, the expected value is confirmed by using a Monte Carlo simulation in SAS. The simulation suggests that the volumes are exponentially distributed.

For an alternative mathematical proof of the result, see Philip, J. (2007) "The Expected Volume of a Random Tetrahedron in a Cube", which is free to view and is very well written.

The post The expected volume of a random tetrahedron in a cube appeared first on The DO Loop.