Wolfram Blog
Oleksandr Pavlyk

A “Trivial” Probability Problem

June 3, 2013 — Oleksandr Pavlyk, Kernel Technology

I am a junkie for a good math problem. Over the weekend, I encountered such a good problem on a favorite subject of mine–probability. It’s the last problem from the article “A Mathematical Trivium” by V. I. Arnol’d, Russian Mathematical Surveys 46(1), 1991, 271–278.

It’s short enough to reproduce in its entirety: “Find the mathematical expectation of the area of the projection of a cube with edge of length 1 onto a plane with an isotropically distributed random direction of projection.” In other words, what is the average area of a cube’s shadow over all possible orientations?

This blog post explores the use of Mathematica to understand and ultimately solve the problem. It recreates how I approached the problem.

Projection of a square on a line

Before tackling the case of the cuboid, I started with a unit square, randomly rotated around its center of mass, with the intent to find the average length of its projection on a horizontal axis.

For the sake of simplicity, I placed the center of the mass of the square at the origin.

rectangleWireStyle = Directive[Thick, Darker[Orange]]; rectangleWireFrame =    Line[{{-1/2, -1/2}, {1/2, -1/2}, {1/2, 1/2}, {-1/2,       1/2}, {-1/2, -1/2}}];

I found the left and right boundaries of the projection as the smallest and largest x coordinates of the vertices of the square, rotated around the origin through an angle of α degrees:

projlen[\[Theta]_] = Block[{xcoords, cubeVertices},   cubeVertices =     RotationTransform[\[Theta]] /@ Tuples[{{-1/2, 1/2}, {-1/2, 1/2}}];   xcoords = First /@ cubeVertices;   Max[xcoords] - Min[xcoords]]

Max[-(Cos[θ]/2) - Sin[θ]/2, Cos[θ]/2 - Sin[θ]/2, -(Cos[θ]/2) + Sin[θ]/2, 
  Cos[θ]/2 + Sin[θ]/2] - 
 Min[-(Cos[θ]/2) - Sin[θ]/2, Cos[θ]/2 - Sin[θ]/2, -(Cos[θ]/2) + Sin[θ]/2, 
  Cos[θ]/2 + Sin[θ]/2]

I then combined these functions within the Manipulate to be able to dynamically change the angle of rotation:

Manipulate to be able to dynamically change the angle of rotation

Here’s the plot of the length of the projection as the function of the rotation angle:

Plot of the length of the projection as the function of the rotation angle

Projection length of the rotated square

Assuming the rotation angle α is uniformly distributed in the interval 0 ≤ α < 360, I readily found the expected length of the projection, that is, the average length of the square’s shadow:

Expectation[projlen[\[Alpha] Degree], \[Alpha] \[Distributed] UniformDistribution[{0, 360}]]

4/π

Change the frame of reference of the rectangle

That was easy enough, but in order to take the victory to 3D, I needed to change the point of view. Instead of rotating the square, I rotated the line we project on.

Perspective of the wire frame:

I found this change of perspective illuminating, as it made me think that the length of the projection is the sum of the lengths of the projections of the two sides of the square visible from the plane.

The length of the projection of a segment of length &#8467 with unit normal vector n1 onto a line with unit normal vector n2 equals &#8467 Dot[n1,n2].

The projection of each side of the square only contributes if Dot[n1,n2] is positive (i.e. the side is visible); otherwise it is hidden behind other sides. The length of the shadow is thus the sum of the contributions of the east, west, top, and bottom sides of the square:

The sum of the contributions of the east, west, top, and bottom sides of the square

Thus the length of the projection is the sum of the absolute values of coordinate components of the normal vector nvec. I implemented this way of computing the length of the shadow in a function:

projlen2[θ_] := 
 Block[{nvec = {-Sin[θ], -Cos[θ]}}, Abs[nvec.{1, 0}] + Abs[nvec.{0, 1}] ]

And, of course, it agrees with the earlier way of computing the projection length:

Computing the projection length

Plotting the projection length

Naturally, the expectation is the same:

Expectation[projlen2[α Degree], α ⎡ UniformDistribution[{0, 360}]]

4/π

The case of the cube

Guided by the insight gained by considering the square, I adopted the reference frame of the cuboid, whose center of mass is situated at the origin. The cuboid casts its shadow on a plane whose orientation is parametrized by its perpendicular (i.e. normal) vector nvec.

I started by straightforwardly building projections of each face of the cuboid and drawing them together. First I defined this helper function to project a vector xvec onto a plane with normal nvec:

ProjectVector[xvec_, nvec_] := xvec - nvec Dot[nvec, xvec]/Dot[nvec, nvec]

The following function gives vertex coordinates of a face in the (i,j) plane with the other coordinate being z0. Here i and j range over {1,2,3}, which designates the x-, y-, and z- directions, and z0 ranges over {-1/2,1/2}:

Vertex coordinates of a face in the (i,j) plane with the other coordinate being Subscript[z, 0]

I then defined a function to project each face onto a plane with normal vector nvec and produce the corresponding 3D polygon:

ProjectionShadow[nvec_] :=
 
 Flatten[Table[
   Polygon[
    Table[ProjectVector[v, nvec], {v, face[ij, z0]}]], {z0, {-1/2, 
     1/2}}, {ij, {{1, 2}, {1, 3}, {2, 3}}}]]

The area of a polygon that happens to be a quadrangle projected onto a plane with normal nvec is computed as the sum of the areas of two triangles that the quadrangle is split into by a diagonal:

PolygonArea[Polygon[{v1_, v2_, v3_, v4_, v1_}], 
  nvec_] := (Abs[Cross[v1 - v2, v3 - v2].nvec] + 
    Abs[Cross[v3 - v4, v1 - v4].nvec])/2

With these utilities in place, I was ready to visualize the projection of the cuboid with edge of length one, whose center of mass is situated at the origin and whose sides are aligned with the coordinate axis. I parametrized the unit normal vector to the projection plane using spherical coordinates θ and φ: {sin(θ) sin(φ),sin(θ) cos(φ),cos(θ)}.

Visualize the projection of the cuboid

Since at any one time, only three of the cuboid’s faces are visible, and since the area contributed by the invisible faces is the same (just imagine a parallel plane on another side of the cuboid), I sum over all the six faces and divide the result by two. For a particular orientation nvec=={-1,2,3}√14 of the plane, the area of the projection therefore is computed as follows:

With[{nvec = {-1, 2, 3}/Sqrt[14]}, 
 1/2 Sum[PolygonArea[poly, nvec], {poly, ProjectionShadow[nvec]}]]

3 Sqrt[2/7]

Applying the insight learned from solving the 2D case, I checked if the area of the projection was again the sum of the absolute values of dot products of the normal vector nvec with axis vectors:

With[{nvec = {-1, 2, 3}/Sqrt[14]},  Abs[Dot[nvec, {1, 0, 0}]] +    Abs[Dot[nvec, {0, 1, 0}] + Abs[Dot[nvec, {0, 0, 1}]]]]

3 Sqrt[2/7]

Bingo! This makes sense, I thought. Each term corresponds to the area of one of at most three visible faces. Indeed, consider a patch of area S on a plane with unit normal vector n1. The area of the projection of the patch onto another plane with unit normal vector n2 equals S Abs[Dot[n1,n2]].

ShadowArea[θ_, ϕ_] := 
 Block[{nvec = {Sin[θ] Cos[ϕ], Sin[θ ] Sin[ϕ], Cos[θ]}}, Total[Abs[nvec]] ]

Here is the spherical plot of the area of the shadow cast by the cuboid onto a plane with unit normal vector {sin(θ) sin(φ),sin(θ) cos(φ),cos(θ)} as the function of Euler’s angles θ and φ:

SphericalPlot3D of the area of the shadow

Spherical plot of the area of the shadow

The minimum of the area function corresponds to the area of one face, which equals 1, and the maximum of √3 corresponds to the projection on the plane, whose normal vector aligns with the cuboid’s diagonal:

ShadowArea[0, Pi/2]

1

ShadowArea[ArcCos[1/Sqrt[3]], Pi/4]

Sqrt[3]

I was then almost ready to tackle the expected value of the area. For the normal unit vector {nx,ny,nz}, the expected projection area is:

The expected projection area

The surface area of an infinitesimal patch (θ,θ+ ⅆθ)×(φ,φ+ⅆφ) of the unit sphere is well known to be sin(θ) ⅆθφ. Dividing the infinitesimal area by the total surface area of the unit sphere, I obtain the infinitesimal probability measure, corresponding to the uniform distribution on the unit sphere: ⅆ S(θ,φ)==1/(4 π)sin(θ) ⅆθφ.

And, finally, the expected projection area equals:

Expectation[ 
 Total[Abs[{Sin[θ] Sin[ϕ], Sin[θ] Cos[ϕ], Cos[θ]}]], {θ, ϕ} ⎡ 
  ProbabilityDistribution[1/(4 π) Sin[θ], {θ, 0, Pi}, {ϕ, 0, 2 Pi}]]

3/2

Of course, I could simplify the computation, using the symmetry argument:

The symmetry argument

It is a well-known fact (also see this relevant question on math.stackexchange.com) that each individual component of a random vector, uniformly distributed on the unit sphere, follows a uniform distribution on the interval (-1,1). With this insight, the answer can be worked out in one’s head, explaining why this problem was deemed by Arnol’d a “trivium”:

3 Expectation[ Abs[nz], nz ⎡ UniformDistribution[{-1, 1}]]

3/2

Going to higher dimensions

The insight I gained allowed me to easily construct the answer for the case of the d-dimensional hypercube, projected on the randomly oriented hyperplane:

The case of the d-dimensional hypercube

I simply needed to know the distribution of a component of the d-dimensional random unit vector.

The computations are simple, and are based on hyperspherical coordinates (e.g. see Wikipedia: n-sphere). The infinitesimal hyperspherical area also factors as (Sin[θ1]d-2θ1)(Sin[θ2]d-3θ2)⋯(Sin[θd-2]ⅆθd-2)ⅆθd-1. Since the nd==cos(θ1), I needed to find the normalization constant for the quasi-density Sin[θ1]d-2:

Calculate the normalization constant for the quasi-density Sin[Subscript[θ, 1]]d-2

(Sqrt[π] Gamma[1/2 (-1 + d)])/Gamma[d/2]

Therefore the expected area of the projection of d-hypercube is:

The expected area of the projection of d-hypercube

(d Gamma[d/2])/(Sqrt[π] Gamma[(1 + d)/2])

Alternatively, I could use the closed-form result for the probability density function (PDF) of the distribution from here:

The closed-form result for the probability density function

True

I concluded this rewarding exploration by reproducing the results obtained earlier and tabulating results for higher dimensions:

Table[area[d], {d, 2, 10}]

{4/π, 3/2, 16/(3 π), 15/8, 32/(5 π), 35/16, 256/(35 π), 315/128, 512/(63 π)}

The average area of the shadow grows boundlessly with space dimension, as the number of faces that contribute to its area also increases:

DiscretePlot[area[d], {d, 2, 48}]

Plot of DiscretePlot[area[d], {d, 2, 48}]

In the high-dimensional space, the area scales as a square root of the dimension d:

Series[area[d], {d, Infinity, 4}]

Result of a square root of the dimension d

This concluded my exploration of Arnol’d's trivium problem with Mathematica. The use of Mathematica led to many “Aha!” moments, and enabled me to readily answer various “What if…” questions. It is my hope that I was able to convey the excitement of the discovery process largely accelerated with the use of Mathematica, and to inspire you to begin explorations of your own.

Download this post as a Computable Document Format (CDF) file.

Posted in: Mathematics
Leave a Comment

8 Comments


Carlo

Great post! It would be interesting to investigate the relationship between this problem and Buffon’s needle.

Posted by Carlo    June 3, 2013 at 11:23 am
Andrew Post

The answer is just one fourth of the surface of the cube. That’s a nice result of geometric probability, and a generalisation of Buffon’s needle experiment, which works for any convex body, not only the cube. So the solution in higher dimension is juste as simple, although the constant

In even dimension n = 2k, the ratio expected value of the shadow / surface of the body is:
[; \frac{2**{2k}}{2k \pi C(2k, k)} ;]
in odd dimension n = 2k+1, the ratio expected value of the shadow / surface of the body is:
[; \frac{C(2k, k)}{2**{2k+1}} ;]
Both expressions are equivalent to 1/ \sqrt{2 \pi n} for large n.

Posted by Andrew Post    June 3, 2013 at 4:22 pm
Sam Blake

Very nice, Sasha!

Posted by Sam Blake    June 4, 2013 at 3:02 am
Aidong Chen

Can you help to verify if we change the cube to a dodecahedron with side as 1, the expectation area is 3/(2 Tan[Pi/5])?

Posted by Aidong Chen    June 11, 2013 at 3:51 pm
    Aidong Chen

    Sorry, the result I got is 15/(2 Tan[Pi/5]). Aidong

    Posted by Aidong Chen    June 11, 2013 at 3:54 pm
      Oleksandr Pavlyk

      Per geometric probability theory the expected area of a shadow of any 3d convex body equals 1/4 of its surface area. Therefore the area of the shadow is 15/4 * Cot[Pi/5]. Try Simplify[PolyhedronData["Dodecahedron", "SurfaceArea"]/
      4 == (15/(4 Tan[Pi/5]))] in Mathematica.

      Posted by Oleksandr Pavlyk    June 12, 2013 at 7:56 pm
topdiablo3

Can you help to verify if we change the cube to a dodecahedron with side as 1, the expectation area is 3/(2 Tan[Pi/5])?

Posted by topdiablo3    June 26, 2013 at 3:09 am


Leave a comment

Loading...

Or continue as a guest (your comment will be held for moderation):